Shock waves in non-Newtonian bubbly liquids

International Journal of Multiphase Flow 27 (2001) 635±655

www.elsevier.com/locate/ijmul¯ow

Shock waves in non-Newtonian bubbly liquids q A.A. Gubaidullin *, O.Sh. Beregova, S.A. Bekishev Tyumen Institute of Mechanics of Multiphase Systems, Siberian Branch of the Russian Academy of Sciences, Taimyrskaya St. 74, Tyumen 625000, Russian Federation Received 5 September 1999; received in revised form 5 January 2000

Abstract Propagation of non-stationary shock waves in a bubbly liquid with a non-Newtonian carrier phase was investigated numerically. The research is carried out within the one-velocity two-temperature two-pressure model of multiphase media mechanics. An aqueous solution of polymer was used as a non-Newtonian carrier phase. It was concluded that the behavior of shock waves in Newtonian and non-Newtonian bubbly liquids with the same Newtonian viscosity may have fundamental di€erences. The in¯uence of de®ning parameters of the two-phase mixture and wave on the structure of shock waves was examined. The numerical analysis of evolution and attenuation of pulse perturbations was carried out. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Bubbly liquid; Non-Newtonian; Polymer solution; Shock wave; Viscosity; Relaxation time

1. Introduction The wave dynamics of bubbly systems with viscous and non-viscous Newtonian carrier phase is developed well enough for the present time (Gubaidullin et al., 1978; Nakoryakov et al., 1990; Gubaidullin, 1991; Nigmatulin, 1991). However, wave ¯ows of those systems based on nonNewtonian liquids with complicated rheology (polymeric solutions and melts, suspensions, highparanaceous, waxy and resinous, tar-content oils, etc.) have not been studied adequately in spite of their wide incidence in practice. Some problems of the behavior of a bubble in viscoelastic relaxing polymeric liquid and propagation of acoustic waves in the same liquid with bubbles are

q

This is an extended version of a paper that was ®rst presented at ICMF'98 (Lyon, France) and subsequently selected by its Scienti®c Committee. * Corresponding author. E-mail address: [email protected] (A.A. Gubaidullin). 0301-9322/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 3 0 1 - 9 3 2 2 ( 0 0 ) 0 0 0 3 0 - 6

636

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

researched in this domain. These problems were best described by Levitskiy and Shulman (1995). The problems of dynamics of non-stationary nonlinear waves are not being studied. This paper investigates the features of propagation of non-stationary shock waves in a nonNewtonian liquid with gas bubbles. Our intent is to compare the wave behavior in bubbly mixtures with non-Newtonian and viscous Newtonian carrier phases and to reveal the in¯uence of de®ning parameters of two-phase mixture (Newtonian viscosity of solution and density of liquid, void fraction and size of bubbles, kind of gas, relaxation time, etc.) on the wave evolution process.

2. Basic equations In this section we are going to obtain a system of equations, describing the dynamic behavior of the two-phase mixture mentioned above, taking into account the following assumptions (Gubaidullin et al., 1978): the radius of bubbles is many times larger than molecular-kinetic dimensions and many times smaller than the distances at which averaged or macroscopic parameters of mixture or phases vary essentially; the mixture is monodispersive; interaction, collisions of bubbles to one another may be neglected; processes of bubble fragmentation, coagulation and formation of new bubbles are absent; the pressure of gas in bubbles is homogeneous, velocities of macroscopic motion of phases are equal, the temperature and the density of carrier liquid are constant, phase transitions are absent. 2.1. The equations of motion of phases and of heat transfer The continuity equations for gas, liquid and unit bubble and momentum equation at the assumptions mentioned above for non-stationary one-dimensional ¯ow along axis x are as follows: oq1 oq1 v ˆ 0; ‡ ox ot dv op ˆ0 q ‡ dt ox

oq2 oq2 v ˆ 0; ‡ ox ot 

d 0 3
q a ˆ 0; dt 2

 d o o ˆ ‡v : dt ot ox

…1†

…2†

Here q; p; v are, respectively, the density, the pressure and the velocity of mixture, q0i ; qi ; ai ; pi are the true and reduced densities, the void fraction, and the pressure of the ith phase, the subscripts i ˆ 1 and 2 refer to liquid and gas parameters, a is the bubble radius. The heat in¯ux equation for gas may be written as q2

du2 dt

a2 p2 dq02 ˆ nq2 ; q02 dt

T1 ˆ const:;

…3†

where Ti and q2 are the temperature and heat in¯ux to a bubble from liquid and n is the number of bubbles per unit volume of mixture. Within the two-temperature scheme (Gubaidullin et al., 1978) q2 is given by

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

637

k2 Nu2 …T1 T2 †; 2a where k2 is the gas thermal conductivity and Nu2 is the Nusselt number, and it can be expressed (Nigmatulin, 1991) as ( 10; Pe2 6 100; Nu2 ˆ p Pe2 ; Pe2 > 100; q2 ˆ 4pa3

Pe2 ˆ 12…c2 …T†

m2 ˆ

k2 0 q2 cp2

1† ;

T1 jT1

ajwj T2 j m…T† 2

;

c2 ˆ cp2 =cv2 :

…T†

Here m2 is the gas thermal di€usivity, w the radial velocity of bubble interface, cp2 ; cv2 the speci®c heats of gas at constant pressure and volume and Pe2 is the Peclet number. The gas is assumed to be ideal and the liquid is incompressible p2 ˆ q02 …c2

1†cv2 T2 ;

u2 ˆ cv2 T2 ;

q01 ˆ const:

There are, by de®nition q1 ˆ a1 q01 ; q2 ˆ a2 q02 ; q ˆ q1 ‡ q2 ; a1 ‡ a2 ˆ 1; 4 a2 ˆ pa3 n; p ˆ a1 p1 ‡ a2 …p2 2R=a†; 3 where R is the surface tension of liquid. 2.2. Conditions of simultaneous strain of phases As a condition of simultaneous strain of phases, a Rayleigh±Lamb equation for radial oscillations of a bubble in liquid is commonly used. The one generalized to the case of non-Newtonian liquid was obtained by Levitskiy and Shulman (1995)). It has the following form in a spherical coordinate system …r; u; w†:   Z 1 …rr† dw 3 2 2R s s…uu† 0 ‡ w ‡ p1 p2 ‡ ˆ2 q1 a dr; dt 2 a r a …4† da wˆ ; dt where s…rr† s…uu† is the normal stress di€erence. The stress tensor s^ is de®ned by the rheological state equation of Maxwell type with the upper convective derivative and the single relaxation time ts (Astarita and Marrucci, 1974; Levitskiy and Shulman, 1995) s^ ˆ s^s ‡ s^p ; " d^ sp s^p ‡ ts dt

s^s ˆ 2gs e^; gs ‡ gp ˆ g0 ; #  1 s^p
e^ ‡ e^
s^p ˆ 2gp e^; 2

638

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

where e^ is the tensor of strain rate; the point means the internal multiplication of tensors of the second rank (multiplication of their matrices), g0 ; gs ; gp are Newtonian viscosity of solution (or zero-shear viscosity), the viscosities of solvent and of polymeric network in solution, respectively; the subscripts p and s refer to polymer and solvent. The equations for normal components s…rr† ; s…uu† are written in the form (Levitskiy and Shulman, 1995) ! 2 dsp…rr† a2 w …rr† …rr† a w ‡ 2sp ˆ 4g ; sp ‡ ts p r3 r3 dt …5† ! …uu† 2 2 ds aw aw p ˆ 2gp 3 ; ‡ ts s…uu† sp…uu† 3 p r r dt ˆ s…rr† s

4gs

a2 w ; r3

s…uu† ˆ 2gs s

a2 w : r3

The expressions for s…rr† ; s…uu† can be found from (5) and substituted into (4). They give   dw 3 2 2R 0 ˆ p2 p1 ‡ S; q1 a ‡ w ‡ dt 2 a

…6†

…7†

S ˆ Sp ‡ Ss ; 4gp Sp …t† ˆ 2 ts a …t † w Ss ˆ 4gs : a

Z 0

t

e

…n t†=ts

a…n†w…n† dn;

…8† …9†

2.3. The constitutive equation of liquid for radial ¯ow around a bubble The integral Sp can be transformed to the di€erential equation   gp w dSp 1 2w ‡ Sp ˆ 4 : ‡ ts a dt ts a

…10†

One can see from comparison of (10) and (9) with (5) and (6) that when r ˆ a; Sp ˆ …rr† s…rr† p ; Ss ˆ ss , i.e. the parameters Sp and Ss characterize the stresses on a bubble interface in polymer network and solvent, correspondingly (Levitskiy and Shulman, 1995). Let us analyze, ®rst, the limiting cases of Eq. (10). It will be assumed that t is the characteristic time of bubble oscillations, which determines the time scale of the process, Sp is the scale of stress can be oscillations in polymer. Timept  estimated, for instance, as a bubble collapse time in inertial Rayleigh regime t  t0 ˆ a0 q010 =p0 . This estimation is valid when the dissipation is less than the inertia of radial movement of liquid around a bubble, i.e. oscillations of a bubble are possible. The estimates of terms of Eq. (10) are obtained:

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

1 Sp dSp Sp  ; Sp  ; dt ts t ts   w S da p since a  a0 ; and w ˆ  a0 =t ; 2 Sp  a dt t 4 gp w  gp : t a t t s

639

…11†

s 

Small relaxation times. When ts  t the ®rst and third terms in the left-hand side of Eq. (10), which contain time derivatives, may be neglected in comparison with Sp =ts . On both the righthand side and the left-hand side of Eq. (10) there remains one term. These terms must be of the g same order, i.e. Sp  tp . The following is the approximate expression for stress in polymer w Sp  4gp ; ts  t ; …12† a where Eqs. (9) and (12) allow one to write the expression for stress in solution as w S ˆ Sp ‡ Ss  4g0 ; ts  t : a

…13†

Hence (12) and (13) have the same form as for a Newtonian liquid with the viscosity being equal to Newtonian viscosity of a polymer solution. Large relaxation times. In the opposite case when ts  t , it follows from (11) that one may ignore the second term Sp =ts in the left-hand side of (10) in comparison with the ®rst and third terms, which involve derivatives. Thus, the following equation is obtained as dSp w ‡ 2 Sp ˆ a dt

4

gp w ; ts a

ts  t ;

…14†

which taking into account the initial conditions …t ˆ 0 : Sp ˆ 0; a ˆ a0 ; † has the solution   gp  a0 2 Sp ˆ 2 1 ; ts  t ; …15† ts a It corresponds to the nonlinear-elastic returning force. When the bubble radius deviates to a small extent, (15) may be reduced to the following form: gp Da ; ts  t ; ts a0 a ˆ a0 ‡ Da; jDja  a0 :

Sp 

4

…16†

Note that e…rr† ˆ 2Da=a0 is the longitudinal strain in the liquid at bubble interface, so gp …rr† …rr† Sp ˆ sp  2 e : rˆa ts rˆa Expressions (15) and (16) show that in this limiting case, the polymeric network in solution behaves as on elastic medium, with the elastic modulus G ˆ gp =ts determining the scale of stress oscillations Sp  gp =ts . The total stress in solution consists of the elastic component due to polymer and the viscous one due to solvent:

640

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

 gp  a0 2 S ˆ Sp ‡ Ss  2 ts a



1

4gs

w ; a

ts  t :

…17†

E€ective gas elastic modulus. It is obvious that the elastic stress in polymer will in¯uence the behavior of a bubble if Sp exceeds the increment in gas pressure Dp2 in a bubble due to Da. Let us analyze this case. The simple polytropic law is as follows:
j p2 =p20 ˆ q02 =q020 ; where j is the polytropic index. Let p2 ˆ p20 ‡ Dp2 ;

a ˆ a0 ‡ Da;

where p20 ˆ p1 ;

jDp2 j  p20 ;

jDaj  a0 :

Since q02 =q020 ˆ a30 =a3 , then Dp2 

3jp20

3 Gg ˆ jp20 ; 4

Da ˆ a0

4Gg

Da ; a0

…18†

where Gg is the e€ective elastic modulus of the gas in a bubble. The elastic stress in polymer must be taken into account if Sp P jDp2 j; that is G P Gg : Note that the equilibrium value of pressure pe behind the front of the wave should be used instead of p20 in (18) if applied to a shock wave of ``step'' type. On the base of the performed analysis it may be argued that a polymeric liquid with bubbles for two limiting cases mentioned above behaves as a Newtonian bubbly liquid with viscosities g0 and gs , correspondingly, i.e. a shock wave might have two limiting structures at varying relaxation time: a shock wave traveling as in a Newtonian bubbly liquid with the viscosity being equal to Newtonian viscosity of solution g0 …ts  t † and a shock wave traveling as in a Newtonian bubbly liquid with the viscosity being equal to that of solvent gs …ts  t †. Modeling of polymeric solution by a Newtonian liquid is possible in these cases if the elastic modulus of polymeric network is less than the e€ective elastic modulus of the gas in a bubble.

3. Numerical results The system of Eqs. (1)±(3) and (7)±(10) described above is solved numerically. The structure of the system di€ers from that for Newtonian bubbly liquid in one additional di€erential equation (10) for the stress in polymer. The latter involves only time derivatives, so the procedure proposed by Gubaidullin et al. (1978) for numerical modeling of Newtonian bubbly liquids can be applied. Following this procedure the system after conversion to the Lagrange coordinates is transformed so that the resulting system consists of one second-order di€erential equation on a space variable only, whereas the others are the ®rst-order equations with respect to time. The modi®ed Euler± Cauchy method is applied to solve the equations on time and the sweep method to solve the

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

641

boundary value problem for the second-order equation at each step on time. A more detailed description of the method numerical modeling is represented in Appendix A. 3.1. Simulation of stepwise shock wave This research examines special features of those processes and the in¯uence of de®ning parameters of the medium and initial parameters of the wave. The wave of step type is initiated by a drastic increase in pressure at the boundary of the medium and then the evolution of compression wave is considered. Varying relaxation time. One of the most interesting simulation results is shown in Fig. 1, where the in¯uence of relaxation time ts on structure of shock waves is examined. The pro®les of nondimensional pressure at x ˆ 20 cm for various values of relaxation time ts ˆ 0:001; 0:1 and 1 ms are presented. Parameters of the mixture and the wave are q010 ˆ 998 kg=m3 ; g0 ˆ 1:48 Pa s; gs ˆ 1 mPa s; a0 ˆ 1 mm; a20 ˆ 2%; T0 ˆ 293 K; p0 ˆ 0:1 MPa; pe ˆ 0:3 MPa. One can see from Fig. 1 that the amplitudes of oscillations signi®cantly increase when the relaxation

Fig. 1. Pro®les of pressure in polymer liquid with air bubbles in shock wave for the various values of relaxation time ts .

642

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

time increases. A certain increase of oscillation period is also observed. The calculations show that the pro®les of pressure at ts < 0:001 ms practically coincide with each other and with the corresponding curve for the Newtonian carrier liquid with viscosity being equal to Newtonian viscosity of solution. In the other limiting case, ts P 1 ms, the pro®les of pressure coincide as well with the curve for the bubbly liquid with constant viscosity being equal to the viscosity of solvent. It fully corresponds to the previous analysispof limiting cases of Eq. (10). Characteristic Rayleigh time for  the calculations discussed is t0 ˆ ae q010 =pe  0.04 ms. Characteristic time can be estimated also from the plots as reciprocal average cyclic frequency of oscillations t  1=x  0:03 ms, which is reasonably close to Rayleigh time. Time t separates the limiting pro®les in Fig. 1 on ts scale by nearly equal spaces in 1.5 decimal orders at each side of t (0.001, 0.03, 1 ms). Newtonian and non-Newtonian carrying ¯uids. The behavior of shock waves in bubbly liquids with viscous Newtonian and non-Newtonian carrier phases is compared. Such computer analysis is carried out by the examples of glycerin and aqueous solution of polymer, the average molecular weight and concentration of which can be selected so that Newtonian viscosity of solution is equal to the viscosity of glycerin. The de®ning parameters of polymer solution are q010 ˆ 998 kg=m3 ; R ˆ 0:073 kg=s2 ; g0 ˆ 1:48 Pa s; gs ˆ 1 mPa s and ts ˆ 1 ms. We assume that the density and the surface tension of solution are close to those of the solvent (Ishiguro and Hartnett, 1992). The parameters of glycerin are q010 ˆ 1260 kg=m3 ; gN ˆ 1:48 Pa s. The evolution of ``step'' type shock wave of intensity pe ˆ 0:15 MPa in these liquids with air bubbles (a0 ˆ 1 mm, a20 ˆ 2%; T0 ˆ 293 K ; p0 ˆ 0:1 MPa) is presented in Fig. 2. It can be seen that the behavior of

Fig. 2. Pro®les of pressure in step-type shock wave in polymer solution and in glycerin with air bubbles. (a) Polymer. (b) Glycerine.

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

643

shock wave in these two cases is fundamentally di€erent: the wave has a monotonous structure in the Newtonian bubbly liquid and oscillatory structure in the non-Newtonian one. As shown below, it occurs because in the non-Newtonian liquid, the e€ective viscosity becomes signi®cantly lower than Newtonian viscosity of solution due to the oscillations of bubbles. It should be noted that the velocity of the shock wave in the ®rst case is lower, because the density of glycerin is greater than that of the polymeric solution. 3.2. Simulation of nonlinear compression pulses The character of damping of nonlinear compression pulses is studied. The results of computational investigation are illustrated by the example of compression waves with the initial semisinusoidal form. The picture of damping the same initial waves (intensity pe ˆ 0:3 MPa, initial duration of pulse is 10 ls) in the non-Newtonian liquid (polymer solution, ts ˆ 0:1 ms, other parameters are as in Fig. 1) and in the Newtonian one (glycerin) with air bubbles is presented in Fig. 3. It may be noted that the character of damping varied, namely, in the non-Newtonian bubbly liquid it comes less intensively. It occurs due to the above-mentioned decrease of e€ective viscosity in the non-Newtonian liquid, which becomes much lower than its Newtonian viscosity, that is the viscosity of glycerin. Let us inspect now how the value of relaxation time in polymer can in¯uence the damping of pulses. We will consider the propagation of semi-sinusoidal pulses of di€erent duration but of the ®xed intensity 0.3 MPa in bubbly mixtures with various carrying phases: water, glycerin and aqueous polymer solutions (g0 ˆ 1:48 Pa s; gs ˆ 1 mPa s; q010 ˆ 998 kg=m3 ) with di€erent relaxation times. Gas is air, a0 ˆ 1 mm; a20 ˆ 0:5%; p0 ˆ 0:1 MPa; T0 ˆ 293 K. In Fig. 4 there are the integral damping curves of pulse disturbance of initial duration 1 ms in aqueous polymer solutions. One can see that with the relaxation time growing, the damping decreases. At the beginning of propagation, the pressure in the medium is higher than in initial pulse ± this ampli®cation is due to the e€ect of radial inertia. The analogous damping curve calculated for

Fig. 3. Damping of compression pulse in non-Newtonian and Newtonian liquids with air bubbles.

644

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

Fig. 4. The integral damping curves for semi-sinusoidal pulse in aqueous polymer solution with air bubbles for various relaxation times ts .

air±water mixture fully coincides with the dotted line (ts ˆ 1 ms) and that one for air±glycerin mixture does this with the dashed line (ts ˆ 0:01 ms). Thus, for pulse disturbances as for the stepwise shock waves, the conclusion remains true that in the limiting cases, when ts  t or ts  t , the polymeric carrying phase behaves as a Newtonian liquid with the viscosities gs and g0 , respectively. This is true for the pulse disturbances independently on their initial duration. Then we studied the in¯uence of de®ning parameters of mixture on the pulse shock wave. Fig. 5 shows pressure pro®les of the same damping pulse as in Fig. 4 but in di€erent medium (g0 ˆ 1:48 Pa s; gs ˆ 1 mPa s; ts ˆ 0:1 ms) with air bubbles (a0 ˆ 0:1 mm, a20 ˆ 1%; p0 ˆ 0:1

Fig. 5. Damping of semi-sinusoidal pulses in polymer solutions with air bubbles with di€erent densities of carrying phase.

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

645

MPa) at two values of carrying phase density q010 ˆ 998 kg=m3 as in water (dotted lines) and q010 ˆ 13546 kg=m3 as in mercury (solid lines). It can be seen that the damping increases with the increase of density, which can be explained by the growth of thermal dissipation. In addition, the form of propagating pulses changes: they have a look of oscillating wave in one but a soliton look with an oscillating tail in the other mixture with a heavier liquid. Further calculations have shown that the in¯uence of de®ning parameters of mixture and initial pulse (bubble radius, density and viscosity of carrying phase, duration of pulse) on the wave propagation in a polymeric bubbly liquid is qualitatively the same as in a Newtonian one. Non-Newtonian properties show themselves when the values of relaxation time are close to the characteristic time of bubble oscillations.

4. Analysis and discussion 4.1. Analysis of linear bubble oscillations As it is known (Nigmatulin, 1991), the oscillatory pro®les in the shock waves propagating through bubbly liquids appear due to bubble oscillations. To understand better the processes which might occur, a more thorough harmonic analysis of dynamic characteristics of a nonNewtonian liquid around an oscillating bubble will be carried out. Expressions for complex amplitudes can be obtained from Eqs. (9) and (10) in case of harmonic oscillations of bubbles with small real amplitude Da > 0 a ˆ a0 ‡ Da
exp…ixt†;

Da  a0 :

…19†

The stresses Sp ; Ss will be found as Sk ˆ Sk exp…ixt†;

k ˆ p; s:

…20†

The velocity of the bubble interface w may be written as W ˆ ixDa exp…ixt†; i.e. the amplitude is w ˆ ixDa;

…21†

then w ixDa
exp…iwt† Da exp…iwt†;  ix ˆ a a0 ‡ Da
exp…iwt† a0 i.e.  w  a

 ix

Da : a0

…22†

By substituting (19)±(22) into (10) and (9), one can obtain the expressions for the amplitude of stress in polymer Sp and in solvent Ss

646

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

Sp ˆ Ss ˆ

Da ix ; a0 1 ‡ ixts Da 4gs ix : a0 4gp

…23† …24†

The last expression is a purely imaginary number that re¯ects the viscous response of the Newtonian solvent. The amplitude of the total stress is   gp Da    ‡ gs : …25† S ˆ Sp ‡ Ss ˆ 4ix a0 1 ‡ ixts The real and imaginary parts of Sp ˆ Sp0 ‡ iSp00 have the following form: Sp0 ˆ Sp00 ˆ

gp Da …ts x†2 ; ts a0 1 ‡ …ts x†2 xDa 1 4gp ; a0 1 ‡ …ts x†2

4

…26† …27†

which describe the elastic and viscous responses of the liquid to the strain, respectively. Their ratio is Sp0 ˆ ts x  ts =t ; Sp00 i.e. when ts  t …ts x  1† the imaginary part predominates xDa  ; Sp  Sp00  4gp a0 and when ts  t …ts x  1†, the real part predominates gp Da  : Sp  Sp0  4 ts a0

…28†

…29†

…30†

Note that (29) and (30) are agreed with viscous and elastic limits (12) and (16). 4.2. Complex viscosities and moduli The coecient of the complex dynamic viscosity at small oscillations of a bubble in nonNewtonian liquid can be introduced in the same manner as for small shear oscillations (Ferry, 1980). In the Rayleigh±Lamb equation, the viscosity of a Newtonian liquid gN is provided by the term SN ˆ

4gN w=a:

And the complex amplitude is SN ˆ

4gN ix

Da ; a0

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

647

thus SN : 4ixDa=a0

gN ˆ

We may introduce the complex dynamic viscosity in a way similar to that in the last expression g ˆ g0

ig00 ˆ

S : 4ixDa=a0

…31†

Here g0 is the real component of the complex viscosity corresponding to the viscous response and g00 is the imaginary component corresponding to the elastic response. Note that the strain rate of an incompressible liquid for the radial ¯ow around a bubble is err ˆ

ovr ˆ or

2vr ; r

and its complex amplitude at bubble interface due to (22) is err rˆa ˆ

2ix

Da : a0

Thus, since S ˆ s…rr† tjrˆa ; s…rr†  g ˆ ; 2err rˆa i.e. such introduction of the complex viscosity is natural. For a polymeric liquid with single relaxation time from (24)±(27), we shall have gp ‡ gs ; g ˆ 1 ‡ its x gp gp ts x g0 ˆ ‡ gs ; g00 ˆ : 2 1 ‡ …ts x†2 1 ‡ …ts x†

…32† …33†

The expressions (33) agree with the expressions for components of the complex dynamic viscosity at periodic shear in Maxwell liquid with single relaxation time, except the additive to the real component corresponding to the solvent viscosity. The plot of the dynamic viscosity components g0 and g00 normalized to g0 against ts x at gs =g0 ˆ 0:001=1:48 is presented in Fig. 6 on a logarithmic scale. It is seen that the real component of viscosity g0 monotonously decreases from Newtonian viscosity of solution g0 to the viscosity of solvent gs . Curve g0 approaches to the horizontal asymptote g ˆ gs proportionally by 1=…ts x†2 while increasing ts x. The imaginary component of the complex dynamic viscosity g00 increases proportionally to ts x from zero to maximum g00 ˆ gp =2 at ts x ˆ 1, then decreases approaching zero as 1=…ts x†. It does not in¯uence the structure of the wave, since the corresponding dynamic elastic modulus G0 ˆ xg00 …G  G0 ‡ iG00 ˆ ixg † is signi®cantly lower than the reduced elastic modulus of gas in bubbles Gg (18), the dependence of which on ts x is shown in Fig. 6 via reduced imaginary viscosity of gas

648

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

Fig. 6. The components of complex dynamic viscosity at bubble oscillations in Maxwell liquid vs. frequency x multiplied by relaxation time ts .

g00g …ts x† ˆ

Gg 3 ˆ c2 pe ts =…ts x† 4 x

at various values ts ˆ 0:001, 0.1 and 1 ms. From considered ®gure and formula (32), one can see that frequency dependencies of dynamic modules (viscosities) on ts x at the radial oscillation of bubbles within the Maxwell model have the same form as at periodic shear. Thus, the data about the spectra of dynamic modules or viscosities obtained by routine methods at periodical shear may be used for approximate analysis and estimations of parameters of bubble oscillations in real polymeric solutions. 4.3. Discussion of simulation results The value of the viscosity realizing at bubble oscillations in polymeric solution in the shock wave shown in Fig. 2(a) is marked on the curve g0 by point D0 . This value is an order of two less than Newtonian viscosity of solution. It causes the appearance of oscillations in the shock wave traveling in the bubbly liquid with non-Newtonian carrier phase, unlike a Newtonian liquid with viscosity g0 (Fig. 2(b)). The plot of viscosity g0 allows to interpret qualitatively the results of calculations illustrated in Fig. 1. The points A0 ; B0 ; C 0 on curve g0 correspond to relaxation times and frequencies of bubble oscillations for the pro®les a; b; c (Fig. 1). It should be noted that for small relaxation times

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

649

ts x ! 0 (point A0 ), the viscosity realized is g0 and the waveform coincides with the waveform in the Newtonian liquid with viscosity g0 . For high relaxation times ts x ! 1 (point C 0 ), the viscosity of the non-Newtonian liquid is close to the viscosity of solvent gs and the waveform as in the Newtonian liquid with viscosity gs is realized. The point B0 corresponds to curve b (Fig. 1) and e€ective viscosity g0 lies between g0 and gs . Since this case is rheologically nonlinear, the pro®le of the wave in mixture with Newtonian carrier phase with the same viscosity g0 does not fully coincide with curve b. The simulations show that even if in some moment of evolution, the con®gurations of waves in non-Newtonian and Newtonian bubbly liquids with ``e€ective'' viscosity coincide, in the other moment of time they may be di€erent. The values of imaginary component of the viscosity corresponding to the wave pro®les a; b; c (Fig. 1) are marked on curve g00 (Fig. 6) by points A00 ; B00 ; C 00 . The points Ag ; Bg ; Cg show corresponding values of g00g and as one can see from Fig. 6 they lie signi®cantly higher than points A00 ; B00 ; C 00 . They have approximately equal ordinates, that is the consequence of weak dependence of frequency of bubble oscillations in a shock wave on time of stress relaxation. E€ective thermal viscosity. It is known (Gubaidullin et al., 1978) that for liquids of low viscosity with fairly large bubbles, the main dissipative mechanism at bubble oscillations, which determines the evolution of the shock wave, is heat transfer between gas and liquid. To compare the contribution of this mechanism with dissipation due to viscosity of liquid one can use the e€ective thermal viscosity (Nigmatulin, 1991). For the results in Figs. 1 and 2, the e€ective thermal viscosity estimated is approximately l…T† e  0:06 g0 . This was calculated with the formulas: v u …T† u …T† ae txe m2e …T† ; le ˆ l0 a0 x0 m…T† 20 s q 3 … c 1 † 1 3c2 p0 …T† 2 p q01 a0 m…T† l0 ˆ x0 ˆ ; 20 x0 ; a0 q01 4 2 …T†

where ae and m2e were used as adiabatic estimations of values behind the wave front ae ˆ a0 =…pe =p0 †1=…3c† ; …T†

…T†

m2e ˆ m20 =…pe =p0 †1=c ; and xe  3:2  104 1/s taken from Fig. 1. From Fig. 6, one can see that when the relaxation time is great, the evolution of wave is determined by thermal dissipation because g0 …C 0 †  le…T† . In the other cases, the role of viscosity is comparable to the e€ect of inter-phase heat transfer (point A0 and B0 ). In general, this relation varies. It depends on the time of stress relaxation and Newtonian viscosity of solution. This is con®rmed by computer simulations made for mixtures with a different kind of gas in bubbles (air, carbon dioxide, helium, etc.). 5. Conclusion The propagation of non-stationary shock waves in a bubbly liquid with a non-Newtonian carrier phase was investigated. A mathematical model describing the dynamic behavior of a

650

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

bubbly mixture with the non-Newtonian carrier phase was developed. The main conclusions are as follows: 1. The behavior of the shock wave in bubbly liquids with Newtonian and non-Newtonian carrier phase is fundamentally di€erent: the wave may have a monotonous structure in a Newtonian bubbly liquid but an oscillatory structure in a non-Newtonian one with the same Newtonian viscosity. It occurs because the reological losses may become much lower than the Newtonian viscosity due to the oscillations of bubbles. 2. Modeling of a polymeric bubbly liquid by a Newtonian one with ``e€ective'' viscosity so that the behavior of shock waves in these bubbly liquids is absolutely the same, is impossible in general case. But the polymeric bubbly solution can behave as a Newtonian bubbly liquid in two limiting cases: (a) If the time of stress relaxation is much less than the time scale of bubble oscillations, then the behavior is like that of the Newtonian liquid with the viscosity being equal to Newtonian viscosity of solution. (b) In the opposite limit, i.e. when the relaxation time is much greater than the time scale, the behavior is like that of the Newtonian liquid with the viscosity being equal to that of the solvent. This limit is purely realized only if the elastic response of polymer network is negligibly small as compared with that of the gas in bubbles; otherwise the polymer elasticity can change the behavior. 3. The attenuation of shock pulse in a non-Newtonian bubbly liquid is less than in a Newtonian one with viscosity being equal to Newtonian viscosity of solution. 4. The relation between thermal dissipation and dissipation due to viscosity of the carrier liquid is a variable value and it is de®ned by the time of stress relaxation and Newtonian viscosity of solution. 5. The frequency dependencies of dynamic modules (or complex viscosities) when radial oscillation of bubbles in the non-Newtonian liquid takes place, have the same form as in the case with periodical shear. Acknowledgements This work was supported by the Council on Program for State Support of Leading Scienti®c Schools (Grant No. 96-15-96001). Appendix A. Method of numerical modeling The problem of a one-dimensional non-stationary ¯ow of a single-velocity medium may be reasonably solved in Lagrangian coordinates …r; t† (Gubaidullin et al., 1978; Nigmatulin, 1991), where r is the distance from a particle to the origin at initial time t ˆ 0. The parameter values at t ˆ 0 are denoted by a subscript 0. In particular, q0 ˆ q0 …r† is the mixture density at t ˆ 0. The current position of the medium particle is characterized by the Eulerian coordinate x …r; t†, so that ox q0 ˆ ; or q

ox ˆ v: ot

…A:1†

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

651

As a result of the assumed simpli®cations, the system of Eqs. (1)±(3), (7)±(10) of motion for a monodisperse mixture of an incompressible liquid with gas bubbles in the absence of phase transitions in the framework of a single-velocity two-temperature approximation assumes the following form in Lagrangian variables: oq1 q1 q ov ‡ ˆ 0; ot q0 or wˆ a

oa ; ot

oq2 q2 q ov ‡ ˆ 0 …q ˆ q1 ‡ q2 †; ot q0 or

o 0 3
q a ˆ 0; ot 2

ow p2 ˆ ot

p1

a2 ˆ

2R=a ‡ S q01

S ˆ Sp ‡ Ss ;

Ss ˆ

ov 1 op ‡ ˆ 0; ot q0 or

4gs

w ; a

q2 ˆ1 q02

3 2 w; 2

oT2 3…c2 1† ˆ 3 ak2 Nu2 …T1 ot 2a0 cv2 q020

T2 †

q01 ˆ const:;

q1 ˆ a1 q01 ; …A:2†

oSp ‡ ot

p ˆ p1 a1 ‡ a2 …p2

q1 ; q01



 1 2w ‡ Sp ˆ ts a

4

gp w ; ts a

2R=a†; 3…c2

1†

w ; a

T1 ˆ T0 ˆ const:;

p2 ˆ q02 R2 T2 : The generalized Rayleigh±Lamb equation (the ®fth di€erential equation) may be transformed identically to a form in which, instead of average pressure p1 , the reduced pressure p appears a

ow p2 ˆ ot

p 2R=a q1

3 2 S w ‡ 0: 2 q1

…A:3†

With the prescribed physical properties of liquid and those of gas (adiabatic index c2 ,, gas constant R2 , coecient of heat conductivity k2 ), and also with the given parameter of interphase heat exchange Nu2 , and in the presence of initial and boundary conditions, the represented system of equations is closed. In connection with the condition of incompressibility, we shall perform some identical transformations. The ®rst two mass-conservation Eq. (A.2) of the system may be represented as oa1 a1 q ov ˆ 0; ‡ q0 or ot

oa2 a2 oq02 a2 q ov ˆ 0; ‡ 0 ‡ q0 or ot q2 ot

whence, upon summing these two equations, we obtain a2

oq02 q02 q ov ˆ 0: ‡ q0 or ot

Using the mass-conservation equation for the bubble (the fourth Eq. (A.2) in the form a3

oq02 ‡ 3a2 q02 w ˆ 0; ot

652

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

we obtain that the longitudinal deformation of mixture with an incompressible continuous phase takes place only due to the bubble radial deformation ov 3q0 a2 w ˆ : or q a

…A:4†

In this case, equations governing variation of both mixture density and bubble volumetric concentration may be written as q2 ov ˆ q0 or

oq ˆ ot

3q

a2 w ; a

oa2 a1 q ov 3a1 a2 w ˆ ˆ : ot q0 or a

The relation between the Lagrangian and Eulerian variables is given as 1 ov 1 ov d2 q02 oq02 : ˆ ˆ ; dt ot r q ox q0 or

…A:5†

…A:6†

Let us di€erentiate both the momentum equation with respect to r (recall that, q0 ˆ q0 …r††, and the continuity Eq. (A.4) with respect to t o2 v ˆ or@t

1 o2 p 1 dq0 op ‡ ; q0 or2 q20 dr or  o2 v w oa2 a2 ow a2 w oq ˆ 3q0 ‡ ot@r qa ot qa ot q2 a ot

 a2 w oa : qa2 ot

…A:7†

These two expressions in the domain of continuous motion may be set equal to each other. If, in this case, the time derivatives in the last equation are replaced by their values in accordance with both (A.5) and equations of radial motion, we obtain o2 p or2

K

op or

3q2 a2 M ˆ 02 qa

Lp ˆ M; 

p2

K…r† ˆ

1 dq0 ; q0 dr

Lˆ

 2R=a S w2 2 ‡ 0 ‡ 3 ‡ 2w : 2 a1 q01 q1

3q20 a2 ; q01 a1 a2 q01

…A:8†

Thus, the condition of incompressibility of continuous phase, as well as in hydrodynamics of a single-phase incompressible liquid, led to a second-order di€erential equation which contains derivatives with respect only to the spatial coordinate. This equation makes it possible to determine the pressure distribution at any point in time provided the remaining parameters (S; a2 ; a; w, and p2 ) are prescribed at the same point in time. It should be kept in mind that pressure p distribution and velocity v distribution at each point in time, including t ˆ 0, are not independent variables because of the incompressibility of the continuous liquid. The p distribution is determined from the above-indicated boundary-value problem in terms of distributions of a2 ; a; w; S and p2 and also in terms of boundary conditions on the ends r ˆ …0; L†, and the distribution of velocity …v† is found from Eq. (A.4), or from the momentum equation. The variation of a2 ; a; w; S and p2 in time is de®ned by di€erential equations (containing only time derivatives) following from Eqs. (A.2) and (A.5).

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

Let us introduce dimensionless variables denoted by bars pi a v w x  ˆ ; x ˆ ; Ti ˆ Ti =T0 ; pi ˆ ; a ˆ ; v ˆ ; w a C w La p0  0 r   C p C a0 1 0   ; Ca ˆ p ; La ˆ p ; C ˆ C10 ˆ ; a10 a20 a10 a20 Ca q01 w ˆ a0 =t ;

t ˆ a0 =C ;

k2 ˆ

k2 ; q020 c2 a0 C

g0 ˆ

653

q0i ˆ q0i =q0i0 ; …A:9†

g0 ; a0 q010 C

…A:10† R ; S ˆ p0 : a0 p0 To facilitate both the analysis and solution of the system of equations being discussed, we shall use dimensionless variables (A.9) together with parameters (A.10) denoted, as before, by a bar  q; q02 ; a); the latter are determined as a ratio of parameters corresponding to their characteristic (p; values (scales) designated by a subscript 0. Some ®xed values (P0 ; q01 ; q020 ; a0 ), being characteristics for distributions of p; q; ; q02 , and a at the initial time t ˆ 0, are chosen as scale factors. The latter, for convenience, are selected so that they satisfy the conditions of equilibrium R0 ˆ

p20 ˆ p10 ‡ 2R=a0 ;

q020 ˆ p20 =…R2 T0 †:

…A:11†

As before (see (A.9)), the equilibrium speed of sound Ce ; ˆ Ca is chosen as the inherent velocity of the longitudinal motion, and the speed C is chosen as the inherent velocity of radial motion near bubbles. The linear scale La is adopted the same as in (A.9). Accordingly, the time scale factor t0 is selected as   t r La a0 t ˆ ; r ˆ t0 ˆ ˆ : t0 La Ca C As a consequence, the system of equations may be reduced to the following form (let us drop the sign ``-'' above dimensionless variables) o2 p op Lp ˆ M; …A:12† K or2 or oUi ˆ Pi …i ˆ 1; 2; 3†; ot U1 ˆ w; U2 ˆ a; U3 ˆ T2 ;   1 p2 p 2R=a w 2 P1 ˆ 1:5w ‡ S ; S ˆ SS ‡ Sp ; SS ˆ 4gS ; a a1 a   gp w oSp 1 2w ; ‡ ‡ Sp ˆ 4 ts a ot ts a 3…c2 1† 3k2 P2 ˆ w; P3 ˆ aNu2 …1 T2 †; wT2 ‡ 0 a 2q20 cv2 1 dq0 3q2 a2 ; Lˆ 0 Kˆ ; q0 dr q a1 a2   3q20 a2 p2 2R=a w2 2 ‡ S ‡ 3 ‡ 2w : Mˆ a1 qa2 2

654

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

Note that in the case of one-dimensional ¯ow with plane waves, the velocity of the medium appears only in the second Eq. (A.4); therefore, to calculate it on each step of integration with respect to time is not necessary; and, in cases when the v distribution must be calculated, it is better to do it for each time interval of integration following the computation of p distribution based upon the momentum equation which in dimensionless variables is written in the form ov ˆ ot

a10 a20 op : q0 or

The obtained system consists of six di€erential equations, each of which contains only a single derivative with respect to one of the coordinates (r or t). The ®rst equation of the system is intended for determining the reduced pressure at an arbitrary time using the known ®elds of the rest of the parameters; the remaining equations describe the laws of variation of Lagrangian particles of the medium in time. For numerical integration of the obtained system of equations, we shall divide the selected volume of the medium by points r ˆ ri …i ˆ 1; 2; . . . ; n† into n material particles: the values of all sought functions we determine at points r ˆ ri …i ˆ 1; 2; . . . ; n†. Then, the last ®ve di€erential equations in partial time-derivatives of variables a; w; Sp ; T2 ; v are transformed into 5n ordinary di€erential equations with respect to time, for whose numerical integration, the modi®ed Euler± Cauchy method may be conveniently used. Corresponding di€erence equation can be written in the form  n ~ i ˆ …Ui †n ‡ …Pi †n Dt; U m m m h n i . n‡1 n ~ i ‡ …Pi †n Dt 2: …Ui †m ˆ …Ui †m ‡ P m m

Here n and m superscript and subscript correspond to the parameters in m ± sites of di€erence ~ j †; Dt is step of integration on time. ~ i ˆ Pi …U network at time moment n; P To ®nd the pressure p values at points r ˆ ri at each ®xed point in time, a linear (for p) boundary ± value problem must be solved for the ®rst di€erential equation (with respect to r) of the second order with boundary conditions o2 p or2

K

op or

Lp ˆ

M;

0 < r < L;

op …0; t† ˆ w0 …t†; or op d0 p…l; t† ‡ d1 …l; t† ˆ wl …t†: or l0 p…0; t† ‡ l1

To facilitate the solution of this problem the sweep method is advisable, in accordance to which the boundary-value problem is reduced to the solution of set of linear algebraic equations (the superscript n related to the network function will be dropped) by means of approximating the derivatives with respect to r by central di€erences

A.A. Gubaidullin et al. / International Journal of Multiphase Flow 27 (2001) 635±655

pm

1

Am pm ‡ Bm pm‡1 ˆ Fm ;

m ˆ 1; 2; . . . ; s

655

1;

p0 ˆ j1 p1 ‡ g1 ; ps ˆ j2 ps

1

‡ g2 :

The sweep method consists of two steps ± the direct and reverse ones. The direct step consists in calculating sweep coecients dm and bm : Bm b ; bm‡1 ˆ m Am dm Am d1 ˆ j1 ; b1 ˆ g1 : dm‡1 ˆ

Fm ; dm

m ˆ 1; 2; . . . ; s

1;

The reverse step consists in solving the problem j2 bs ‡ g2 ; pm ˆ dm‡1 pm‡1 ‡ bm‡1 ; 1 j2 ds m ˆ s 1; s 2; . . . ; 0:

ps ˆ

The advantage of the sweep method is the small number of algebraic operations and the low sensitivity to the error of calculation. It should be noted, for comparison, that the shorting method is inapplicable in this case because of its instability. References Astarita, G., Marrucci, G., 1974. Principles of Non-Newtonian Fluid Mechanics. McGraw-Hill, London. Ferry, J.D., 1980. Viscoelastic Properties of Polymers. Wiley, New York. Gubaidullin, A.A., Ivandaev, A.I., Nigmatulin, R.I., 1978. Nonsteady shock waves in gas-liquid mixtures of bubbly structure. J. Appl. Mech. Tech. Phys. 19, 204±210. Gubaidullin, A.A., 1991. The peculiarity of nonlinear waves evolution in bubbly liquids. In: Physical Acoustics. Fundamental and Applications, Plenum, New York, pp. 347±351. Ishiguro, S., Hartnett, J.P., 1992. Surface tension of aqueous polymer solutions. Int. Comm. Heat Mass Transfer 19, 285±295. Levitskiy, S.P., Shulman, Z.P., 1995. Bubbles in Polymeric Liquids: Dynamics, Heat and Mass Transfer. Technomic Publishing A G, Bassel, Switzerland. Nakoryakov, V.E., Pokusaev, B.G., Shreiber, I.R., 1990. Wave Propagation in a Gas±Liquid Media. CRC Press, Boca Raton, USA. Nigmatulin, R.I., 1991. Dynamics of Multiphase Media. Hemisphere, New York.