Solutions of the buoyancy-drag equation with a time-dependent acceleration
Serge E. Bouquet, Robert Conte, Vincent Kelsch, Fabien Louvet
aa r X i v : . [ n li n . S I] A ug Solutions of the buoyancy-drag equation with atime-dependent acceleration
Serge E. Bouquet , , Robert Conte , , , Vincent Kelsch and Fabien Louvet
1. CEA/DAM/DIF, Bruy`eres-le-Chˆatel, F–91297 Arpajon Cedex, France2. Laboratoire univers et th´eories (LUTH), Observatoire de Paris,Universit´e de recherche Paris sciences et lettres - PSL Research university,CNRS, Universit´e Paris-Diderot, Sorbonne Paris Cit´e,5, place Jules Janssen, F–92190 Meudon, FranceE-mail [email protected]. Centre de math´ematiques et de leurs applications,´Ecole normale sup´erieure de Cachan, CNRS,Universit´e Paris-Saclay,61, avenue du Pr´esident Wilson, F–94235 Cachan, France.4. Department of mathematics, The University of Hong Kong,Pokfulam Road, Hong Kong.E-mail [email protected] July 2017
Abstract
We perform the analytic study of the the buoyancy-drag equation with a time-dependentacceleration γ ( t ) by two methods. We first determine its equivalence class under the pointtransformations of Roger Liouville, and thus for some values of γ ( t ) define a time-dependentHamiltonian from which the buoyancy-drag equation can be derived. We then determinethe Lie point symmetries of the buoyancy-drag equation, which only exist for values of γ ( t )including the previous ones, plus additional classes of accelerations for which the equation isreducible to an Abel equation. This allows us to exhibit two r´egimes for the asymptotic (largetime t ) solution of the buoyancy-drag equation. It is shown that they describe a mixing zonedriven by the Rayleigh–Taylor instability and the Richtmyer–Meshkov instability, respectively. Contents
Interface instabilities such as Richtmyer–Meshkov instability, Kelvin–Helmholtz instability, orRayleigh–Taylor instability, are known to produce the mixing of two fluids living on each sideof a common interface [5, 23]. The Rayleigh–Taylor instability (RTI) plays, however, a special rolefor at least two reasons. First, from the physical viewpoint, RTI has been evidenced to occur invarious objects/processes with a spatial scale ranging from one millimeter for inertial confinementfusion laser targets to 10 meters for supernovae and typically 10 meters ( ≈
30 parsec) in su-pernova remnants. The second reason is related to the mathematical description. In the linearstage approximation, the formalism of the Richtmyer–Meshkov instability can be recovered fromthe RTI approach where the acceleration γ ( t ) ( t is time) experienced by the interface is restrictedto a kick (impulsive acceleration), and, moreover, the dispersion relation of the RTI can be eas-ily derived from the Kelvin-Helmholtz instability in a gravitational acceleration g provided g isformally replaced by γ ( t ) [5]. Because of the broad application of its formalism, we have decidedin this paper to study the mixing zone produced by the RTI. Considering a vertical downwardacceleration field, we take a two-fluid configuration where the high mass density material (heavyfluid), labeled by the subscript “2” lies above the low mass density medium (light fluid), labeledby the subscript “1” (see [4]). The fluids are assumed to extend vertically from −∞ to + ∞ with avelocity of the flow having only a vertical component depending on the position in the horizontalplane (two-dimensional configuration). Because of the RTI, spikes of “2” with length h ( t ) (thislength is measured from the initial position of the interface) drop into “1” while bubbles of “1”rise into “2” with elevation h ( t ) (measured also from the initial position of the interface betweenthe two fluids) in the nonlinear regime. The mixing zone corresponds to a dynamical region thatis bounded at its bottom by a surface intersecting the spike tips, and its upper border is given bythe surface joining the tops of the bubbles. Above (resp. below) the upper (resp. lower) frontier,the location of which depends upon time, the fluid “2” (resp. “1” ) is pure. In a zero-dimensionalmodeling, a spatial average is performed horizontally and the upper and lower boundaries reduceto two horizontal straight lines with respective position h ( t ) above and h ( t ) below the interfaceand the inner spatial structures of the flow cannot be described in between. Nevertheless, theheight h ( t ) of the mixing zone satisfies h ( t ) = h ( t ) + h ( t ) and a nonlinear ordinary differen-tial equation (ODE) governing the evolution of h ( t ) can be derived. This equation is called thebuoyancy-drag equation (BDE). Similarly to an equation of motion, the BDE is a second orderordinary differential equation that provides the length h ( t ) as a function of timed h d t = Bγ ( t ) − Ch (cid:18) d h d t (cid:19) (1)where B (buoyancy coefficient) and C (drag coefficient) are two positive constants (see [4]), andwhere γ ( t ) is the time-dependent acceleration experienced by the interface. More precisely, onehas 0 ≤ B ≤ ≤ C ≤ C and, more importantly, Bγ ( t ) allowingat least a partial integration of (1). The case of a time-independent (constant) acceleration hasbeen studied earlier in [8] and [4]. Although only the asymptotic (large value of time) form of thesolution was found by Cheng, Glimm and Sharp [8], the general solution of the BDE for a constantacceleration has been derived for the first time a few years later [4]. In the present work, our goalis not to obtain the general solution of (1) for any time-dependent acceleration but, instead, tofind classes of accelerations γ ( t ) for which the BDE can be fully or partially integrated.There exist three general methods to investigate a given second order ordinary differentialequation. The first one corresponds to the Painlev´e analysis [9]. It is based on those singularitiesof the general solution which depend on the initial conditions. The second approach applies tothe class of equations, to which the BDE belongs, investigated by Roger Liouville [20], and in thethird method the differential order of the BDE is reduced by using Lie symmetries.In our case, the first method cannot apply, since a necessary condition of its applicability is theexistence of a relative integer n such that C = − /n [9], but this condition never happens.However a useful property is the form invariance of equation (1) [24]. Indeed, the one-parameterpoint transformation ϕ a , T = att + a , h = (cid:18) tT (cid:19) / (1+ C ) H, γ = (cid:18) tT (cid:19) − (3+4 C ) / (1+ C ) Γ , a arbitrary , (2)only changes the acceleration, d H d T = B Γ( T ) − CH (cid:18) d H d T (cid:19) . (3)The inverse of ϕ a is ϕ − a and the composition of two such transformations does not generatea new transformation, since ϕ b ϕ a = ϕ ab/ ( a + b ) . In particular, a power-law accelation γ ∼ t n ismapped to another power-law γ ∼ t − α − n , with α = 3 + 4 C C = 4 −
11 + C , (4)the fixed point of such a map being γ ∼ t − α/ , which will indeed be encountered below.The paper is organized as follows. In section 2, we first recall the existing analytic solutionsof the BDE. In section 3, we build a Hamiltonian description of the dissipative equation (1) and,for specific t -dependences of the acceleration γ ( t ), find an invariant of the resulting Hamiltoniansystem. In section 4, for selected accelerations γ ( t ), we lower the differential order by one unit. Some solutions of Eq. (1) for the RTI are already known in four cases, the first three have beenderived using Lie symmetries [4, 24, 18]. Throughout the paper, t denotes some time unit.1. For a constant acceleration γ = γ , the BDE is autonomous and admits the first integral γ = γ : K = h C "(cid:18) d h d t (cid:19) − Bγ C h . (5)The BDE is then invariant under the two symmetries ∂ t and t∂ t + 2 h∂ h (the subscriptsrepresent the partial derivatives), which allows us to obtain its general solution in the implicitform [4] [14, formula 3.194] γ = γ , t = t + h C (1 + C ) √ εK F (cid:18) , C C , C C , − Bγ h C (1 + 2 C ) εK (cid:19) , ε = sign( K ) , (6)in which F is the hypergeometric function of Gauss, and the two arbitrary constants are t and K . 3his expression yields the asymptotic behavior of the height as t → + ∞ h ∼ Bγ C ) t . (7)Such a time-dependence has already been found earlier by Neuvazhaev [26].This quadratic law is consistent with the fact that in a constant acceleration field the freefall distance scales as t . This result is not really surprising because the friction term − ( C/h )(d h/ d t ) scales as t and the BDE becomes a free fall equation.In a constant acceleration field γ without friction ( C = 0) , the general solution of Eq. (1)behaves like h ∼ ( Bγ / t , and the comparison with (7) shows, however, that the actualsolution accounts for friction effects. We could have expected a balance between the ac-celeration term Bγ and the friction term − ( C/h )(d h/ d t ) , i.e. h ∼ [ Bγ / (4 C )] t , in theasymptotic evolution. Nevertheless, Eq. (7) proves that this condition is not fulfilled and,instead, the balance between the three terms (acceleration, friction and inertia d h/ d t ) isachieved in the BDE.2. For the power law acceleration γ ∼ t − α , the general solution is simply obtained from (6) byaction of the transformation (2) to generate the case Γ = γ in (3).As expected, for large time t , the behavior h ∼ t is no longer valid. The actual asymptoticsolution is obtained by noting that lim t → + ∞ T = a . As a consequence, T and H are boundedwith H → H ( a ). Using (2), we get h ∼ ( t/a ) / (1+ C ) H ( a ), i.e. h ∼ t / (1+ C ) , t → + ∞ . (8)This behavior is actually very surprising for the Rayleigh–Taylor mixing since such a powerlaw is relevant to the Richtmyer–Meshkov instability (no acceleration or impulsive acceler-ation when a shock wave front reaches an interface). Indeed, for a zero acceleration in (1),the BDE reduces to d h d t = − Ch (cid:18) d h d t (cid:19) , (9)and the general solution is γ = 0 : h RM = K ( t − t ) / (1+ C ) (10)where the subscript RM stands for Richtmyer–Meshkov and where K and t are the twoconstants of integration. Asymptotically (large time t ), we get h RM ∼ t / (1+ C ) , t → + ∞ . (11)As a consequence, the physical interpretation of (8) is as follows: for n = − α , the accelera-tion γ ( t ) decreases very quickly with time and the mixing regime is not driven by the RTIat large time. Instead, the system experiences the Richtmyer–Meshkov instability (RMI)asymptotically because γ → t → + ∞ .Finally, as in the linear stage of the RMI we have a ballistic growth of the mixing zone, weget h ∼ t . However, (8) shows that in the nonlinear regime, h increases more slowly thanlinearly in t (the exponent 1 / (1 + C ) is always smaller than 1) and the observed decelerationis produced by the friction term − ( C/h )(d h/ d t ) .3. For a general power law acceleration γ ∼ t n with an arbitrary real exponent n , the BDE (1)only admits the symmetry t∂ t + (2 + n ) h∂ h , which allows its reduction to an Abel equation[4], γ = γ (cid:18) tt (cid:19) n , I = (cid:18) tt (cid:19) − n − h, J = (1 + C ) (cid:18) t d hh d t − n − (cid:19) ,IJ d J + (1 + C ) "(cid:18) J + (1 + C ) n + 3 + 4 C (cid:19) − − (1 + C ) Bγ t I d I = 0 . (12)4ince this Abel ODE is in general not integrable [7, 25, 30, 32], this only defines thezero-parameter scaling solution, γ = γ (cid:18) tt (cid:19) n ,h = h (cid:18) tt (cid:19) n +2 , h = Bγ t ( n + 2)[ n (1 + C ) + 1 + 2 C ] , ( n + 2) (cid:18) n + 1 + 2 C C (cid:19) = 0 . (13)In these relations, the exponent n must satisfy two major constraints. First, as the mixingzone thickness increases with time, the exponent n should obey n + 2 >
0. Second, thecoefficient h has to be positive and, since n > −
2, we get n > − C C = − C . (14)As a consequence, one expects the solution h ∼ t n +2 to be valid for n either positive ornegative, provided its value is larger than the lower bound given by Eq. (14). For n >
0, theacceleration grows with time and the thickness of the mixing zone increases. For the case n <
0, although the acceleration vanishes for large time, the thickness of the mixing zoneincreases too. In both cases, the mixing is driven by the RTI.Moreover, for the lower bound (14), we have n +2 = 1 / (1+ C ) and the time dependence in (8)and (13) coincide. According to (11), this specific time variation of the acceleration producesthe same asymptotic mixing rate for the RTI than for the RMI. Finally, for n = 0 we have h = Bγ t / [2(1+2 C )] and therefore h = Bγ t / [2(1+2 C )]. This solution corresponds to (7)and we conclude that the particular solution (13) provides the correct asymptotic behaviorof the RTI for a constant acceleration.If the acceleration decreases with t faster than prescribed by the lower bound value (14),the solution h ∼ t n +2 does not happen anymore. According to our explanations, we mightexpect that the asymptotic solution of the BDE be given by (11). This prediction actuallyarises for γ ∼ t − α [see item (2) above]: the value n = − α = − − [(1 + 2 C ) / (1 + C )] isalways smaller than − h ( t ) behaves like t / (1+ C ) instead of t n +2 = t − (1+2 C ) / (1+ C ) that decreases with time.In the next item, this prediction is also proven to hold for another negative value of n .4. For the value of n which leaves the behaviour γ ∼ t n invariant under the map ϕ a , i.e. n = − α/
2, the Abel equation (12) becomes linear in J and there exists a first integral n = − α , γ = γ (cid:18) tt (cid:19) n , K = (cid:20) J − − C ) C Bγ t I (cid:21) I C ) . (15)Moreover, when K vanishes, the first order ODE for h ( t ) can be integrated. Indeed, thechange of function ( h, t ) → ( H, T ) : h = T H, T = (cid:18) tt (cid:19) n +2 , (16)maps the first integral to T (cid:18) d H d T (cid:19) − H ( H − H ) = 0 , H = − Bγ t (1 + C ) C , (17)an equation linearizable by derivation. This ODE has two kinds of solutions [6]: the so-calledsingular solution H = H , h = H ( t/t ) − α/ , which must be rejected because it is neversolution of the BDE, and the general solution [24], H = ( c T − / + c T / ) , c c = H c c arbitrary , (18)5 .e. n = − α , γ = γ (cid:18) tt (cid:19) n , h = " c + c (cid:18) tt (cid:19) n +2 , c c = H · (19)The latter solution depends on the arbitrary parameter c /c but it never reduces to theprevious solution (13).The exponent n + 2 in (19) is always positive and hence, for large time t , the constant c canbe neglected in (19) showing that h grows like t n +2) , which is, h ∼ t / (1+ C ) . (20)Interestingly, the behavior (8) is recovered.This result is not surprising because the exponent n = − α/
2, which can be written as thesum n = − / (1 + C ) − / [2(1 + C )] is therefore smaller than the lower bound (14).We conclude that, although the acceleration decreases with t more slowly for n = − α/ n = − α , the function γ ( t ) vanishes rapidly enough for the system to exhibit a mixingzone driven by the RMI instead of the RTI. The approach developed by R. Liouville [20] applies to the class of equationsd h d t + a ( h, t ) (cid:18) d h d t (cid:19) + 3 a ( h, t ) (cid:18) d h d t (cid:19) + 3 a ( h, t ) d h d t + a ( h, t ) = 0 , (21)whose property is to be form invariant under the point transformation( t, h ) → ( T, H ) : t = F ( T, H ) , h = G ( T, H ) . (22)By determining the invariants of (21) under the transformation (22), one may be able to integrate.A nice account of this method can be found in [3], to which we refer for the notation. In the caseof (1), one obtains L = CBγ ( t ) h − , L = 0 , ν = 0 , w = 0 , i = C (2 + C ) Bγ ( t ) h − , (23)therefore (see e.g. [3, Lemma 1 page 458]) the ordinary differential equation (1) can be mapped toan ODE (21) in which a = a = a = 0, t = t, h = H / (1+ C ) , d H d t − (1 + C ) Bγ ( t ) H C/ (1+ C ) = 0 . (24)Since this equation is independent of d H/ d t , it can be interpreted as the Hamilton equation ofa classical time-dependent Hamiltonian H defined by H ( q, p, t ) = p V ( q, t ) , q = H, p = d q d t , V = − (1 + C ) C Bγ ( t ) q C C (25)where q , p and V are respectively the position, linear momentum and potential. Under thistransformation, equation (24) becomes a standard equation of motiond q d t + ∂∂q V ( q, t ) = 0 , (26)for which we are going to look for an invariant of motion I ( q, p, t ).The invariant obeys the equation ∂ I /∂t + p ( ∂ I /∂q ) − ( ∂V /∂q )( ∂ I /∂p ) = 0 and since H isquadratic in p , it is natural to look for I ( q, p, t ) also quadratic in p . Such an invariant does exist, I ( q, p, t ) = ( d t + 2 d t + d ) H ( q, p, t ) − ( d t + d ) qp + d q d , d , d ) arbitrary , (27)6owever only for specific accelerations γ ( t ) defined by, (cid:0) d t + 2 d t + d (cid:1) d γ d t + ( d t + d ) αγ = 0 , α = 3 + 4 C C · (28)According to the values of the d i ’s, three types of t -dependence for γ ( t ) come out, γ = γ (cid:18) t − t t (cid:19) − α/ ( d = 0) ,γ (cid:18) tt (cid:19) − α/ ( d = 0 , d = 0) ,γ ( d = d = 0 , d = 0) , (29)where t is an arbitrary real constant of any sign.The first expression in (29) has been obtained by performing a time translation t → t − d / (2 d )(which leaves the BDE invariant), and t = ( d ) / (4 d ) − d /d .This expression is an extension of the second case listed in section 2 and it is recovered for t = 0.The last two values correspond to the first and fourth cases listed in section 2 (the second valuerequires the shift t → t − d /d ).Le us now show that one can recover these three values by another approach, and even findmore general ones. Given any partial differential equation E ( x, t, u ( x, t ) , u x , u t , . . . ) = 0, a Lie point symmetry is atransformation( x, t, u ) → ( X, T, U ) : X = F ( x, t, u ) , T = G ( x, t, u ) , U = H ( x, t, u ) , mapping a solution u ( x, t ) to another solution U ( X, T ).Practically, instead of this finite transformation, one computes the infinitesimal transformation X = x + εζ ( x, t, u ) , T = t + εξ ( x, t, u ) , U = u + εη ( x, t, u ) , (30)associated to the infinitesimal point symmetry S = ζ∂ x + ξ∂ t + η∂ u where the subscripts stand forthe partial derivatives [28, 29].In the case of (1), the assumption T = t + εξ ( t, h ) , H = h + εη ( t, h ) (31)yields the set of determining equations for ξ ( t, h ) , η ( t, h ) [18] ξ hh − ( C/h ) ξ h = 0 ,η hh + ( C/h ) η h − ( C/h ) η − ξ th = 0 ,ξ tt + 3 Bγ ( t ) ξ h − C/h ) η t − η th = 0 ,η tt − Bγ ( t ) ξ t + Bγ ( t ) η h − Bγ ′ ( t ) ξ = 0 . (32)where γ ′ ( t ) is the time derivative of γ ( t ) . Remark . The above assumption (31) for Lie point symmetries does not act on γ ( t ), thereforeit cannot detect the form invariance (2).The system (32) is a linear overdetermined system for ξ ( t, h ) , η ( t, h ), solved as follows.The first equation is an ODE for ξ ( h ) (with t a parameter) having the type of Fuchs, whosegeneral solution depends on two arbitrary functions of t , ξ = C ( t ) h C C + C ( t ) . (33)7he second equation is then an ODE for η ( h ) (with t a parameter) of the type of Fuchs, whosegeneral solution introduces two more arbitrary functions of t , η = C ( t ) h + C ( t ) h − C + C ′ ( t ) h C (1 + C ) (34)where C ′ ( t ) stands for the time derivative of C ( t ).Inserting these two expressions in the last two equations (32) puts constraints on C j ( t ) and γ ( t ).Finally, ξ and η depend on four arbitrary constants d i [18] ξ = d t + 2 d t + d , η = d t + d + d C h, and γ ( t ) must obey a first order linear ODE, (cid:0) d t + 2 d t + d (cid:1) d γ d t + [( d t + d + d ) α − d ] γ = 0 (35)whose only fixed parameter is the positive constant α . This ODE is an extension of (28).This four-dimensional algebra, generated by X = t ∂ t + th∂ h , X = 2 t∂ t + h∂ h , X = ∂ t , X = h∂ h , (36)is decomposable into { X , X , X } ⊕ X and its nonzero commutators are,[ X , X ] = − X , [ X , X ] = − X , [ X , X ] = − X . (37)Since the ODE (35) contains one more parameter than the similar ODE (28) resulting fromthe Hamiltonian structure, the set of values of γ ( t ) is now larger, γ = γ (cid:18) t − t t (cid:19) n (cid:18) t + t t (cid:19) n , n + n = − α, n − n = α t t ( d = 0) ,γ (cid:18) tt (cid:19) n ( d = 0 , d = 0) ,γ e αt/t ( d = d = 0 , d = 0 , d = 0) ,γ ( d = d = d = 0 , d = 0) , (38)in which t , t , t and n are arbitrary real constants, and time has been shifted like in (29) in orderto derive the first three expressions.The second and fourth values of γ ( t ) yield the solutions already mentioned in section 2, butthe first and the third ones are new.The first value defines a new case of reduction to an Abel equation in which the invariants I and J are given by d γγ d t = − α t − t t − t , t and t = arbitrary real constants ,γ = γ (cid:18) t − t t (cid:19) n (cid:18) t + t t (cid:19) n , n + n = − α, n − n = αt /t ,I = (cid:18) t − t t (cid:19) − n − (cid:18) t + t t (cid:19) − n − h,J = 1 t (cid:20) (1 + C )( t − t ) d hh d t − t − (3 + 4 C ) t (cid:21) ,IJ d J d I + (1 + C ) "(cid:18) J + (3 + 4 C ) t t (cid:19) − (cid:18) t t (cid:19) − (1 + C ) Bγ t I = 0 . (39)This solution contains the first expression of the list (29) as a special case for n = n = − α/ h = (1 + C ) Bγ t (3 + 4 C ) t − t (cid:18) t − t t (cid:19) n +2 (cid:18) t + t t (cid:19) n +2 , (40)8nd it depends on the two arbitrary constants t , t .Let us examine the asymptotic form of this solution. For large time ( t ≫ t ), the accelerationand the height are γ ∼ t n + n , h ∼ t n + n +4 , t → + ∞ . (41)Obviously, this behavior does not belong to the class (13) with n = n + n . However, with (4)and the constraint n + n = − α , we obtain h ∼ t / (1+ C ) which not surprisingly is identical to (8).The RMI growth is again recovered.Similarly, the third value in the list (38) also defines a reduction to an Abel equation γ = γ e αt/t , t = 0 I = e − αt/t h,J = (1 + C ) (cid:18) t d hh d t − α (cid:19) ,IJ d J d I + (1 + C ) (cid:20) ( J + 3 + 4 C ) − (1 + C ) Bγ t I (cid:21) = 0 (42)and the particular solution h = Bγ t α (1 + C ) e αt/t . (43)It is clear that this solution arises only for t >
0: if t <
0, the acceleration decreases exponentiallywith time and after a transient phase, h will follow the RMI law (8) instead of (43).The three ODEs for J ( I ) defined in (12), (39) and (42) are Abel equations of the second kind, IJ d J d I + a (cid:20) ( J − b ) − c − dI (cid:21) = 0 , (44)for the respective values a = 1 + C, d = (1 + C ) Bγ t , c = 12 , b = − (1 + C ) n − C , (45) a = 1 + C, d = (1 + C ) Bγ t , c = t t , b = − (3 + 4 C ) t t , (46) a = 1 + C, d = (1 + C ) Bγ t , c = 0 , b = − (3 + 4 C ) · (47)The classical method to investigate the integrability of an Abel ODE is recalled in the A. Inour case (44), the first two relative invariants of Roger Liouville evaluate to s = 2 a b I [ a ( b − c ) I − d ( a − ,s = 2 a b I (cid:2) a ( b − c )( b + 3 c ) I − d ( a − b a + 18 c a − b + 3 c ) I + 3 d (3 a − (cid:3) . (48)Since the absolute invariant s /s is generically nonconstant, there exists no change of variablesmaking Eq. (44) separable.The conditions s = s = 0 define two nongeneric cases, b = 0 and ( a = 1 , b = 9 c ), the secondone being nonphysical. Moreover, the Maple package developed in [7], which knows all the AbelODEs which have been integrated before the year 2000, fails to map (44) to one of the integratedequivalence classes, except in the unphysical case d = 0.9o summarize, the only cases when a first integral K is known to Eq. (44) are, b = 0 , a = 1 / K = log I + log[(2 a − J − c ) − ad/I ]2 a ,b = 0 , a = 1 / K = log I + ( J − c ) Id ,a = 1 , c = b : K = ( c + b ) log[ d + ( b − c ) I ( J − b − c )] + ( c − b ) log[ d + ( b + c ) I ( J − b + c )] ,a = 1 , c = b : K = log[2 bIJ + d ] − bd I ( J − b ) ,d = 0 , c = 0 : K = log I + ( c + b ) log[ J − b − c ] + ( c − b ) log[ J − b + c ]2 ac ,d = 0 , c = 0 : K = log I + 1 a log( J − b ) − ba ( J − b ) . (49)The first case mentioned in the list (49) corresponds to the solution (15).Various assumptions extrapolating the three cases b ( a − d = 0, such as K = log f ( I ) + k + log[ J + f + ( I )] + k − log[ J + f − ( I )] , (50)have failed to provide any new integrable case. Six closed form solutions ( γ ( t ) , h ( t )) of the BDE have been presented in this article: (6) and itssister solution for γ = γ ( t/t ) − α , (13), (19), (40), (43).If one excludes the unit of time t and the unit of acceleration γ , the number of their arbitraryparameters is respectively two ( t and K ), one ( n ), one ( c /c ), two ( t , t ) and one ( t ).As far as we know, sister solution of (6) and solutions (19), (40) and (43) are new.Moreover, there exist two values of γ ( t ) allowing the BDE to reduce to a first order ODE of thetype of Abel. The only hope to integrate this Abel equation (44) is to guess an integrating factor,necessarily outside the classes (57) already examined by Abel, Liouville and Appell [2].In spite of its simplicity, solution (13) has been found to play a key role. It has helped to exhibittwo families of solutions for the spatial extension h ( t ) of the mixing zone when the accelerationbehaves asymptotically ( t → + ∞ ) like γ ( t ) ∼ t n where n is an arbitrary positive or negativeexponent. For n >
0, the growth of the acceleration is monotonic with time and h increaseslike h ( t ) ∼ t n +2 for t → + ∞ . For n <
0, the acceleration decreases asymptotically with time,and for large negative values of n , equation (9) shows that the asymptotic solution does notdepend anymore on n but is given by h RM [see Eq. (11)]. In this paper, we claim that h ( t ) ∼ t n +2 describes an acceleration-driven mixing, i.e. a mixing of Rayleigh–Taylor type while h ( t ) ∼ t / (1+ C ) corresponds to an acceleration decreasing too fast with time in order to be able to drive theevolution of the mixing zone and the solution of the BDE corresponds to a Richtmyer–Meskhovmixing. The BDE exhibits therefore two leading behaviors and the way each of them arises isexplained as follows: first, we notice that for n ∈ [ − , − (1 + 2 C ) / (1 + C )], solution (13) is notphysically valid since h < n ∈ ] − , − (1 + 2 C ) / (1 + C )[ and h diverges for the twovalues n = − n = − (1 + 2 C ) / (1 + C ). Second, we show that for n = − α (this value isalways below the range [ − , − (1 + 2 C ) / (1 + C )]) and n = − α/ − , − (1 + 2 C ) / (1 + C )]), the asymptotic solution is h ( t ) ∼ t / (1+ C ) .As a consequence, we conclude that the threshold value for n is n th = − (1 + 2 C ) / (1 + C ) (lowerbound value (14)) : for n < n th , the acceleration decreases quickly with time and the systemexperiences the RMI with an asymptotic solution given by h ( t ) ∼ t / (1+ C ) ; for n > n th , the systemis driven by the RTI with the asymptotic solution h ( t ) ∼ t n +2 .After submission of the present manuscript, we noticed that these two regimes had been alreadyevidenced by Pandian, Swisher and Abarzhi under the name “acceleration-driven mixing” and“dissipation-driven mixing” for n > n th and n < n th , respectively [31].We notice that the RTI mixing occurs therefore also for n th < n < γ ( t ) decreaseswith time. Finally, for n = n th , the two regimes collapse in a single one since t n +2 = t / (1+ C ) : thegrowth rates of the RTI and the RMI are the same. According to us, this claim is an importantissue that would be worth being investigated through numerical simulations.10he solutions exhibited in this work provide an extension of the self-similar variable accelerationRayleigh–Taylor (SSVART) flows studied by Llor [22, 23] where the author uses an acceleration γ ( t ) ∼ t n for t > n > −
2, and γ ( t ) ∼ ( − t ) n for t < n < − n > − n = − n , the direct application of (13) leads to h ∼ t n +2 ∼ t which of courseis physically wrong although the BDE is satisfied. Actually, this special value of n satisfies n < n th and the correct asymptotic behavior is h ∼ t / (1+ C ) .The case n = 0 describes the standard RTI ( i.e. constant acceleration) and, as expected, (7) isrecovered from the particular solution (13). For a constant acceleration, h ( t ) is usually written as h ( t ) ∼ αγ t [13] and the comparison with (7) leads to the analytical value α = ( B/ / (1 + 2 C ).For B = 1, the case C = 4 leads to the smallest value 1 / ≈ .
055 for α in agreement withexperiments whereas the value derived from numerical simulations is about twice smaller [13].Finally, let us examine the case n = −
1. This value is above n th and one expects the asymptoticsolution h ∼ h ( t/t ) where the constant h is given by (13). Indeed, d h/ d t = 0 for this solutionbut the balance between Bγ ( t ) and − ( C/h )( dh/dt ) is satisfied in the BDE. Acknowledgments
RC was partially supported by the Laboratoire de recherche conventionn´e LRC-M´eso.
A Abel equation
We recall in this Appendix how to integrate an Abel equation.The most general Abel equation (which can always be assumed to be of the first kind), − d u d x + a ( x ) u + a ( x ) u + a ( x ) u + a ( x ) = 0 , (51)is form-invariant under the mapping( u, x ) → ( U, X ) : x = F ( X ) , u = P ( X ) U ( X ) + Q ( X ) , F ′ P = 0 . (52)The classical method to decide whether it is integrable or not has been introduced by RogerLiouville [21]. This first requires to compute its invariants under this mapping, which are rationalfunctions of the a j ’s and their derivatives. The relative invariants of (51) have weight 2 m + 1, with m a positive integer, s = a a + 13 (cid:18) a − a a a + a a ′ − a a ′ (cid:19) , (53) s m +1 = a s ′ m − − (2 m − s m − (cid:18) a ′ + a a − a (cid:19) , (54)and the absolute invariants I n are I = s s , I = s s s , I = s s , . . . (55)If I is constant, all other absolute invariants are also constant and there exists a mapping (52)making the transformed equation separable.If I is not constant, integrating is equivalent to finding an integrating factor µ ( u, x ) to (51).According to Abel, it is then more convenient to first put the equation under the form u d u d x + p ( x ) + q ′ ( x ) u = 0 , (56)(whose advantage is to minimize the global degree in u and u ′ ), then to make various assumptionsfor µ , such as those considered by Abel,log µ = ( b ( x ) u + b ( x )) − , µ = ( b ( x ) u + b ( x )) n , µ = ( u + b ( x )) a ( u + b ( x )) b , log µ = b ( x ) u + b ( x ) u + b ( x ) u + b ( x ) , µ = [7, p. 215 Eq. (48)]) . (57)11ut the primary use of the invariants is to decide whether the equation belongs to the sameequivalence class than one of the Abel ODEs previously integrated. An outstanding presentationof the current achievements can be found in [7], where the authors reviewed all the cases listedin the books [17], [25] [32], added a few ones and ordered them in different classes of equivalence modulo (52).Let us mention a different approach [30] to try to integrate Abel ODEs. The authors split thesingle ODE into a system made of one Riccati ODE and another condition. If the Riccati ODEcan be integrated, then this may lead to the integration of the Abel ODE. However, this methodhas not yet produced new integrated Abel equations. References [1] U. Alon, J. Hecht, D. Hofer and D. Shvarts, Power laws and similarity of Rayleigh-Taylor andRichtmyer-Meshkov mixing fronts at all density ratios, Physical review letters (4) (1995)534–537.[2] P. Appell, Sur les invariants de quelques ´equations diff´erentielles, Journal de math´ematiquespures et appliqu´ees (1889) 360-424.[3] M.V. Babich and L.A. Bordag, Projective differential geometrical structure of the Painlev´eequations, J. differential equations (1999) 452–485.[4] S. Bouquet, P. Gandeboeuf and P. Pailhori`es, Analytic study of the buoyancy-drag equation,Math. meth. appl. sci. (2007) 2027–2035.[5] S. Chandrasekhar, Hydrodynamic and hydromagnetic stability (Oxford university press, Ox-ford, 1961).[6] J. Chazy, Sur les ´equations diff´erentielles du troisi`eme ordre et d’ordre sup´erieur dontl’int´egrale g´en´erale a ses points critiques fixes, Acta Math. (1911) 317–385.[7] E.S. Cheb-Terrab and A.D. Roche, Abel ODEs: equivalence and integrable classes, Computerphysics communications (2000) 204–231.[8] B. Cheng, J. Glimm and D.H. Sharp, Dynamical evolution of Rayleigh-Taylor and Richtmyer-Meshkov mixing fronts, Physical review E (3) (2002) 036312 (7 pages).[9] R. Conte and M. Musette, The Painlev´e handbook (Springer, Berlin, 2008). (Regular andchaotic dynamics, Moscow, 2011).[10] R.M. Davies and F.R.S. Sir Geoffroy Taylor, The mechanics of large bubbles rising throughextended liquids and through liquids in tubes, Proceedings of the royal society of London,Ser. A (1950) 375–390[11] G. Dimonte, Spanwise homogeneous buoyancy-drag model for Rayleigh-Taylor mixing andexperimental evaluation, Physics of plasmas (6) (2000) 2255–2269.[12] G. Dimonte and M. Schneider, Density ratio dependence of Rayleigh–Taylor mixing for sus-tained and impulsive acceleration histories, Physics of fluids (2) (2000) 304–321.[13] G. Dimonte, D.L. Youngs, A. Dimits, S. Weber, M. Marinak, S. Wunsch, C. Garasi, A. Robin-son, M.J. Andrews, P. Ramaprabhu, A.C. Calder, B. Fryxell, J. Biello, L. Dursi, P. MacNeice,K. Olson, P. Ricker, R. Rosner, F. Timmes, H. Tufo, Y.-N. Young and M. Zingale, A compara-tive study of the turbulent Rayleigh–Taylor instability using high-resolution three-dimensionalnumerical simulations: the Alpha-group collaboration, Physics of fluids (5) (2004) 1668–1693.[14] I.S. Gradshteyn and I.M. Ryzhik, Tables of integrals, series, and products (Academic press,New York, 1980).[15] B.J. Gr´ea, The rapid acceleration model and the growth rate of a turbulent mixing zoneinduced by Rayleigh-Taylor instability, Physics of fluids (2013) 015118, 1–20.1216] J. Hecht, U. Alon and D. Shvarts, Potential flow models of Rayleigh–Taylor and Richtmyer–Meshkov bubble fronts, Physics of fluids (12) (1994) 4019–4030.[17] E. Kamke, Differentialgleichungen: L¨osungsmethoden und L¨osungen , Vol. 1, 243 pages; vol. 2,668 pages. Akademische Verlagsgesellschaft, Geest & Portig k.-G., Leipzig 1947. Chelsea, NewYork, 1948.[18] V. Kelsch, S. Bouquet et R. Conte, Calcul formel des sym´etries de Lie et application `al’´equation de pouss´ee-traˆın´ee, Rapport interne CEA-DIF (to appear).[19] D. Layzer, On the instability of superposed fluids in a gravitational field, The astrophysicaljournal (1) (1955) 1–12.[20] R. Liouville, Sur les invariants de certaines ´equations diff´erentielles et sur leurs applications,Journal de l’´Ecole polytechnique (1889) 7–76.[21] R. Liouville, Sur une ´equation diff´erentielle du premier ordre, Acta mathematica (1902)55–78.[22] A. Llor, Bulk turbulent transport and structure in Rayleigh–Taylor, Richtmyer–Meshkov, andvariable acceleration instabilities, Laser and particle beams (2003) 305–310.[23] A. Llor, Analytical “0D” evaluation criteria, and comparison of single-and two-phase flowapproaches, Statistical hydrodynamic models for developed mixing instability flows , Lecturenotes in physics (Springer, Berlin, 2005).[24] F. Louvet et S. Bouquet, L’´equation “pouss´ee-traˆın´ee” avec acc´el´eration variable avec le temps,Rapport interne CEA-DIF (2010).[25] G.M. Murphy,
Ordinary differential equations and their solutions (Van Nostrand, Princeton,1960).[26] V.E. Neuvazhaev, Theory of turbulent mixing, Soviet physics doklady (6) (1975) 398–400[27] V.E. Neuvazhaev, Properties of a model for the turbulent mixing of the boundary betweenaccelerated liquids differing in density, Journal applied mechanics technical physics (5)(1983) 680–687.[28] P.J. Olver, Applications of Lie groups to differential equations (Springer, Berlin, 1986).[29] L.V. Ovsiannikov,
Group properties of differential equations , (Siberian section of the Academyof Sciences of the USSR, Novosibirsk, 1962) in Russian. Translated by G.W. Bluman (1967),
Group analysis of differential equations (Academic press, New York, 1982).[30] D.E. Panayotounakos and T.I. Zarmpoutis, Construction of exact parametric or closed formsolutions of some unsolvable classes of nonlinear ODEs (Abel’s nonlinear ODEs of the first kindand relative degenerate equations), International journal of mathematics and mathematicalsciences (2011) 387429, 13 pages. doi:10.1155/2011/387429.[31] A. Pandian, N.C. Swisher and S.I. Abarzhi, Deterministic and stochastic dynamics ofRayleigh–Taylor mixing with a power-law time-dependent acceleration, Physica scripta (2017) 014002 (13 pages)[32] A.D. Polyanin, V.F. Zaitsev, Handbook of exact solutions for ordinary differential equations (CRC Press, Boca Raton, 1995).[33] J.D. Ramshaw, Simple model for linear and nonlinear mixing at unstable fluid surfaces withvariable acceleration, Physical review E (5) (1998) 5834–5840.[34] Y. Srebro, Y. Elbaz, O. Sadot, L. Arazi and D. Shvarts, A general buoyancy-drag model forthe evolution of the Rayleigh–Taylor and Richtmyer–Meshkov instability, Laser particle beams21