Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 45, 4, pp. 801-817, Warsaw 2007 PROPAGATION OF SPHERICAL SHOCK WAVES IN A DUSTY GAS WITH RADIATION HEAT-FLUX K.K. Singh North Eastern Hill University, Department of Mathematics, Shillong, India J.P. Vishwakarma D.D.U. Gorakhpur University, Department of Mathematics and Statistics, Gorakhpur, India e-mail : jpv univgkp@yahoo.com The propagation of spherical shock waves in a dusty gas with radiation heat-flux and exponentially varying density is investigated in the paper. The equilibrium flow conditions are assumed to be maintained, and the radiation is considered tobe of adiffusion type for anoptically thick grey gas model. The shock wave moves with variable velocity and the total energy of the wave is non-constant. Non-similar solutions are obtained, and the effects of variation of the radiation parameter and time are investigated. The effects of an increase in (i) the mass concentration of solid particles in the mixture and (ii) of the ratio of the density of solid particles to the initial density of gas on the flow variables in the region behind the shock are also investigated. Key words: shock waves, dusty gas, radiation heat-flux 1. Introduction Grover and Hardy (1966), Hayes (1968), Ray and Bhowmick (1974), Laum- bach and Probstein (1969), Verma andVishwakarma (1980) andmany others have discussed the propagation of shock waves in a medium where density varies exponentially. These authors have not taken radiation effects into ac- count. Laumbach and Probstein (1970a,b), Bhowmick (1981) and Singh and Srivastava (1982) obtained similarity or non-similarity solutions for the shock propagation in an exponentialmediumwith radiation heat transfer effects. Pai et al. (1980) obtained similarity solutions for a strong shockwave propagation 802 K.K. Singh, J.P. Vishwakarma in a dusty gas with constant density. Vishwakarma and Nath (2006) found similarity solutions for an unsteady flow behind a strong exponential shock driven out by a piston in a dusty gas in both cases, when the flow between the shock and the piston was isothermal or adiabatic. Vishwakarma (2000) stu- died the propagation of shockwaves in a dusty gas with exponentially varying density, using a non-similarity method. In the present work, we generalize the solution given by Vishwakarma (2000) for a strong explosion in a dusty gas (mixture of a perfect gas and small solid particles) taking radiation flux into account. It is assumed that the dusty gas is grey and opaque, and the shock is isothermal. The assumption that the shock is isothermal is a result of themathematical approximation inwhich the radiation flux is taken to be proportional to the temperature gradient, which excludes the possibility of a temperature jump (see, e.g. Zel’dovich andRaizer (1967), Bhowmick (1981), Singh and Srivastava (1982)). Radiation pressure and radiation energyare considered tobevery small in comparision tomaterial pressure and energy, respectively, and therefore only the radiation flux is taken into account. In order to get some essential features of shock propagation, small solid particles are considered as a pseudo-fluid, and it is assumed that the equilibrium flow condition is maintained in the flow field, and that the viscous stress and heat conduction of the mixture are negligible (Suzuki et al. (1976), Pai et al. (1980)). Although density of the mixture is assumed to be increasing exponentially, the volume occupied by the solid particles may be very small under ordinary conditions owing to large density of the particle material. Hence, for simplicity, the initial volume fraction of solid particlesZ1 is assumed to be a small constant (Vishwakarma (2000)). The effects of variation of the radiation parameter at different times on the flow variables behind the shock are obtained. The effects of an increase in (i) themass concentration of solid particles in themixture and (ii) the ratio of the density of solid particles to the initial density of gas on the flow variables behind the shock are also investigated. 2. Fundamental equations and boundary conditions The fundamental equations for one dimensional, spherically symmetric and unsteadyflowof amixture of gas and small solid particles taking radiation flux into account, can bewritten as (c.f. Vishwakarma, 2000; Singh andSrivastava, 1982) Propagation of spherical shock waves... 803 ∂u ∂t +u ∂u ∂r + 1 ρ ∂p ∂r =0 ∂ρ ∂t +u ∂ρ ∂r +ρ ∂u ∂r + 2ρu r =0 (2.1) ∂Um ∂t +u ∂Um ∂r − p ρ2 (∂ρ ∂t +u ∂ρ ∂r ) + 1 ρr2 ∂ ∂r (Fr2)= 0 where ρ is density of themixture, u – flowvelocity in the radial direction, p – pressure of the mixture, Um – internal energy per unit mass of the mixture, F – radiation heat flux, r – distance, and t – time. Assuming local thermodynamic equilibrium, and taking Rosseland’s diffu- sion approximation, we have F =− cµ 3 ∂ ∂r (aT4) (2.2) where ac/4 is the Stefan-Boltzmann constant; c – velocity of light; and µ – mean free path of radiation, which is a function of density ρ and absolute temperature T . Following Wang (1966), we have µ=µ0ρ α∗Tβ ∗ (2.3) where α∗ and β∗ are constants. The equation of state of a mixture of gas and small solid particles can be written as (Pai, 1977) p= (1−Kp) 1−Z ρR∗T (2.4) where R∗ is the gas constant, Z – volume fraction of solid particles in the mixture and Kp – mass concentration of solid particles. The relation between Kp and Z is given by Kp = Zρsp ρ (2.5) where ρsp is the species density of solid particles. In the equilibrium flow, Kp is a constant in the whole flow field. The internal energy of the mixture may be written as follows Um = [KpCsp+(1−Kp)Cv]T =CvmT (2.6) 804 K.K. Singh, J.P. Vishwakarma where Csp is the specific heat of the solid particles, Cv – specific heat of the gas at constant volume, and Cvm – specific heat of the mixture at constant volume. The specific heat of the mixture at constant pressure process is Cpm =KpCsp+(1−Kp)Cp (2.7) where Cp is the specific heat of the gas at the constant pressure process. The ratio of the specific heats of the mixture is given by (Marble, 1970; Pai, 1977) Γ = Cpm Cvm = γ ( 1+ δβ′ γ ) 1+ δβ′ (2.8) where γ = Cp Cv δ= Kp 1−Kp β′ = Csp Cv The internal energy Um is therefore, given by Um = p(1−Z) ρ(Γ −1) (2.9) We consider that a spherical shock wave is propagated into the medium, at rest, with small constant counter pressure. Also, the initial density of the medium is assumed to obey the exponential law ρ=Keαr (2.10) where α and K are suitable constants. The shock is assumed to be isothermal (formation of the isothermal shock is a result of themathematical approximation in which the flux is taken to be proportional to the temperature gradient. This excludes the possibility of a temperature jump, see for example Zel’dovich and Raizer (1967), Bhowmick (1981), Singh and Srivastava (1982)) and, hence, the conditions across it are ρ2(V −u2)= ρ1V =ms(say) p2+ρ2(V −u2)2 = p1+ρ1V 2 (2.11) Um2 + p2 ρ2 + 1 2 (V −u2)2− F2 ms =Um1 + p1 ρ1 + 1 2 V 2 Z2 ρ2 = Z1 ρ1 T2 =T1 Propagation of spherical shock waves... 805 where V = dR/dt denotes the velocity of the shock at r = R(t), indices 1 and 2 refer to the values just ahead and just behind the shock surface, and F1 =0 (Laumbach and Probstein, 1970). From equations (2.11), we get u2 =(1−β)V ρ2 = ρ1 β p2 =(1−Z1)ρ1V 2 Z2 = Z1 β (2.12) F2 =(1−β) [(1+Γ)β+(1−Γ)−2Z1 2(Γ −1) − 1−Z1 (Γ −1)M2e ] ρ1V 3 where β is given by β=Z1+ 1−Z1 ΓM2e (2.13) and M2e = V 2 a21 a21 = Γp1 ρ1(1−Z1) Me being the shock-Mach number referred to the speed of sound a1 in the dusty gas. The initial volume fraction of the solid particles Z1 is, in general, not constant. But the volume occupied by the solid particles is very small because density of the solid particles is much larger than that of the gas (Miura and Glass, 1985), hence Z1 may be assumed as a small constant. The expression for Z1 is (Pai, 1977; Naidu et al., 1985) Z1 = Kp G(1−Kp)+Kp (2.14) where G = ρsp/ρg is the ratio of the density of solid particles to the initial density of the gas. Values of Z1 for some typical values of Kp and G are given in Table 1. Table 1.Values of Z1 for some typical values of Kp and G Kp G Z1 0.2 10 0.02439 50 0.00498 100 0.00249 0.4 10 0.06250 50 0.01316 100 0.00662 806 K.K. Singh, J.P. Vishwakarma Let the solution to equations (2.1) and (2.2) be of the form (Ray and Bhowmick, 1974; Verma andVishwakarma, 1976; Singh and Srivastava, 1982; Vishwakarma, 2000) u= t−1U(η) ρ= tΩD(η) p= tΩ−2P(η) F = tΩ−3Q(η) (2.15) where η= teλr λ 6=0 (2.16) and the constants Ω and λ are to be determined subsequently.We choose the shock surface to be given by η0 = const (2.17) so that its velocity is given by V =− 1 λt (2.18) which represents the outgoing shock surface, if λ< 0. The solution toe equations (2.1) and (2.2) in form (2.15) are compatible with the shock conditions if Ω=2 λ=− α 2 α∗ =1 β∗ =− 5 2 (2.19) Since necessarily λ < 0, relation (2.19) shows that α > 0, meaning there- by that the shock expands outwardly in an exponentially increasing medium (Hayes, 1968; Vishwakarma, 2000; Yousaf, 1987). The strength of the shock, under these conditions, remains constant, for M2e = V 2 a21 = V 2 Γp1 ρ1(1−Z1) = 4(1−Z1)K Γp1α2η 2 0 = const From equations (2.18) and (2.19), we obtain R= 2 α log t τ (2.20) where τ is the duration of the initial impulse. Propagation of spherical shock waves... 807 3. Solution to the equations The flow variables in the flow-field behind the shock front will be obtained by solving equations (2.1) and (2.2). From equations (2.15), (2.18) and (2.19), we obtain ∂u ∂t =uλV −V ∂u ∂r ∂ρ ∂t =Vρα−V ∂ρ ∂r (3.1) ∂p ∂t =−V ∂p ∂r Using equations (3.1)-(3.3) and the transformations r′ = r R u′ = u V p′ = p p2 ρ′ = ρ ρ2 F ′ = F F2 (3.2) in fundamental equations (2.1) and (2.2), we obtain dρ′ dr′ = ρ′ 1−u′ [ 2log t τ + du′ dr′ + 2u′ r′ ] dp′ dr′ = ρ′ (1−Z1)β [ (1−u′) du′ dr′ +u′ log t τ ] dF ′ dr′ = (1−Z1) (1−β) [ (1+Γ)β+(1−Γ)−2Z1 2 − 1−Z1 M2 e ] · · {[(β−Z1ρ′)(1−u′)2ρ′ (1−Z1)β2 −Γp′ ]du′ dr′ + (3.3) + (1−u′)(β−Z1ρ′)ρ′u′ log tτ (1−Z1)β2 − 2Γp′u′ r′ } − 2F ′ r′ du′ dr′ = 1 p′β2(1−Z1)− (β−Z1ρ′)ρ′(1−u′)2 · · [F ′(1−Z1)β(1−u′) √ ρ′ log t τ NL √ p′ √ β−Z1ρ′ +(1−u′)(β−Z1ρ′)ρ′u′ log t τ + −2p′β2(1−Z1)log t τ − 2p′u′β2(1−Z1) r′ ] where N = 4acµ0α 3 √ R∗3 (3.4) 808 K.K. Singh, J.P. Vishwakarma is a non-dimensional radiation parameter and L= (Γ −1) √ (1−Z1)3 2(1−β) [ (1+Γ)β+(1−Γ)−2Z1 2 − 1−Z1 M2 e ] β √ (1−Kp)3 (3.5) Also, the total energy of the flow field behind the shock front is given by E = 16πKR3 βα2τ2 1 ∫ 0 [p′(1−Z1)(β−Z1ρ′) Γ −1 + 1 2 ρ′u′2 ] r′2 dr′ (3.6) Hence, the total energy of the shock wave is non-constant and varies with R3. In terms of the dimensionless variables r′, p′, ρ′, u′ and F ′, the shock conditions take the form r′ =1 p′ =1 ρ′ =1 u′ =1−β F ′ =1 (3.7) Equations (3.3) alongwith the boundary conditions (3.7) give the solution to our problem. The solution so obtained is a non-similar one, since motion behind the shock can be determined only when a definite value for time is prescribed. 4. Results and discussion The distribution of the flow variables behind the shock front is obtained by numerical integration of equations (3.3) with boundary conditions (3.7). For the purpose of numerical integration, the values of the parameters are taken to be (Pai et al., 1980; Miura andGlass, 1985; Vishwakarma, 2000; Singh and Srivastava, 1982), γ = 1.4; Kp = 0, 0.2, 0.4; G= 10, 50; β ′ = 1; M2e = 20; N =0.6, 0.8, 10; and t/τ =2, 4. Starting from the shock front, the numerical integration is carried out until the singularity of the solution p′β2(1−Z1)− (β−Z1ρ′)ρ′(1−u′)2 =0 is reached. This marks the inner boundary of the disturbance, and at this surface the value of r′ remains constant. The inner boundary is the position in the flow-field behind the shock front at which the Chapman-Jouget condition Propagation of spherical shock waves... 809 Fig. 1. Variation of reduced density ρ′ in the region behind the shock front for Kp =0.2 and G=50 Fig. 2. Variation of reduced pressure p′ in the region behind the shock front for Kp =0.2 and G=50 is satisfied, i.e., the position at which the line r/R(t)= const coincides with an isothermal characteristic. Figures 1 and 5, 2 and 6, 3 and 7, and 4 and 8 show variation of the reduced flow variables ρ′, p′,F ′ andu′, respectively, with reduced distance r′. 810 K.K. Singh, J.P. Vishwakarma Figures 1, 2, 5 and 6 show that, as we move inwards from the shock front, reduced density and pressure decrease, while Figures 3, 4, 7 and 8 show that the reduced radiation heat flux and fluid velocity increase. Fig. 3. Variation of reduced radiation heat flux F ′ in the region behind the shock front for Kp =0.2 and G=50 Fig. 4. Variation of reduced flow velocity u′ in the region behind the shock front for Kp =0.2 and G=50 Propagation of spherical shock waves... 811 Tables 2, 3 and 4 display the density ratio 1/β across the shock and the position of the inner boundary surface r′p (say) for various values of constant parameters. Table 2. Density ratio ρ2/ρ1 = 1/β across the shock front for different values of Kp and G Kp G ρ2 ρ1 = 1 β 0 27.99998 0.2 10 16.30122 50 23.43556 100 24.82961 0.4 10 9.96989 50 18.88506 100 21.42447 Table 3.Position of the inner boundary surface for different values of the radiation parameter N and time t/τ with Kp =0.2 and G=50 t τ N Position of the inner boundary surface (r′p) 2 0.6 0.95844 0.8 0.95828 10 0.95769 4 0.6 0.97232 0.8 0.97212 10 0.97127 It is found, fromFig.1 to Fig.4 andTable 3 that the effects of an increase in the value of the radiation parameter N, which depends on the mean free path of radiation, are: • to decrease the density ρ′, and to increase the pressure p′ at anypoint in the flow-field behind the shock. The decrease of density and the increase of pressure become significant near the inner boundary surface, • to increase the radiation heat flux F ′ and the velocity u′ near the inner boundary surface, and • to increase, slightly, the distance of the inner boundary surface from the shock surface. 812 K.K. Singh, J.P. Vishwakarma The effects of an increase in the time (t/τ) are (see Figs1-4): • to decrease the density ρ′ and the pressure p′, • to increase the radiation heat flux F ′ and the velocity u′, and • to decrease the distance of the inner surface from the shock front (see Table 3). Figures 5 to 8 show that for given values of N and G, the effects of an increase in themass concentration of the solid particles Kp at a given instant are • to increase the density ρ′, the pressure p′, the radiation heat flux F ′ and to decrease the flow velocity u′, • to decrease the slopes of the density, pressure, flow velocity profiles and to increase the slopeof the radiationheatfluxprofile in the regionbehind the shock front, and • to increase the distance between the inner contact surface and the shock front (see Table 4). This means that an increase in the mass concentra- tion of the solid particles has an effect to decrease the shock strength. Table 4.Position of the inner boundary surface for different values of Kp and Gwith N =10 and t/τ =2 Kp G Position of the inner boundary surface (r′p) 0 0.96241 0.2 10 0.94595 50 0.95769 100 0.95923 0.4 10 0.92257 50 0.95101 100 0.95494 Also Figs 5 to 8 show that for given values of N and Kp, the effects of an increase in the ratio of the density of the solid particles to the initial density of the gas G at a given instant are • to decrease the density ρ′, the pressure p′, the radiation heat flux F ′ and to increase the flow velocity u′, • to increase the slopes of thedensity, pressure,flowvelocity profiles and to decrease the slope of the radiation heat flux profile in the region behind the shock front, and Propagation of spherical shock waves... 813 • to decrease the distance between the inner contact surface and the shock front (seeTable 4).Thismeans that an increase in the ratio of thedensity of the solid particles to the initial density of gas has an effect to increase the shock strength. Fig. 5. Variation of reduced density ρ′ in the region behind the shock front for N =10 and t/τ =2 Fig. 6. Variation of reduced pressure p′ in the region behind the shock front for N =10 and t/τ =2 814 K.K. Singh, J.P. Vishwakarma Fig. 7. Variation of reduced radiation heat flux F ′ in the region behind the shock front for N =10 and t/τ =2 Fig. 8. Variation of reduced flow velocity u′ in the region behind the shock front for N =10 and t/τ =2 The effects of an increase in Kp or G on the shock strengthmay be expla- ined with the help of the compressibility of the medium as follows. The adiabatic compressibility of the mixture of the gas and small solid particles may be calculated as (c.f. Moelwyn-Hughes, 1961) Propagation of spherical shock waves... 815 C = 1 ρ (∂ρ ∂p ) S = 1−Z Γp where (∂ρ/∂p)S denotes the derivative of ρ with respect to p at a constant entropy. The non-dimensional compressibility C′ =C/C2 can be expressed as C′ = ( 1− Z1 β ρ′ ) p′ ( 1− Z1 β ) Fig. 9. Variation of non-dimensional compressibility C′ in the region behind the shock front for N =10 and t/τ =2 It is plotted against r′ in Fig.9. This figure shows that the compressibility decreases as the value of Kp increases, whereas it increases as the value of G increases. The decrease in the compressibility causes weaker compression of the gas behind the shock and, hence, a decrease in the shock strength. The increase in the compressibility causes stronger compression of the gas behind the shock and, hence, an increase in the shock strength. References 1. Bhowmick J.B., 1981,An exact analytical solution in radiation gas dynamics, Astrophys. Space Sci., 74, 2, 481-485 816 K.K. Singh, J.P. Vishwakarma 2. Grover R., Hardy J.W., 1966, The propagation of shocks in exponentially decreasing atmospheres,Astrophys. J., 143, 48-60 3. Hayes W.D., 1968, A self-similar strong shocks in an exponential medium, J. Fluid Mech., 32, 305-315 4. Laumbach D.D., Probstein R.F., 1969, A point explosion in a cold expo- nential atmosphere Part I, J. Fluid Mech., 35, 53-75 5. Laumbach D.D., Probstein R.F., 1970a, A point explosion in a cold expo- nential atmosphere Part II, Radiating flow, J. Fluid Mech., 40, 833-858 6. Laumbach D.D., Probstein R.F., 1970b, Self-similar strong shocks with radiation in a decreasing exponential atmosphere,Phys. Fluids, 13, 1178-1183 7. Marble F.E., 1970,Dynamics of dusty gases,A.Rev. FluidMech., 2, 397-446 8. Miura H., Glass I.I., 1985, Development of the flow induced by a piston moving impulsively in a dusty gas,Proc. Roy. Soc., London A, 397, 295-309 9. Moelwyn-Hughes E.A., 1961,Physical Chemistry, PergamonPress, London 10. Naidu G.N., Venkatanandam K., Ranga Rao M.P., 1985, Approximate analytical solutions for self-similar flows of a dusty gas with variable energy, Int. J. Eng. Sci., 23, 39-49 11. Pai S.I., 1977,Two phase flows, Chap. V, Vieweg Tracts in PureAppl. Phys., 3, Vieweg Verlag, Braunschweig, Germany 12. Pai S.I., Menon S., Fan Z.Q., 1980, Similarity solution of a strong shock wave propagation in a mixture of a gas and dusty particles, Int. J. Eng. Sci., 18, 1365-1373 13. Ray G.D., Bhowmick J.B., 1974, Propagation of cylindrical and spherical explosion waves in an exponential medium,Def. Sci. J., 24, 9-14 14. Singh J.B., Srivastava S.K., 1982, Propagation of spherical shock waves in an exponential medium with radiation heat flux,Astrophys. Space Sci., 88, 2, 277-282 15. Suzuki T., Ohyagi S., Higashino F., Takano A., 1976, The propagation of reacting blast waves through inert particle clouds, Acta Astronautica, 3, 517-529 16. VermaB.G., Vishwakarma J.P., 1976,Propagation ofmagnetogasdynamic plane shockwaves in an exponentialmedium,NuovoCimento,32B, 2, 267-272 17. Verma B.G., Vishwakarma J.P., 1980,Axially symmetric explosion inma- gnetogasdynamics,Astrophys. Space Sci., 69, 177-188 18. Vishwakarma J.P., 2000, Propagation of shock waves in a dusty gas with exponentially varying density,Eur. Phys. J. B, 16, 369-372 Propagation of spherical shock waves... 817 19. Vishwakarma J.P., Nath G., 2006, Similarity solutions for unsteady flow behind an exponential shock in a dusty gas,Phys. Scr., 74, 493-498 20. WangK.C., 1966,Approximate solutionofaplane radiating”pistonproblem”, Phys. Fluids, 9, 10, 1922-1928 21. Yousaf M., 1987, Strong shocks in an exponential atmosphere,Phys. Fluids, 30, 12, 3669-3672 22. Zel’dovich Ya.B., Raizer Yu.P., 1967, Physics of Shock Waves and High Temperature Hydrodynamics Phenomena, Vol. II, Academic Press, NewYork Propagacja sferycznych fal uderzeniowych w zanieczyszczonym gazie z uwzględnieniem radiacyjnej wymiany ciepła Streszczenie W pracy zajęto się problemem propagacji sferycznych fal uderzeniowych o wy- kładniczym rozkładzie gęstości w zanieczyszczonym gazie z uwzględnieniem radia- cyjnej wymiany ciepła. Założono równowagowewarunki przepływu czynnika, a samą radiację przyjęto typu dyfuzyjnegowmodelu optycznie nieprzezroczystego gazu. Fala uderzeniowa przemieszcza się ze zmienną prędkością, a całkowita energia fali również się zmienia. W analizie otrzymano rozwiązania niepodobne przy rozważaniu wpływ czasu i zmienności strumienia radiacji. Zbadano ponadto efekt wzrostu koncentracji masy cząstek stałych zanieczyszczenia oraz stosunku gęstości tych cząstek do począt- kowej gęstości gazu na parametry przepływuw obszarze bezpośrednio za czołem fali. Manuscript received March 9, 2007; accepted for print May 23, 2007