1 1 0 2 Combined effect of viscosity and vorticity on single mode n a J 0 2 Rayleigh -Taylor instability bubble growth ] h p - m Rahul Banerjee, Labakanta Mandal, S. Roy, M. Khan, M. R. Gupta ∗ † ‡ § ¶ s a l p Deptt. of Instrumentation Science & Centre for Plasma Studies . s c i s Jadavpur University, Kolkata-700032, India y h p [ 1 v Abstract 3 5 8 The combined effect of viscosity and vorticity on the growth rate of the bubble associated 3 . 1 0 with single mode Rayleigh -Taylor instability is investigated. It is shown that the effect of 1 1 viscosity on the motion of the lighter fluid associated with vorticity accumulated inside the : v i X bubble due to mass ablation may be such as to reduce the net viscous drag on the bubble r a exerted by the upper heavier fluid as the former rises through it. ∗e-mail: [email protected] †e-mail: [email protected] ‡e-mail: [email protected] §e-mail: mkhan [email protected] − ¶e-mail: mrgupta [email protected] − 1 I INTRODUCTION Recently obtained simulation results[1] regarding single mode Rayleigh Taylor instability(RTI) prob- lem indicates that the bubble developed at the two fluid interface is accelerated to a velocity quite above the classical model value[2] . This is ascribed to vorticity formation on the bubble-spike inter- face that is suggested[1] to diminish the friction drag due to viscosity. The theoretical calculations and the consequent numerical results regarding Rayleigh Taylor instability in a viscous fluid[3] also indicates growth of vorticity inside the bubble in the neighbourhood of its tip, where the curvature is maximum. Recent experimental and simulation results[4] show that the influence of viscosity on RTI may strongly suppress the growth rate. It is suggested that high pressure and strain rate condi- tions give rise to a phonon drag mechanism resulting in the lattice viscosity provided a solid state is maintained[4]. Moreover, it was also suggested by Betti and Sanz[5] that at the spike, mass ablation induces a transverse component of velocity and thus gives rise to vorticity generation. It has been shown by the same authors that the instability growth rate is augmented by vorticity accumulation inside the bubble resulting from mass ablation. As the bubble rises through the denser fluid it encounters viscous resistance. In this paper we have incorporated the joint effect of viscosity and vorticity. It is shown how the vorticity induced contribution to bubble velocity as derived by Betti and Sanz[5] is affected by the viscous drag exerted by the lighter fluid inside the bubble. It is found that the combined effect of vorticity and viscosity of the bubble fluid tends to diminish the friction drag exerted by the upper heavier fluid as the bubble penetrates through it. In this respect it is similar to the friction drag reduction effect due to vorticity 2 as suggested by the phenomenological model of Ramaprabhu et al.[1] . Section II contains the mathematical formulation as also the two fluid boundary conditions. The results and discussions are presented is section III. II MATHEMATICAL FORMULATION AND INTERFACIAL BOUNDARY CONDITION Weareconsidering 2Dsingle modeRayleigh-Taylor instability problem. They axis istaken vertically upwardoppositetothedirectionofgravityandxaxistakenalongthehorizontalplaneofunperturbed two fluid interface plane. The equation to the perturbed surface of bubble which is assumed to be parabolic is y η(x,t) = η (t)+η (t)x2; y > 0 (1) 0 2 ≡ with η (t) > 0 and η (t) < 0 (terms O(x2+i)(i 1) are neglected as we are interested only in the 0 2 ≥ motion very close to the tip of the bubble x 1)[6]. | | ≪ The kinematic boundary conditions satisfied on the perturbed interface given by Eq.(1) are ∂η ∂η +v = v (2) hx hy ∂t ∂x ∂η (v v ) = v v (3) hx lx hy ly ∂x − − where (v ) are the velocity components of the heavier (density ρ ) and lighter (density ρ ) fluids h,l x,y h l occupying the regions y > 0 and y < 0 respectively. 3 Assume single mode potential flow for the upper fluid governed by the potential φ (x,y,t) = a(t)cos(kx)e−k(y−η0(t)) (4) h with v = ∂φh and v = ∂φh hx − ∂x hy − ∂y Since both 2v = 0 and 2v = 0 for potential flow, the equation of motion of the upper hx hy ∇ ∇ incompressible fluid of constant density ρ and coefficient of kinematic viscosity ν (µ = ρ ν ) h h h h h ∂v~ ρ [ h +(v~ .~ )v~ ]+ ~(p +gρ y)+µ 2v~ =~0 (5) h h h h h h h ∂t ∇ ∇ ∇ leads to the following integral for the irrotational motion. ∂φ 1 ρ [ h + (~ φ )2 +gy]+p = f (t) (6) h h h h − ∂t 2 ∇ For the lighter density fluid the flow inside the bubble is assumed rotational[5] with vorticity ω~ = (∂vly ∂vlx)zˆ. The motion is described by the stream function Ψ(x,y,t); ∂x − ∂y ∂Ψ ∂Ψ v = and v = (7) lx ly −∂y ∂x so that 2Ψ = ω (8) ∇ − Let χ(x,y,t) be a function such that 2χ = ω (9) ∇ − 4 Hence (Ψ χ) is a harmonic function as 2(Ψ χ) = 0. Let Φ(x,y,t) be its conjugate function − ∇ − ∂Φ ∂Ψ ∂χ = ∂x ∂y − ∂y ∂Φ ∂Ψ ∂χ = + (10) ∂y −∂x ∂x i.e. that velocity components of the lighter fluid are ∂Ψ ∂Φ ∂χ v = = lx −∂y −∂x − ∂y ∂Ψ ∂Φ ∂χ v = = + (11) ly ∂x −∂y ∂x The following choice is made for the stream function according to [5] Ψ(x,y,t) = b (t)x+[b (t)ek(y−η0) +ω (t)/k2]sin(kx) (12) 0 1 0 with χ(x,y,t) = ω (t)sin(kx)/k2 (13) 0 Eq.(10) now gives Φ(x,y,t) = b (t)y +b (t)cos(kx)ek(y−η0) (14) 0 1 − Substituting for the velocity components as determined by the velocity potential φ and the stream h function Ψ in the kinematic boundary conditions Eq.(1) and Eq.(2) and equating coefficients of xi for i = 0 and 2 we obtain the following four equations for nondimensionalized bubble height ξ = kη , 1 0 nondimensionalized curvature ξ = η /k and b , b 2 2 0 1 dξ 1 = ξ (15) 3 dτ 5 dξ 1 2 = (6ξ +1)ξ (16) 2 3 dτ −2 kb 6ξ (2ξ ωˆ) 0 2 3 = − (17) √kg (6ξ 1) 2 − k2b (6ξ +1)ξ ωˆ 1 2 3 = − (18) √kg − (6ξ 1) 2 − Here k2a ξ = (19) 3 √kg is the nondimensionalized velocity of the tip of the bubble, τ = t kg (20) q is the nondimensionalized time and ω 0 ωˆ = (21) √kg is the nondimensionalized vorticity. We need one more equation, i,e., the one for ξ (τ) to complete the set of five necessary equations 3 for the five functions describing the time evolution of bubble. This is obtained from the dynamical boundary condition as done below. Starting from the equation of motion of the lighter (lower) fluid with coefficient of kinematic viscosity ν (µ = ρ ν ) l l l l ∂~v 1 ρ [ l + ~(~v 2)+ωzˆ ~v ]+ ~(p +ρ gy) µ 2~v = 0 (22) l l l l l l l ∂t 2∇ × ∇ − ∇ 6 and using Eqs.(7)-(11) we get ∂Φ 1 p ∂ω ∂χ˙ ∂ω ∂χ˙ ρ [ + (~v )2 ωΨ+ l +gy]+ ρ [(Ψ )dx+(Ψ + )dy] l l l − ∂t 2 − ρ ∂y − ∂y ∂x ∂x l Z ∂ω ∂ω + µ ( dx dy) = 0 (23) l ∂y − ∂x Z The last term on the L.H.S being the contribution from viscous drag. For rotational motion Bernoulli’s equation does not exist except for special cases[7]. However for velocity potential and stream function as given by Eqs.(12)-(14) the integrations in Eq.(23) can be expressed in closed form[5]. Taking the difference of Eq.(6) and Eq.(23) at the interface of the two fluids and applying the dy- namical boundary condition, i.e, the pressure difference at the interface is balanced by the difference of viscous stress and surface tension[8] we obtain the following: ∂φ 1 ∂Φ 1 ∂ω ∂χ˙ ∂ω ∂χ˙ ρ [ h + (~ φ )2 +gy] ρ [ + (~ Φ)2 ωΨ+gy] ρ [(Ψ )dx+(Ψ + )dy] h h l l − ∂t 2 ∇ − − ∂t 2 ∇ − − ∂y − ∂y ∂x ∂x Z ∂ω ∂ω ˜ µ ( dx dy) = [p p ]+f(t) l h l − ∂y − ∂x − − Z ∂v ∂v T hy ly ˜ = (2µ 2µ )+ +f(t) h l ∂y − ∂y R ∂2φ ∂2Ψ T h ˜ = 2µ +2µ + +f(t) (24) h l − ∂y2 ∂x∂y R where T is the surface tension and R is the radius of curvature at the interface. Substituting for φ ,Ψ,χ,Φ, using Layzer’s approximation method[9], expanding in power of x h to O(xi)(i 2) and equating coefficient of x2 we obtain after some laborious but straightforward ≤ 7 manipulation the time dependence equation for ξ (τ). The set of time evolution equations of the 3 bubble arethus given by [together withEq.(17) andEq.(18)] the following threedifferential equations for ξ ,ξ and ξ 1 2 3 dξ 1 = ξ (25) 3 dτ dξ 1 2 = ξ (6ξ +1) (26) 3 2 dτ −2 dξ ξ2 k2 1 3 = [ N(ξ ,r) 3 +2(r 1)(6ξ 1)ξ (1 12ξ2 )] dτ − 2 (6ξ 1) − 2− 2 − 2k2 D(ξ ,r) 2− c 2 ωˆ2 5(6ξ +1)ωˆξ 1 2 3 ˙ +[ − +ωˆ] (1 6ξ ) D(ξ ,r) 2 2 − 2ξ 2c rsωˆ(1+2ξ )(1 3ξ ) c r[(s+1)(1 12ξ2)+4ξ (s 1)] 3 + h 2 − 2 (27) − h − 2 2 − D(ξ ,r) D(ξ ,r) 2 2 where N(ξ ,r) = 36(1 r)ξ2 +12(4+r)ξ +(7 r) (28) 2 − 2 2 − D(ξ ,r) = 12(1 r)ξ2 +4(1 r)ξ +(r +1) (29) 2 − 2 − 2 (ρ ρ )g µ µ ρ ν k2 k2 = h − l ; s = l; ν = h; r = h; c = h (30) c T µ h ρ ρ h √kg h h l 8 III NUMERICAL RESULTS AND DISCUSSIONS The time evolution of bubble is to be worked out by numerical integration of Eqs.(25)-(27) by employing Runge-Kutta-Fehlberg scheme. To this end it is necessary to know the dependence of the vorticity ωˆ(τ) on τ. We choose ωˆ(τ) in the following form so as to have a time dependence having close resemblance with the simulation results[5]. ωˆ c ωˆ(τ) = [tanh(τ )(1+tanh(τ))+tanh(τ τ )] (31) 0 0 1+2tanh(τ ) − 0 Note that ωˆ(τ) increases from ωˆ(0) = 0 and gradually builds up to a saturation value ωˆ . The c constants τ and ωˆ are adjusted so that the time dependence has approximate qualitative agrement 0 c with simulation results given by Fig.4 of Ref.[5]. This may be seen from the graph of ωˆ(τ) plotted against τ in Fig.1. In Fig.2 are shown the nondimensionalized height (ξ ), curvature (ξ ) and velocity (ξ ) of the tip 1 2 3 of the bubble and plotted against time τ = √kg. (a). The first two terms [ N(ξ ,r) ξ32 +2(r 1)(6ξ 1)ξ (1 12ξ2k2)] 1 give the classical − 2 (6ξ2−1) − 2− 2 − 2kc2 D(ξ2,r) bubble velocity based on Goncharov model[2] modified by effect of surface tension. (b). The third term [ωˆ2−5(6ξ2+1)ωˆξ3 + ωˆ˙] 1 is the contribution from vorticity accumulation (1−6ξ2) D(ξ2,r) inside the bubble[5] causing enhancement of the bubble velocity. (c). The fourth term 2c r[(s+1)(1 12ξ2)+4ξ (s 1)] ξ3 accounts for the lowering of the − h − 2 2 − D(ξ2,r) bubble velocity due to viscous drag exerted by the denser fluid. (d). The last (fifth) term 2c rsωˆ(1+2ξ )(1 3ξ )/D(ξ ,r) gives the vorticity dependent contri- h 2 2 2 − bution that diminishes the friction drag; such a possibility was suggested by Ramaprabhu et al.[1] 9 . Eq.(26) and Eq.(27) show that as τ the bubble velocity reaches the following (nondimen- → ∞ sionalized) asymptotic steady value with ξ 1 (A= Atwood number), ωˆ(τ ) = ωˆ and 2 → −6 → ∞ c ˙ ωˆ(τ ) = 0, → ∞ 2 A ωˆ21 A 4 2 (v ) = + c − + c2 +sc ωˆ c (32) b asymp s31+A 4 1+A 9 h h c − 3 h It is be noted that this result agrees with the saturation value of Betti and Sanz[5] when viscosity is neglected(c = 0) and with that of Sung-Ik-Sohn[8] when vorticity is neglected. h In Fig.2 the physical process(a) is represented by plot I. It is the Goncharov (classical) model result 2 A (v )Gon = (33) b asymp s31+A In absence of viscosity the vorticity aided time development of the bubble velocity is given by plot II in Fig.2. This represents the combined effect (a)+(b). As ωˆ(τ) increases gradually from ωˆ(τ = 0) = 0, plot II indicates that so does the bubble velocity ξ (τ) from the initial value ξ (τ = 0) 3 3 remaining greater than the [ξ (τ)] but close to it as long as the vorticity does not increase 3 classical appreciably. This continues as τ increases till ωˆ(τ) builds up sufficiently as indicated in Fig.1 with a consequent rapid growth of ξ (τ) above the classical saturation value [ξ (τ)] as given in plot I 3 3 classical of Fig.2. As τ the vorticity aided bubble velocity approches the asymptotic value → ∞ 2 A ωˆ21 A (v )rot = + c − (34) b asymp s31+A 4 1+A 10