Abstract
For evaluating the motion of a solid body in a gaseous medium, one has to know the drag constant of the body. It is therefore not surprising that this subject was extensively investigated in the past. While accurate knowledge is available for the drag coefficient of a sphere in a steady flow condition, the case where the flow is time dependent is still under investigation. In the present work the drag coefficient of a sphere placed in a shock tube is evaluated numerically. For checking the validity of the used flow model and its numerical solution, the present numerical results are compared with available experimental findings. The good agreement between present simulations and experimental findings allows usage of the present scheme in nonstationary flows.
1. Theoretical Background
When a solid particle is exposed to a postshock gas flow, its response depends on the relative velocity that exists between the particle and the flow. Until the particle reaches a steady postshock flow velocity, the relative velocity between the particle and the gas flow changes and the particle motion is nonstationary. In shock tube experiments, the particle trajectory could be recorded accurately and its drag coefficient is evaluated from the recorded trajectory as follows. The equation of motion of a solid particle accelerated by the gas flow readswhere , , and are the solid sphere velocity, diameter, and material density, respectively. and are the gas velocity and density, respectively. It was shown in Igra and Takayama [1] that based on (1) the particle drag coefficient () and the appropriate Reynolds number can be expressed as follows:where and are components of the velocity vector in and directions, respectively, is the gravity acceleration, and is the gas viscosity. Similarly, the sphere’s Mach number, based on the relative velocity, readswhere and are the gas specific heat ratio and the gas constant, respectively. In the following, first the experimental results of Tanno et al. [2] are simulated. In their experiments a 80 mm sphere was kept on a strut in a 150 × 500 mm cross section shock tube. The case of free moving sphere will be considered in a separate paper; for such cases (1)–(3) are relevant. When the sphere is kept stationary throughout its collision with the oncoming shock wave, the drag coefficient is directly deduced from the computed pressure distribution around the sphere. In simulating the shock tube experiments the planar incident shock wave is initially situated at some distance upstream of the sphere. The main parameters controlling the flow are the sphere’s diameter, the cross section of the shock tube where the sphere is placed, the mass of the sphere, and the strength of the incident shock wave. As mentioned, in the following computations the sphere drag coefficient is deduced from the computed, nondimensional pressure coefficient around the sphere.
2. Numerical Scheme
System of mass, momentum, and energy conservation equations for ideal nonviscous gas for a moving finite volume can be written in the following form:where , , is the vector of conserved variables, is vector of fluxes normal to the boundary of the control volume, , , , and is the velocity vector of the control volume boundary.
The explicit time step operator that represents a finite difference algorithm for approximating the system of (4) can be factored into a symmetric product of time step operators in directions. Vector of fluxes in a normal direction to the boundary of grid cell is determined on a basis of the finite difference scheme suggested by Yee et al. [3]. A slightly improved version of Harten’s scheme was suggested by Martyushov [4] and was used in the present calculations. These improvements are briefly outlined in the following. Vinokur [5] showed that Roe’s procedure employs the following formula for the sound velocity.
This expression may yield results outside the interval . In this case the procedure for calculating eigenvectors is incorrect and the sign of the eigenvalues can be wrong. In the present version of Harten scheme, this situation was taken into account and in the case when is outside the interval Roe interpolation for the values at boundary between cells is replaced by a half sum of the corresponding variables in the considered cells. For calculation of characteristic variables at the cell boundary using Harten’s scheme geometric and gas dynamics variables is used in five cells: . In order to reduce the influence of the geometric parameters of distant cells instead of the characteristic variables , , and in the present computation, pseudocharacteristic variables , , and are used. Detailed investigation of the algorithm and its validation for one-dimensional flows can be found in Il’in and Timofeev [6].
The numerical grid used in the present calculations was generated using the Thompson-type algorithm based on solving a system of three Poisson equations (see Thompson et al. [7] and Martyushov [8]). Example of the employed grid is shown in Figure 1 where every fifth coordinate line is plotted.
For calculation of the shock wave diffraction over the sphere placed inside a shock tube it is convenient to separate the flow domain into two subdomains: the flow domain upstream of the sphere and the flow domain downstream the sphere. The calculations in this case are performed using the block-structured grid consisting of two blocks.
For simultaneous calculation of the flow parameters in both subdomains it is necessary to determine formulas for coupling numerical solutions in both subdomains. This can be done in different ways. One method is to determine the boundary conditions which couple numerical solutions in both subdomains. For the difference scheme of Yee et al. [3] it is sufficient to calculate the flows at the boundary using formulas similar to the cells main scheme ones, where gas dynamics variables at the boundary between the two subdomains are calculated using Roe’s procedure using gasdynamic variables in the right and in the left regions.
3. Results and Discussion
In the past two basically different shock tube experiments were conducted aimed at studying the drag force (and drag coefficient) acting on a sphere as a result of its head-on collision with a planar shock wave. In one type of experiment the sphere was free to move as a result of its collision with the incident shock wave. Its speed depends on its mass and the shock intensity. Examples of publications discussing such type of shock tube experiment are the work of Igra and Takayama [1], Suzuki et al. [10], and Jourdan et al. [11]. In the second category the sphere is fixed inside the shock tube; it does not move throughout the investigated time. Examples for such experimental and numerical investigation are the work of Tanno et al. [2] and Sun et al. [9]. In the present study the pressure distribution around the surface of a fixed sphere which resulted from its head-on collision with the incident shock wave is calculated and compared with the above-mentioned experimental findings. The drag force acting on the sphere is deduced from the computed pressure distribution. In the case of a free sphere, the sphere starts moving after its collision with the incident shock wave. When simulating this case the grid moves with the sphere. The velocity of grid points at every time step is calculated using the following relation:where for system (4) . Once the sphere velocity is computed the sphere’s trajectory can easily be assessed and compared with experimental findings. This will be done in a future paper.
The following initial values were used in the present computations for a fixed sphere: sphere diameter 80 mm, incident shock wave Mach number 1.22, preshock pressure 101 kPa, and the initial temperature equal to the room temperature, that is, 300 K. The early stage of the shock wave diffraction over the sphere is shown in Figure 2. At this time, = 208 μs, after the incident shock hits the sphere one can see clearly the incident and the reflected shock waves.
From the calculated pressure fields behind these shock waves one can evaluate the drag force acting on the sphere. The obtained dimensionless pressure distributions, , around the sphere at different times during the shock diffraction are shown in Figure 3. Here , where and are the pressures at the stagnation point and the static pressure ahead of the incident shock wave, respectively. These numerical results (solid lines in Figure 3) are compared with the experimental findings of Tanno et al. [2], appearing as circles in Figure 3. These circles were copied from appropriate plot appearing in paper by Tanno et al. [2]. Tanno et al.’s experimental findings were recorded using pressure transducers placed at equidistant distribution along the sphere surface (15 degrees) during the following times: 140 μs, 208 μs, 296 μs, and 380 μs, measured from the time when the incident shock wave reached the sphere.
(a)
(b)
(c)
(d)
While a very good agreement is found at early times between the present numerical results and the appropriate experimental finding (see Figures 3(a) and 3(b) at a later time) a difference is noticed between the two (see Figures 3(c) and 3(d)) (the time, , appearing in Figure 3 is in microseconds). This discrepancy is caused by the strut which is used in the experiments for keeping the sphere fixed inside the test section. This strut is attached to the sphere at its trailing edge. The higher pressures measured in the experiments arise from the shock interaction with this strut. The computed pressure distribution around the sphere can be used for evaluating the sphere’s drag coefficient: , where is the drag force exerted on the sphere, is the radius of the sphere, and and are density and velocity behind the incident shock wave. The next results to be shown are simulations made to the experimental findings of Sun et al. [9] conducted also for stationary spheres of different sizes. The present numerical results (solid lines) and the appropriate experimental findings taken from Sun et al. [9] (open circles) are shown in Figure 4. For comparing our numerical results with findings by Sun et al. [9], the time is normalized by a coefficient , where /s. Results shown in Figure 4 were obtained for the following initial conditions: incident shock wave Mach number of , preshock pressure of 101 kPa, and preshock temperature of 300°K. The results shown in Figure 4 are for a sphere diameter of 80 mm. As is apparent from Figure 4, good agreement exists between the present computational results and the appropriate experimental findings of Sun et al. [9]. As could be expected, a very large drag coefficient is observed in Figure 4 during the shock diffraction over the sphere; that is, for normalized time < 1.6. A peak value of is reached during this time period. The value is significantly higher than that obtained in a similar shock-free flow.
Once the shock diffraction over the sphere is completed, a significant reduction in is evident in Figure 4. For Figure 4 suggests that . This negative drag coefficient is a direct result of the experimental setup, testing a large sphere (80 mm in diameter) in a relatively small shock tube. The shocks reflected form the shock tube walls hit first the sphere’s rear surface to result in a negative drag force. However, as is evident from Figure 4 with progressing time the drag coefficient increases toward the appropriate steady flow value.
4. Appendix
For checking convergence of the present solution, it was repeated using three different grid sizes; the used grids are shown in Figure 5. Results obtained for the sphere drag coefficient while using the three different grids are shown in Figure 6. It is clear from this figure that the solution quickly converged to a unique result.
(a)
(b)
(c)
5. Conclusion
The present study focused on simulating the complex flow field which resulted from the head-on collision of a planar shock wave with a fixed sphere. The good agreement obtained between the present numerical results and experimental findings attested to the validity of the simple physical model used (see (4) describing inviscid flow) and the efficiency of the used numerical scheme. It also prepared the foundation for simulating the case where the sphere is free to move following its interaction with the oncoming shock wave. In the case of a freely moving sphere the experimental flow duration is significantly longer and therefore viscous effects cannot be ignored any longer.
Competing Interests
The authors declare that they have no competing interests.