- Email: [email protected]
The Russo–Smereka kinetic equation: Conservation laws, reductions and numerical solutions
Physica D 303 (2015) 50–58
Contents lists available at ScienceDirect
Physica D journal homepage: www.elsevier.com/locate/physd
The Russo–Smereka kinetic equation: Conservation laws, reductions and numerical solutions Alexander A. Chesnokov a,b,∗ , Maxim V. Pavlov c,d a
Lavrentyev Institute of Hydrodynamics, 15 Lavrentyev Ave., Novosibirsk 630090, Russia
b
Novosibirsk State University, 2 Pirogova Str., Novosibirsk 630090, Russia
c
Lebedev Physical Institute of Russian Academy of Sciences, Leninskij Prospekt 53, 119991 Moscow, Russia
d
National Research Nuclear University MEPHI, Kashirskoe Shosse 31, 115409 Moscow, Russia
highlights • • • •
The Russo–Smereka equation is associated with an integrable hydrodynamic chain. A stability of bubbly flow in terms of hyperbolicity of the model is studied. Infinitely many particular solutions for the kinetic equation are found. A finite-dimensional approximation for the model is proposed.
article
info
Article history: Received 2 February 2014 Received in revised form 30 March 2015 Accepted 31 March 2015 Available online 11 April 2015 Communicated by J. Bronski
abstract The one-dimension Russo–Smereka kinetic equation describing the propagation of nonlinear concentration waves in a rarefied bubbly fluid is considered. Stability of the bubbly flow in terms of hyperbolicity of the kinetic equation is studied. It is proved that a hydrodynamical chain associated with the Russo–Smereka kinetic equation possesses infinitely many conservation laws. Reductions of the model to finite component systems are derived. Conservative form of the kinetic model is proposed and numerical solution of the Cauchy problem with discontinuous initial data is obtained. © 2015 Elsevier B.V. All rights reserved.
Keywords: Kinetic equation Bubbly flow Hyperbolicity Conservation laws Hydrodynamic chain
1. Introduction A kinetic theory based on the statistical description of the interaction of a large number of bubbles has been developed for modelling of nonlinear waves in a rarefied bubbly flow [1–3]. The kinetic models taking into account the effect of collective interaction between bubbles are derived using the system of Hamilton’s ODEs describing the motion of individual bubbles. To obtain this system of equations one needs to know the kinetic energy of the fluid [4]. Assuming that all bubbles are rigid massless spheres
∗ Corresponding author at: Lavrentyev Institute of Hydrodynamics, 15 Lavrentyev Ave., Novosibirsk 630090, Russia. E-mail addresses: [email protected] (A.A. Chesnokov), [email protected] (M.V. Pavlov). http://dx.doi.org/10.1016/j.physd.2015.03.013 0167-2789/© 2015 Elsevier B.V. All rights reserved.
of the same radius Russo and Smereka [1] approximately calculated the energy and Hamiltonian of a bubble motion and obtained the kinetic equation for the evolution of the one-particle distribution function. This model is analogous to the Vlasov equation for plasma flow and describes the collective behaviour of a large number of bubbles subject to long-range interactions, modelled by selfconsistent field. The characteristic properties of the Russo–Smereka kinetic equation for the case of one space variable are studied by Teshukov [5] on the basis of a generalized theory of characteristics and hyperbolicity of integro-differential equations [6,7]. In [5] hyperbolicity conditions of the model are formulated and Riemann invariants and infinite series of conservation laws are found. The exact solutions of the Russo–Smereka kinetic equation in the classes of travelling and simple waves, as well as solutions with linearly dependent Riemann invariants, are obtained and studied in [8,9].
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
In the paper, we study the stability of solutions in terms of hyperbolicity conditions for the Russo–Smereka kinetic equation for bubbly flow considered as an infinite-dimensional quasilinear system. It is known that any solution of the linearized Vlasov equation is stable if it is defined by a distribution function with a single maximum. We show that it is not true for the Russo–Smereka model. We propose a conservative form of the kinetic model for consideration of discontinuous solutions. Differential conservation laws approximating the Russo–Smereka kinetic equation are derived. These laws are used to perform numerical calculations of wave propagation in a rarefied bubbly flow showing the possibility of the kinetic roll-over (the formation of two peaks of the distribution function which originally had a single peak). This effect usually leads to instability of the flow. We consider a hydrodynamical chain associated with the Russo–Smereka kinetic equation and show that this chain is integrable. Some new reductions of the kinetic equation to finite component systems are also obtained and their properties are studied. The paper is organized as follows. In Section 2 we present different formulations of the Russo–Smereka kinetic model. Each formulation is more suitable for some specific application. We also recall hyperbolicity conditions for this integro-differential model and give an example of verification of the hyperbolicity conditions. Conservative formulation of the kinetic model for a bubbly flow is proposed in Section 3, as well as its finite-dimensional approximation in the form of a system of differential conservation laws with a large number of unknown variables. In Section 4 we introduce a hydrodynamic chain associated with the Russo–Smereka kinetic equation. We show that this chain possesses infinitely many conservation laws. In Section 5 we apply the method of hydrodynamic reductions for constructing infinitely many particular solutions of the Russo–Smereka kinetic model. Then we obtain finitecomponent reductions using generalized and step-functions. Numerical results are presented in Section 6. We show that during the evolution the kinetic roll-over of the distribution function is possible. Finally we draw some conclusions. 2. Mathematical model
51
Let us show that Eqs. (4) are derived as a consequence of Russo–Smereka model (1). We use f˜ to denote the distribution function in semi-Lagrangian coordinates: f˜ (t , x, λ) = f (t , x, p(t , x, λ)). Then, the derivatives are represented as
∂f ∂f ∂p ∂ f˜ = + , ∂t ∂t ∂p ∂t
∂ f˜ ∂f ∂f ∂p = + . ∂x ∂x ∂p ∂x
These formulae and Eqs. (3) imply the obvious equalities
∂ f˜ ∂ f˜ ∂f ∂f ∂j ∂f + (p − j) = + (p − j) + p = 0, ∂t ∂x ∂t ∂x ∂x ∂p which lead to the second equation of system (4). In what follows, we assume that the inequality pλ > 0 is satisfied. It provides invertibility of the change of variables. 2.1. Hyperbolicity conditions for the kinetic equation Eqs. (4) belong to the class of systems with operator coefficients of the form Ut + A⟨Ux ⟩ = 0,
(5)
for which a generalization of the hyperbolicity notion was proposed in [6]. Here, U(t , x, λ) is an unknown vector function and A is a non-local operator on a set of functions of λ. Characteristics of system (5) are determined by equation x′ (t ) = k(t , x), where k is an eigenvalue of problem (F, (A − kI )⟨ϕ⟩) = 0. The eigenfunctional F is defined on a set of functions of variable λ, while the values of t and x are considered to be fixed, and is sought in a class of locally integrable or generalized functions (ϕ(λ) is a smooth test function). System of Eqs. (5) is hyperbolic if all eigenvalues k are real and the set of relations on the characteristics (F, Ut + kUx ) = 0 is equivalent to (5). As it follows from [5,7], hyperbolicity conditions for Russo– Smereka equation (1) on a solution f (t , x, p) are formulated in terms of characteristic function
χ (k + j) = 1 − n + (k + j)
2
In the one-dimensional case the Russo–Smereka kinetic equation in dimensionless variables is written as follows [1,5]:
∂f ∂f ∂j ∂f + (p − j) + p = 0, ∂t ∂x ∂x ∂p j(t , x) = pf (t , x, p) dp.
(1)
Here f (t , x, p) ≥ 0 is the distribution function for bubbles in phase space, x and p are the position and momentum of the bubble, t is the time, and j is the first moment of the distribution function. It is assumed that f decreases rapidly at infinity or has compact support in p. The model is suitable for the description of a rarefied bubbly flow in the case of small pressure variations. The condition for a bubbly flow to be rarefied is given by the inequality n(t , x) =
f dp < 1.
(2)
It was shown in [5], to study the properties of kinetic equation (1), it is appropriate to make use of the Eulerian–Lagrangian coordinates x and λ by substitution the variable p = p(t , x, λ), where the function p(t , x, λ) is a solution of the Cauchy problem pt + (p − j)px = pjx ,
p(0, x, λ) = p0 (x, λ),
λ ∈ [0, 1]. (3) As a result, for the new functions p(t , x, λ) and f (t , x, λ), we obtain
f dp
(p − k − j)2
,
(6)
or, more precisely, in terms of its limit values on the real axis from upper and lower complex half-planes
χ ± (p) = 1 − n(t , x) + p2 ± π ip2
∂ f (t , x, p′ ) dp′ ∂ p′ p′ − p
∂ f (t , x, p) . ∂p
(7)
The above formula is obtained from (6) by integration by parts and application of Sokhotski–Plemelj formulae. According to [5,7], kinetic equation (1) is hyperbolic on the rapidly decreasing solution f (t , x, p) if the following conditions are valid
1 arg χ ± (p) = 0.
(8)
The argument increment is calculated when p changes from −∞ to ∞ at fixed values of variables t and x. If supp f is bounded, then conditions (8) take the form
1 arg χ + (p)/χ − (p) = 0, here the argument increment is calculated for p ∈ supp f . Hyperbolicity conditions (8) guarantee that characteristic equation χ (k + j) = 0 has no complex roots, and these conditions are necessary for the flow stability.
the integro-differential system of equations [5,7] 2.2. Example of verification of the hyperbolicity conditions
pt + (p − j)px − pjx = 0, ft + (p − j)fx = 0,
j=
ppλ f dλ.
(4)
Conditions (8) allow one to verify whether Russo–Smereka equation (1) is hyperbolic for given solution f (t , x, p). Following
52
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
satisfied at the single interior point of p = pc at which the distribution function f (p) reaches a local maximum. Since the distributions f = f (p) with one maximum obey the inequality (pc − p)f ′ (p) ≥ 0, one can conclude that the sign of the quantity Z 1 (pc ) = 1 − n + p2c
Fig. 1. Contours C + in the complex plane (Z 1 , Z 2 ) (the solid curve corresponds to a = 1/2; the dashed curve to a = 1).
f ′ (p)dp p − pc
depends on the point pc . Bubbly flow is stable if the extreme point pc is sufficiently close to zero (this guarantees that the inequality Z 1 (pc ) > 0 is fulfilled). Back to the previous example (Figs. 1 √and 2), one can conclude that the flow is stable for |a| < ac = 1/ 2 (unstable for |a| > ac ), where a = ac is the root of the equation Z 1 (a) = 0. 3. Conservative form of the model The evolution of a smooth solution of the hyperbolic system of equations can involve a gradient catastrophe. Therefore, further description of the solutions is possible only for the class of discontinuous functions. It leads to the necessity to formulate the model in a form of conservation laws. To do this, we present (4) in the form Ht + (p − j)H
x
jt + A2 − 3j /2
Fig. 2. Distribution function f (p) with one maximum (notation is the same as in Fig. 1).
[7,10], in the plane (Z 1 , Z 2 ) we construct a closed contour C consisting of the contours C + and C − . The contour C + is given parametrically by the equations Z 1 = Re{χ + (p)},
Z 2 = Im{χ + (p)},
where the complex functions χ ± (p) are defined by (7). A contour C − , which is symmetric about the Z 1 axis to the contour C + , is given by the same equation with the function χ − (p). If the point of Z 1 = 0, Z 2 = 0 lies in the domain bounded by the contour C , then the characteristic equation χ (k) = 0 has complex roots (the function χ (k) is given by (6)). Otherwise, the kinetic equation for the corresponding solution is hyperbolic. Fig. 1 shows the contours C + (the direction of circulation about the contours is positive, i.e. counter-clockwise) corresponding to the distribution function with one maximum f (p) =
1 √ exp −(p − a)2 ,
2 π
p ∈ (−∞, ∞).
This function is shown in Fig. 2 for a = 1/2 (solid curve) and a = 1 (dashed curve). As it can be seen from the graphs, hyperbolicity conditions (8) are fulfilled if the extreme point p = a is closely spaced to zero (solid curve in Fig. 2). Increasing the parameter a leads to an increment in the argument of the functions χ ± (dashed line in Fig. 1). It corresponds to existence of the complex characteristic roots and, consequently, the flow is unstable. In the theory of plasma waves the following result is known [11]: any solution of the one-dimensional linearized Vlasov equation is stable if it is defined by a distribution function with a single maximum. We show that in the kinetic theory for a bubbly flow this is not true. It is mentioned above, hyperbolicity conditions (8) are violated if the point Z 1 = 0, Z 2 = 0 is in the domain bounded 1 by the contour C . Let Z 1 → Z±∞ > 0 for p → ±∞. Then, this is possible only if the conditions Z 1 (p∗ ) < 0 and Z 2 (p∗ ) = 0 are satisfied at some point p∗ ∈ (−∞, ∞). The equality Z 2 (p∗ ) = 0 is
j=
2
= 0,
pλt + (p − j)pλ
x
= 0,
=0 pH dλ, A2 = p2 H dλ, H = pλ f . x
(9)
The first condition in (9) is the local conservation of a number of bubbles, the second condition is equivalent to the conservation of the function f along the trajectory, and the last equation is the conservation law for hydrodynamics momentum j. Note that Eqs. (9) are similar to conservation laws of the shallow water equations for shear flows [12,13]. As a rule, in the case of infinitedimensional systems to substantiate the correctness of the choice of conservation laws is difficult. In what follows, we consider solutions with small amplitude jumps. From conservation laws (9) we obtain the Hugoniot conditions at the shock front x = x(t ) moving with the velocity V = x′ (t ):
(p − j − V )H = 0, (p − j − V )pλ = 0, A2 − 3j2 /2 − Vj = 0.
(10)
Here [ϕ(t , x, λ)] = ϕ(t , x(t ) + 0, λ) − ϕ(t , x(t ) − 0, λ) is a jump of the function ϕ at the shock front. Calculating the ratio of continuous quantities (p − j − V )H and (p − j − V )pλ , we obtain the following jump relation: [f ] = 0. Let us show that systems (4) and (9) are equivalent for smooth solutions. System (9) is a consequence of Eqs. (4) obtained by simple transformations. Indeed differentiating with respect to λ the first equation in (4) gives the second one in (9). The following equality leads to the first equation in (9): 0 = pλt + ((p − j)pλ )x f + ft + (p − j)fx pλ
= Ht + (p − j)H x . Finally we note that
(pH )t + ((p − j)pH )x − pHjx = 0. Integrating this equation with respect to λ, we get the third equation in (9). Now we show that (4) is a consequence of (9). From the first two equations (9) we have ft + (p − j)fx = 0. Integrating with respect to λ the second equation in (9) gives the equality pt + (p − j)px − pjx = G(t , x),
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
53
where G(t , x) is an arbitrary function. Let us multiply this equation by H and integrate with respect to λ. Taking into account the equality
(pt + (p − j)px )H = (pH )t + ((p − j)pH )x , which holds by virtue of (9), after integration we have jt + A2 − 3j2 /2
x
= nG.
The third equation in (9) leads to the identity G ≡ 0 which proves equivalence of (4) and (9) for smooth solutions. 3.1. Approximate model Fig. 3. Graph of the function χ( ¯ k) for M = 5.
To derive differential conservation laws that approximate integro-differential model (9) we divide the segment λ ∈ [0, 1] into M intervals 0 = λ0 < λ1 < · · · < λM −1 < λM = 1 and introduce the variables pi + pi−1 , pi = p(t , x, λi ) hi = pi − pi−1 , p¯ i = 2 λi 1 H dλ. f¯i = hi λi−1 Taking into account the equality H dλ = f dp and using a piecewise constant approximation of the distribution function f (t , x, p) = f¯i (t , x), λi λi−1
A2 =
p ∈ [pi−1 , pi ]
i =1
pi
p2 f dp =
pi−1
M ¯ 3 fi h i
12
i =1
2
+
× j+
hk +
k=1
M
(11)
f¯i hi
−1
i=1 M
1 2 i =1
f¯i h2i −
M
1 qi − k
f¯i = 0
(13)
f¯i hi ,
j=
M
p¯ i f¯i hi .
i =1
i=1
f¯M < f¯M −1 < · · · < f¯l ,
that correspond to an approximation of the distribution function with a single maximum. Then the equation χ¯ (k) = 0 has a unique root k = ki on each interval k ∈ (qi−1 , qi ), i = 1, . . . , l − 1, l + 1, . . . , M. Indeed, function χ¯ (k) is monotonic and varies from −∞ to ∞ on the intervals (qi−1 , qi ) for i = 1, . . . , l − 1, and from ∞ to −∞ for i = l + 1, . . . , M. Consequently, system (11) is hyperbolic if Eq. (13) has two roots in the interval (ql−1 , ql ). Since χ (−j) = 1 − n > 0 and χ¯ (k) tends to −∞ as k → ql−1 (or k → ql ), then the inequalities
(ql−1 < −j < ql )
are the sufficient conditions for hyperbolicity of system (11). The typical form of the function χ¯ (k) for M = 5 is shown in Fig. 3. This dependence is obtained for the piecewise constant approximation of the distribution function with bounded support: f (p) = (cos p + 1)/8, p ∈ (−π , π ). Note that integro-differential model (4) has continuous characteristic spectrum k = p − j, where p ∈ supp f . 4. Local conservation laws
f¯i hi
i
hk .
Introducing infinitely many moments
k=1
To solve differential conservation laws (11) numerically, one can apply Godunov type methods. In this case, due to a large number of equations in system (11), it is convenient to use central schemes [14], which do not require exact or approximate solution of the Riemann problem. System (11) can be written in the following form ut + A(u)ux = 0,
M
pl−1 < 0 < pl
The quantities pci (t , x) included in (11) are given by the formulae p¯ i = −
qi−1 − k
−
where χ¯ (k) is a discrete analogue of the characteristic function (6), qi = pi − j and the moments n and j are given by the formulae
f¯1 < f¯2 < · · · < f¯l ,
+ f¯i hi p¯ 2i .
∂ ∂ hi + (¯pi − j)hi = 0, ∂t ∂x ∂ ¯ ∂ fi hi + p¯ i − j f¯i hi = 0, ∂t ∂x M ¯ 3 3 f i hi ∂ ∂j 2 ¯ + + fi hi p¯ i − j2 = 0. ∂t ∂ x i=1 12 2 i
1
In the hyperbolic case Eq. (13) has M + 1 real roots ki (i = 1, . . . , M + 1) in the interval (p0 − j, pM − j). Let the values f¯i satisfy the inequalities
Further, we integrate the first two equations in (9) with respect to λ over the intervals (λi−1 , λi ) and use the previous formulae. As a result we obtain a system of conservation laws consisting of 2M + 1 differential equations for unknown functions hi (t , x), f¯i (t , x), and j(t , x):
hi
M i=1
i =1
pH dλ = f¯i hi p¯ i , M
χ¯ (k) = 1 − n + (k + j)2
n=
we have the equalities
f¯it + (¯pi − j)f¯ix = 0. This means that the matrix A(u) has the following eigenvalues: k∗i = p¯ i − j (i = 1, . . . , M ). To find the other eigenvalues of A(u), one can use the equation
(12)
where u = (h1 , . . . , hM , f¯1 , . . . , f¯M , j) and A(u) is the corresponding matrix. This form is used to compute the characteristic values. As follows from (11), the functions f¯i satisfy the equations
Ak =
pk f dp (k = 1, 2, . . .)
kinetic model (1) leads to the modified Benney hydrodynamic chain
∂ Ak ∂ Ak+1 ∂ Ak ∂ A1 + − A1 − (k + 1)Ak = 0, ∂t ∂x ∂x ∂x
(14)
which was first derived by Kupershmidt (see details in [15]) in framework of the KP theory. This hydrodynamic chain is integrable. This means that (14) possesses infinitely many conservation laws (18), commuting flows and local Hamiltonian structures.
54
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
Suppose that some particular solution of (14) is found, i.e. Ak (t , x) are determined. Then the distribution function f (t , x, p) can be reconstructed. Indeed in this case one can look for solution in the form f (t , x, p) = F (p, A1 , A2 , . . .). Substituting this ansatz directly into (1) yields infinitely many equations which can be integrated (see details in [16,17]). This dependence can be presented in the Laurent series form (p → ∞) F =
1 p
+
A1 p2
+
A2 p3
+
A3 p4
+ ···.
(15)
Substitution of expansion (15) into kinetic model (1) implies Eqs. (14) again. As we have already mentioned modified Benney chain (14) possesses infinitely many local conservation laws. A generating function of these conservation laws can be obtained from kinetic model (1). The first equation in (1) divided by fp takes the form ft fp
+ (p − j)
fx
+ pjx = 0.
fp
Then taking into account (according to the theorem about differentiation of an implicit function) that pt = −ft /fp and px = −fx /fp , one can obtain
pt +
p2 2
1
= 0,
− jp
(16)
x
− P1 − P2 F − P3 F 2 − · · · ,
(17)
F and all conservation law densities Pk are the polynomials with respect to the moments Am . Explicit expressions can be obtained by substituting expansion (15) into (17) and equating factors of every degree of p. For instance, P1 = A1 , P2 = A2 − A21 , P3 = A3 − 3A1 A2 + 2A31 . Then substituting (17) into (16) and equating to zero the expressions of the same degrees of F leads to the following conservation laws
∂ Qk ∂ Pk + = 0, ∂t ∂x
(k = 1, 2, . . .)
(18)
where the first three polynomials Qk are Q1 = A2 − 3A22 /2, Q2 = A3 − 3A1 A2 + 2A31 , Q3 = A4 − 4A1 A3 − A22 + 8A21 A2 − 4A41 . The existence of an infinite number of conservation laws is a rare property for hydrodynamics models. Benney system of the long-wave theory is a well-known example of the model having an infinite number of conservation laws [18]. Both models, Russo–Smereka and Benney, are generalized hyperbolic in the sense of [6] and can be written in terms of the Riemann invariants conserving along the characteristics. 5. Method of hydrodynamic reductions The method of hydrodynamic reductions was established in [19] for Vlasov (collisionless Boltzmann) equation. In the case of Russo–Smereka kinetic model (1) the same approach is applicable. The method of hydrodynamic reductions is based on a natural observation: all moments Ak (t , x) of hydrodynamic chain (14) can be parameterized by N arbitrary functions am (t , x) only, then hydrodynamic chain (14) becomes N component hydrodynamic type system, where N is an arbitrary natural number. Example. The moment decomposition (εm are arbitrary constants) Ak =
1
N
k + 1 m=1
akt +
εm (am )k+1
(19)
(ak )2 2
− ak A1 (a) = 0,
(20)
x
where a = (a1 , . . . , aN ), the function A1 (a) is given by (19) for k = 1. It should be noted that system (19), (20) is not a sole reduction of hydrodynamic chain (14). Indeed one can compare (20) with (16). Thus reduction (20) follows directly from (16) replacing p by ak (t , x). Let us assume that all moments Ak depend on a finite number of field variables am (t , x) satisfying hydrodynamic type system (20). All equations (infinitely many) of hydrodynamic chain (14) must be consistent with each other. Without loss of generality (see details in [16,17]) we use the first equation from hydrodynamic chain (14) only (we assume summation with respect to the repeated index):
∂ A1 (a) k ∂ A1 (a) k ∂ A2 (a) k a + a − 3A1 (a) a = 0. ∂ ak t ∂ ak x ∂ ak x Taking into account (20) we obtain the following relations
where a generating function of conservation law densities p(F , A1 , A2 , . . .) can be found from (15), i.e. p=
transforms hydrodynamic chain (14) into the so-called ‘‘waterbag’’ reduction
∂ A1 (a) ∂ A1 (a) ∂ A1 (a) − 2A1 (a) m k ∂a ∂a ∂ ak ∂ A1 (a) ∂ A2 (a) k − ak + ax = 0, ∂ ak ∂ ak
am
where each bracket must vanish independently for each index k. The compatibility conditions
∂ ∂ Al ∂ ∂ Al = k i i k ∂a ∂a ∂a ∂a for each pair of distinct indices i and k lead to the overdetermined system for A1 (a(t , x))
( ai − ak )
∂ 2 A1 ∂ A1 ∂ A1 ∂ A1 ∂ A1 + k Rˆ i − i Rˆ k = 0, ∂ ai ∂ ak ∂a ∂a ∂a ∂a
i ̸= k.
(21)
Here Rˆ = Σ am ∂/∂ am is the scaling operator. A direct computation shows that all other higher equations of hydrodynamic chain (14) allow to reconstruct dependences Ak (a) in quadratures (we omit here this technical verification) if the function A1 (a) is a solution of overdetermined system (21). Since i ̸= k all second derivatives are expressible via the first derivatives except the second derivatives with coincided indices. Thus system (21) possesses a general solution parameterized by N arbitrary functions. A similar system for description of N component reductions for the remarkable Benney hydrodynamic chain was derived by Gibbons and Tsarev in [19]. By this reason we call Eqs. (21) the Gibbons–Tsarev system. In a separate paper [20] we proved that chain (14) is nothing but the so-called modified Benney hydrodynamic chain. This means that both chains have the same amount of local conservation laws. Some multi-parametric solutions of Gibbons–Tsarev system (21) were also found in [20]. Instead of substitution (19), one can choose any other particular solution of overdetermined system (21). One can see, for instance, that a family of particular solutions can be extracted from (21) by the natural ansatz A1 (a) =
N
gm (am ),
m=1
where functions gm (am ) can be found by direct substitution into (21). Indeed in this case overdetermined system becomes
ˆ i′ (ai ) = gi′ (ai )Rg ˆ k′ (ak ), gk′ (ak )Rg
i ̸= k.
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
ˆ i′ (ai ) = ai gi′′ (ai ) one can obtain Taking into account that Rg ak gk′′
where km are the roots of characteristic equation (13). This representation allows one to construct exact solutions in the form
(a ) = α gk (a ), ′
k
55
k
ri = ri0 = const (i ̸= s),
where α is an arbitrary constant. This means that gk′ (ak ) = βk (ak )α , where βk are extra arbitrary (integration) constants. Below for simplicity we restrict our consideration on two main interesting and important cases determined by the Dirac deltafunctional and Heaviside step-functional ansatz. We also recall special class of solutions [9] characterized by a linear relationship between the Riemann invariants. This class of solutions will be used for testing the numerical results in the framework of approximate model (11).
where f = f (t , x, λ) and
5.1. Generalized functions and reductions
R(t , x, λ) =
Under certain assumptions kinetic model (1) can be reduced to the system of differential equations. The following representation of the solution
rs = rs (k),
k = ks + j.
These relations define the functions p0 (k), . . . , pM (k), where k(t , x) is a solution of the equation kt + kkx = 0. The existence of Riemann invariants for system (22) follows from the fact that integro-differential equations (4) can be presented in the characteristic form [5] ft + (p − j)fx = 0,
Rt + (p − j)Rx = 0,
n(t , x) − 1 p(t , x, λ)
f (t , x, ν)pν dν
+
p(t , x, ν) − p(t , x, λ)
(23)
are the Riemann invariants (here n is given by formula (2)). 5.2. Special class of solutions
M
f =
ni (t , x)δ(p − pi (t , x))
We also consider solutions of Eq. (1) in the class of functions that are piecewise continuous in the variable p with a bounded support:
i=1
leads to the ‘‘frozen’’ bubbly flow
f = f 1 (t , x, p) θ (p − pl (t , x)) − θ (p − pr (t , x)) .
∂ ni ∂ + (pi − j)ni = 0, ∂t ∂x M ∂ pi ∂ p2i + − pi j = 0, j = ni p i . ∂t ∂x 2 i=1
Here θ is a Heaviside step function, pl and pr are the boundaries of the interval in the variable p, beyond which the distribution function is identically zero, and f 1 (t , x, p) ≥ 0 is a non-negative function which is continuously differentiable on the set {(t , x, p) : t ≥ 0, x ∈ R, p ∈ [pl , pr ]}. Substitution of this ansatz into (1) yields
Here δ(p) is a Dirac delta function. This system of equations has imaginary characteristic roots [21]. The concept of waterbag is known in plasma physics [22]. It allows one to obtain particular solutions of kinetic equation (1) described by closed system of equations. We represent the kth moment in the form [16]
∂f ∂f ∂j ∂f + (p − j) + p = 0, ∂t ∂x ∂x ∂p pr j(t , x) = pf (t , x, p) dp,
Ak =
1
M
k + 1 i=0
M
εi pik+1 ,
εi = 0,
i =0
where εi are arbitrary constants whose sum is equal to zero. Substitution of the above introduced moments into chain (14) yields a closed system of M + 1 equations for unknown functions pi (t , x):
∂ pi ∂ pi ∂j + (pi − j) − pi = 0, ∂t ∂x ∂x
j=
M 1
2 i =0
εi p2i .
(22)
This system coincides with conservation laws (11) only for the class of solutions f¯i = fi = const (in this case (11) is the exact consequence of (1)). Note that due to the following representation of the solution: M
f =
M
fi θ (p − pi−1 ) − θ (p − pi ) ≡
i =1
i=0
εi θ (p − pi ),
i =0
M
i
εi = 0,
fi+1 =
εi
k=0
system (22) can be obtained directly from kinetic model (1) (here θ is a Heaviside step function). This system of equations is hyperbolic under some conditions derived above. Waterbag reduction (22) can be written in Riemann invariants
∂ rm ∂ rm + (km + j) = 0 (m = 1, . . . , M + 1), ∂t ∂x M p −k −j n−1 m i rm = + fi ln , km + j p − k − j i − 1 m i =1
pl
(24)
∂ pl ∂ pl ∂j + (pl − j) − pl = 0, ∂t ∂x ∂x ∂ pr ∂ pr ∂j + (pr − j) − pr = 0, ∂t ∂x ∂x where f = f 1 .
The class of solutions of kinetic model (1) with a bounded support characterized by a linear relationship between the Riemann integral invariants was obtained in [9] on the basis of the following property: if the functions f (t , x, p), pl (t , x), and pr (t , x) are a solution of system (24) then the function R(t , x, p), which has been defined above by (23) in semi-Lagrangian coordinates, satisfies the equation Rt + (p − j)Rx + pjx Rp = 0. Linear relationship between the invariants R and f R + f π cot(µπ ) − b = 0
(25)
(µ and b are constants) leads to a special class of solutions [9]. In this case the distribution function f has the form f =
sin(µπ )
1 − µb(pr − pl )
π
p
+b
p − pl
µ
pr − p
,
(26)
and the functions pl and pr satisfy system (24) with the following first moment
j(pl , pr ) = µ(pr − pl ) 1 +
1+µ 2
bpl +
1−µ 2
bpr .
(27)
Description of discontinuous solutions from this special class is given by the following closed system of conservation laws nt + (1 − n)j
x
= 0,
jt + A2 (n, j) − 3j2 /2
x
= 0.
(28)
56
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
By virtue of (26) the moments n, j, and A2 are expressed in terms of pl and pr by (27) and formulae
n=1−
pl
µ
pr
A2 = µ(pr − pl )
+
1 − µb(pr − pl ) ,
1 − µ2 3
1+µ 2
pr +
1−µ 2
6. Numerical results
pl
b(pr − pl ) + bpl pr . 2
Relations n = n(pl , pr ) and j = j(pl , pr ) can be inverted if the Jacobian determinant D(n, j)
µ
pl µ = −(pr − pl ) 1 + pl − µ(pr − pl ) b D(pl , pr ) pl pr pr × 1 + pr − µ(pr − pl ) b is non-zero. In the case b = 0 this condition reads pl ̸= pr (pl and pr det
2
are not equal to zero). Below these conservation laws will be used to perform numerical simulation of discontinuous bubbly flows from special class of solution as well as to test results obtained on the basis of general model (11). Remark. In general relation (25) is not fulfilled at the shock front. This means we cannot use conservative form (28) for the description of discontinuous flows. Only the general system of conservation laws (11) can be applied for this purpose. Nevertheless, modelling of discontinuous flows with small amplitude jumps δ = [n] is possible on the basis of simplified Eqs. (28). Following [7,12], one can show that [R] = O(δ 2 ) and, consequently, Hugoniot conditions
(p − j − V )H = 0, R + f π cot(µπ ) − b = 0, A2 − 3j2 /2 − Vj = 0, which are performed for a special class of solutions, coincide with relations (10) if the second-order terms in δ are neglected. Notice that the previous jump conditions are consequence of the following conservative form Ht + (p − j)H
= 0, jt + (A2 − 3j2 /2)x = 0, x (R + f π cot(µπ ) − b)H t + (R + f π cot(µπ ) − b)(p − j)H x = 0
n+1/2
uj
= unj − Λf′j /2,
(Λ = 1t /1x)
ujn++11/2 = (unj + unj+1 )/2 + (u′j − u′j+1 )/8
− Λ f(
n+1/2 uj+1
n+1/2 uj
) − f(
(30)
) .
This scheme approximates the system of conservation laws of the form ut + (f(u))x = 0 including models (11) and (28). Here 1x is the spatial grid spacing, while 1t is the time-step satisfying the Courant condition, and unj = u(t n , xj ). The computational domain on the x axis is divided into N cells, the cell centres are denoted by xj . Values u′j /1x and f′j /1x are approximations of the first-order derivatives with respect to x, calculated according to the ‘‘ENO limiter’’ procedure [23]. At t = 0 the initial data u0j are specified. The boundary conditions u(t n , x1−j ) = u(t n , x1 ) and u(t n , xN +j ) = u(t n , xN ) are used, which allow calculations to be performed until the initial perturbations reach the boundaries of the computation domain. Scheme (30) does not require exact or approximate solution of the Riemann problem that is very convenient in our case, due to a large number of equations in system (11). It should be noted that central schemes on non-staggered grids, such as Kurganov–Tadmor scheme [24] or, in general, schemes based on the local Lax–Friedrichs flux are also appropriate. Such non-staggered approach has the advantage that it is usually simpler to assign the boundary conditions. However, in the cases under study we consider trivial boundary conditions. Therefore, the Nessyahu–Tadmor scheme is perfectly appropriated.
Let us perform a comparison between the numerical solutions obtained for differential approximation (11) and for system (28) defining solutions from this special class. At t = 0 we define the function f by formula (26), where pl (x) = 0.5 and pr (x) = 2.25 if x ∈ (−0.5, 0), otherwise pl = 0.25 and pr = 2. Initial distribution function f (x, p) and moments n(x) and j(x) are shown in Figs. 4 and 5 by dashed lines. We calculate the solution in the domain x ∈ [−1, 1] for N = 150 and M = 120. To start calculations on the base of model (11) we specify unknown variables
5.3. Fluid dynamic limit Let the distribution function have the form
(p − p¯ )2 n f = √ exp − . 2T 2π T Calculation of moments of this distribution function A1 = j = np¯ ,
In this section we present some results related to the numerical modelling of bubbly flows. We implement here the Nessyahu– Tadmor second-order central scheme [14]
6.1. Comparison with a special class of solution
of system (4).
A0 = n,
Taking T = h21 /12, p¯ = p¯ 1 , and n = f¯1 h1 it is easy to see that system (29) coincides with the differential approximation (11) for M = 1. Note that the eigenvalues k1,2 are the roots of characteristic equation (13).
A2 = (¯p2 + T )n,
A3 = (¯p2 + 3T )np¯
(h1 , . . . , hM , f¯1 h1 , . . . , f¯M hM , j)
and their substitution into the first three conservation laws (18) leads to the closed system of equations for unknown functions n(t , x), p¯ (t , x) and T (t , x):
at t = 0. Here hi = h ≡ (pr − pl )/M; f¯i is the mean value of the function f (0, x, p) in the interval (pi−1 , pi ), where pi = pl + ih (this condition is imposed at the initial time). The following parameters for the special class of solutions are chosen µ = 0.6, b = 0. A sufficiently detailed resolution in the variable p is needed due to the rapid change of the distribution function in a special class of solutions in the neighbourhood of the points pl and pr . Figs. 4 and 5 (solid lines) show the results of computation at t = 1. This test confirms that the calculations of discontinuous solutions obtained by ‘‘multilayer’’ approximation (11) and by the system for special solutions (28) give close results, at least for small amplitude jumps.
nt + (¯un)x = 0,
(np¯ )t + (nT + np¯ u¯ − n2 p¯ 2 /2)x = 0, (nT + np¯ u¯ )t + (3T + (1 − 2n)¯p2 )nu¯ x = 0,
(29)
where u¯ = (1 − n)¯p. This system of equations, obtained in [21], is hyperbolic for T > np¯ 2 /3. In fact, let us rewrite (29) in the form (12), where u = (n, p¯ , T ). The eigenvalues of A(u) are k0 = (1 − n)¯p,
k1,2 = (1 − 2n)¯p ±
(1 − n)(3T − np¯ 2 ).
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
Fig. 4. Distribution function f at t = 0 (curves 1 and 2) and t = 1 (curves 3 and 4): 1—x ∈ (−0.5, 0); 2—x ̸∈ (−0.5, 0); curves 3 and 4 refer to calculation using models (28) and (11) at x = −0.2.
57
Fig. 6. Distribution function f at t = 0 (curve 1—x < 0, curve 2—x > 0) and at t = 2 (curve 3—x = −0.5, curve 4—x = 0.5).
Fig. 7. Moments n and j at t = 0 (curves 1, 2) and at t = 2 (curves 3 and 4). Fig. 5. Density n (curves 1–3) and moment j (curves 4–6): 1 and 4—t = 0; 2 and 5—t = 1 calculation using model (28); 3 and 6—t = 1 calculation using model (11).
6.2. Kinetic roll-over In the theory of quasineutral collisionless plasma flows the following result is known [10,25]: during the evolution the kinetic roll-over of the distribution function is possible (the formation of two peaks of the distribution function which originally had a single peak). Let us establish a similar property for the considered kinetic model for a bubbly flow. Suppose that at t = 0 the bubbly flow in the half-space x < 0 is defined by a distribution function f = fl (p) and in the halfspace x > 0 by a distribution function f = fr (p). In the interval p ∈ (pl , pr ) we choose the functions fl and fr as follows: fl (p) =
exp(−p2 )
, √ 5 π erf(1)
fr (p) =
exp(−4p2 )
√
5 π erf(2)
(pr = −pl = 1)
(curves 1 and 2 in Fig. 6); outside this interval fl = fr = 0. Distribution function at t = 2 for x = −0.5 and x = 0.5 is shown in Fig. 6 (curves 3 and 4, correspondingly). It can be seen that the distribution function for fluid more saturated by bubbles (curves 1 and 3) has changed in the minimum of p. At the same time the distribution function of the less saturated fluid (curves 2 and 4) has qualitative changes due to kinetic roll-over. In the vicinity of x = 0 (point of discontinuity at t = 0) the distribution function f for t > 0 has the following form: f → fl (p) for p < 0 and f → fr (p) for p > 0. Note that the formation of two local peaks of the distribution function usually leads to loss of the hyperbolicity of the kinetic model. Fig. 7 shows plots of the hydrodynamic moments at t = 2 (solid curves) and t = 0 (dotted). In the calculations we used the following resolution N = 400 and M = 100. Solution of the
Riemann problem with other initialdata (functions fl (p) and fr (p) have one maximum and pfl dp > pfr dp) has a similar form. Increase or decrease in the number of nodes has no significant effect on the numerical results. This remark is valid for all numerical examples presented in the section. 7. Conclusion Different formulations of Russo–Smereka kinetic equation (1) are presented. Semi-Lagrangian formulation allows one to rewrite this model as infinite-dimensional system of quasilinear equations (4) and propose its conservative form (9). This formulation gives possibility to consider discontinuous solutions. However a question about the physically correct conservative form of the model remains open. For this reason we restrict our consideration for discontinuous solutions with small amplitude jumps. The procedure of averaging with respect to the Lagrangian coordinate allows one to obtain finite-dimensional approximation in the form of differential conservation laws (11). Some known fluid dynamic limits (Eqs. (29) obtained under the assumption of local thermodynamic equilibrium and ‘‘waterbag’’ reduction (22)) are special cases of this system of equations. An example of verification of the hyperbolicity conditions of the kinetic model is given. It is established that the distribution function with one peak (Figs. 1 and 2) can lead to instability for a bubbly flow. This property distinguishes the Russo–Smereka model from the Vlasov equation because any solution of the onedimensional linearized Vlasov equation is stable if it is defined by a distribution function with a single maximum. It is established that Russo–Smereka kinetic model (1) is associated with an integrable hydrodynamic chain (14). It has infinitely many conservation laws
58
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
(18). The method of hydrodynamical reductions is used to construct infinitely many particular solutions for the Russo–Smereka kinetic equation. Conservation laws (11) and (28) are used to perform numerical calculations of wave propagation in a bubbly flow initiated by discontinuous Cauchy data. It is shown the correspondence of the numerical results in the framework of the proposed finitedimensional approximation and in the special class of solutions characterized by a linear relationship between the Riemann integral invariants (Figs. 4 and 5). The effect of the kinetic roll-over of the distribution function is demonstrated (Fig. 6). As a rule it leads to instability of the flow. Acknowledgements The authors thank S.L. Gavrilyuk for stimulating discussions, and an anonymous referee for useful comments and suggestions. AAC’s work was supported by RFBR grant 13-01-00249 and Integration Project of SB RAS No. 30. MVP’s work was partially supported by RF Government grant 010-220-01-077, ag. G34.31.0005, by grant of Presidium of RAS ‘‘Fundamental Problems of Nonlinear Dynamics’’, and by RFBR grant 14-01-00389. References [1] G. Russo, P. Smereka, Kinetic theory for bubble flow I: collisionless case, SIAM J. Appl. Math. 56 (1996) 327–357. [2] H. Herrero, B. Lucquin-Desreux, B. Perthame, On the motion of dispersed bubbles in a potential flow, SIAM J. Appl. Math. 60 (1999) 61–83. [3] V.M. Teshukov, S.L. Gavrilyuk, Kinetic model for the motion of compressible bubbles in perfect fluid, Eur. J. Mech. B Fluids 21 (2002) 469–491. [4] L.M. Milne-Thomson, Theoretical Hydrodynamics, Macmillan, London, 1960. [5] V.M. Teshukov, Characteristics, conservation laws, and symmetries of the kinetic equations of motion of bubbles in a fluid, J. Appl. Mech. Tech. Phys. 40 (1999) 263–275.
[6] V.M. Teshukov, Hyperbolicity of long-wave equations, Dokl. Akad. Nauk 284 (1985) 555–559. [7] V.Yu. Liapidevskii, V.M. Teshukov, Mathematical models for long-wave propagation in an inhomogeneous fluid, Izd. Sib. Otd. Ross. Akad. Nauk (2000) (in Russian) Novosibirsk. [8] A.A. Chesnokov, Exact solutions of the one-dimensional Russo–Smereka kinetic equation, J. Appl. Mech. Tech. Phys. 41 (2000) 593–603. [9] G. Russo, V.M. Teshukov, A.A. Chesnokov, Special class of solutions of the kinetic equation of a bubbly fluid, J. Appl. Mech. Tech. Phys. 46 (2005) 176–184. [10] A.K. Khe, A.A. Chesnokov, Propagation of nonlinear perturbations in a quasineutral collisionless plasma, J. Appl. Mech. Tech. Phys. 53 (2012) 657–663. [11] T. Stix, The Theory of Plasma Waves, McGarw-Hill, New York, 1962. [12] V. Teshukov, G. Russo, A. Chesnokov, Analytical and numerical solutions of the shallow water equations for 2-D rotational flows, Math. Models Methods Appl. Sci. 14 (2004) 1451–1479. [13] A.A. Chesnokov, A.K. Khe, Is Landau damping possible in a shear fluid flow? Stud. Appl. Math. 131 (2013) 343–358. [14] H. Nessyahu, E. Tadmor, Non-oscillatory central differencing schemes for hyperbolic conservation laws, J. Comput. Phys. 87 (1990) 408–463. [15] B.A. Kupershmidt, Deformations of integrable systems, Proc. R. Ir. Acad. A 83 (1983) 45–74. [16] M.V. Pavlov, Integrability of the Gibbons–Tsarev system, Amer. Math. Soc. Transl. Ser. 224 (2008) 247–253. [17] M.V. Pavlov, Algebro-geometric approach in the theory of integrable hydrodynamic type systems, Comm. Math. Phys. 272 (2007) 469–505. [18] D.J. Benney, Some properties of long nonlinear waves, Stud. Appl. Math. 52 (1973) 45–50. [19] J. Gibbons, S.P. Tsarev, Conformal maps and reductions of the Benney equations, Phys. Lett. A 258 (1999) 263–271. [20] M.V. Pavlov, Kupershmidt hydrodynamic chains and lattices, Int. Math. Res. Not. (2006). [21] G. Russo, P. Smereka, Kinetic theory for bubble flow II: fluid dynamic limit, SIAM J. Appl. Math. 56 (1996) 358–371. [22] R.C. Davidson, Methods in Nonlinear Plasma Theory, Academic Press, New York, 1972. [23] A. Harten, B. Engquist, S. Osher, S. Chakravarthy, Uniformly high order accurate essentially non-oscillatory schemes III, J. Comput. Phys. 71 (1987) 231–303. [24] A. Kurganov, E. Tadmor, New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations, J. Comput. Phys. 160 (2000) 214–282. [25] A.V. Gurevich, L.P. Pitaevskii, Nonlinear dynamics of a rarefied plasma and ionospheric aerodynamics, in: Reviews of Plasma Physics. Vol. 10, Atomizdat, Moscow, 1980, pp. 3–87 (in Russian).
Contents lists available at ScienceDirect
Physica D journal homepage: www.elsevier.com/locate/physd
The Russo–Smereka kinetic equation: Conservation laws, reductions and numerical solutions Alexander A. Chesnokov a,b,∗ , Maxim V. Pavlov c,d a
Lavrentyev Institute of Hydrodynamics, 15 Lavrentyev Ave., Novosibirsk 630090, Russia
b
Novosibirsk State University, 2 Pirogova Str., Novosibirsk 630090, Russia
c
Lebedev Physical Institute of Russian Academy of Sciences, Leninskij Prospekt 53, 119991 Moscow, Russia
d
National Research Nuclear University MEPHI, Kashirskoe Shosse 31, 115409 Moscow, Russia
highlights • • • •
The Russo–Smereka equation is associated with an integrable hydrodynamic chain. A stability of bubbly flow in terms of hyperbolicity of the model is studied. Infinitely many particular solutions for the kinetic equation are found. A finite-dimensional approximation for the model is proposed.
article
info
Article history: Received 2 February 2014 Received in revised form 30 March 2015 Accepted 31 March 2015 Available online 11 April 2015 Communicated by J. Bronski
abstract The one-dimension Russo–Smereka kinetic equation describing the propagation of nonlinear concentration waves in a rarefied bubbly fluid is considered. Stability of the bubbly flow in terms of hyperbolicity of the kinetic equation is studied. It is proved that a hydrodynamical chain associated with the Russo–Smereka kinetic equation possesses infinitely many conservation laws. Reductions of the model to finite component systems are derived. Conservative form of the kinetic model is proposed and numerical solution of the Cauchy problem with discontinuous initial data is obtained. © 2015 Elsevier B.V. All rights reserved.
Keywords: Kinetic equation Bubbly flow Hyperbolicity Conservation laws Hydrodynamic chain
1. Introduction A kinetic theory based on the statistical description of the interaction of a large number of bubbles has been developed for modelling of nonlinear waves in a rarefied bubbly flow [1–3]. The kinetic models taking into account the effect of collective interaction between bubbles are derived using the system of Hamilton’s ODEs describing the motion of individual bubbles. To obtain this system of equations one needs to know the kinetic energy of the fluid [4]. Assuming that all bubbles are rigid massless spheres
∗ Corresponding author at: Lavrentyev Institute of Hydrodynamics, 15 Lavrentyev Ave., Novosibirsk 630090, Russia. E-mail addresses: [email protected] (A.A. Chesnokov), [email protected] (M.V. Pavlov). http://dx.doi.org/10.1016/j.physd.2015.03.013 0167-2789/© 2015 Elsevier B.V. All rights reserved.
of the same radius Russo and Smereka [1] approximately calculated the energy and Hamiltonian of a bubble motion and obtained the kinetic equation for the evolution of the one-particle distribution function. This model is analogous to the Vlasov equation for plasma flow and describes the collective behaviour of a large number of bubbles subject to long-range interactions, modelled by selfconsistent field. The characteristic properties of the Russo–Smereka kinetic equation for the case of one space variable are studied by Teshukov [5] on the basis of a generalized theory of characteristics and hyperbolicity of integro-differential equations [6,7]. In [5] hyperbolicity conditions of the model are formulated and Riemann invariants and infinite series of conservation laws are found. The exact solutions of the Russo–Smereka kinetic equation in the classes of travelling and simple waves, as well as solutions with linearly dependent Riemann invariants, are obtained and studied in [8,9].
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
In the paper, we study the stability of solutions in terms of hyperbolicity conditions for the Russo–Smereka kinetic equation for bubbly flow considered as an infinite-dimensional quasilinear system. It is known that any solution of the linearized Vlasov equation is stable if it is defined by a distribution function with a single maximum. We show that it is not true for the Russo–Smereka model. We propose a conservative form of the kinetic model for consideration of discontinuous solutions. Differential conservation laws approximating the Russo–Smereka kinetic equation are derived. These laws are used to perform numerical calculations of wave propagation in a rarefied bubbly flow showing the possibility of the kinetic roll-over (the formation of two peaks of the distribution function which originally had a single peak). This effect usually leads to instability of the flow. We consider a hydrodynamical chain associated with the Russo–Smereka kinetic equation and show that this chain is integrable. Some new reductions of the kinetic equation to finite component systems are also obtained and their properties are studied. The paper is organized as follows. In Section 2 we present different formulations of the Russo–Smereka kinetic model. Each formulation is more suitable for some specific application. We also recall hyperbolicity conditions for this integro-differential model and give an example of verification of the hyperbolicity conditions. Conservative formulation of the kinetic model for a bubbly flow is proposed in Section 3, as well as its finite-dimensional approximation in the form of a system of differential conservation laws with a large number of unknown variables. In Section 4 we introduce a hydrodynamic chain associated with the Russo–Smereka kinetic equation. We show that this chain possesses infinitely many conservation laws. In Section 5 we apply the method of hydrodynamic reductions for constructing infinitely many particular solutions of the Russo–Smereka kinetic model. Then we obtain finitecomponent reductions using generalized and step-functions. Numerical results are presented in Section 6. We show that during the evolution the kinetic roll-over of the distribution function is possible. Finally we draw some conclusions. 2. Mathematical model
51
Let us show that Eqs. (4) are derived as a consequence of Russo–Smereka model (1). We use f˜ to denote the distribution function in semi-Lagrangian coordinates: f˜ (t , x, λ) = f (t , x, p(t , x, λ)). Then, the derivatives are represented as
∂f ∂f ∂p ∂ f˜ = + , ∂t ∂t ∂p ∂t
∂ f˜ ∂f ∂f ∂p = + . ∂x ∂x ∂p ∂x
These formulae and Eqs. (3) imply the obvious equalities
∂ f˜ ∂ f˜ ∂f ∂f ∂j ∂f + (p − j) = + (p − j) + p = 0, ∂t ∂x ∂t ∂x ∂x ∂p which lead to the second equation of system (4). In what follows, we assume that the inequality pλ > 0 is satisfied. It provides invertibility of the change of variables. 2.1. Hyperbolicity conditions for the kinetic equation Eqs. (4) belong to the class of systems with operator coefficients of the form Ut + A⟨Ux ⟩ = 0,
(5)
for which a generalization of the hyperbolicity notion was proposed in [6]. Here, U(t , x, λ) is an unknown vector function and A is a non-local operator on a set of functions of λ. Characteristics of system (5) are determined by equation x′ (t ) = k(t , x), where k is an eigenvalue of problem (F, (A − kI )⟨ϕ⟩) = 0. The eigenfunctional F is defined on a set of functions of variable λ, while the values of t and x are considered to be fixed, and is sought in a class of locally integrable or generalized functions (ϕ(λ) is a smooth test function). System of Eqs. (5) is hyperbolic if all eigenvalues k are real and the set of relations on the characteristics (F, Ut + kUx ) = 0 is equivalent to (5). As it follows from [5,7], hyperbolicity conditions for Russo– Smereka equation (1) on a solution f (t , x, p) are formulated in terms of characteristic function
χ (k + j) = 1 − n + (k + j)
2
In the one-dimensional case the Russo–Smereka kinetic equation in dimensionless variables is written as follows [1,5]:
∂f ∂f ∂j ∂f + (p − j) + p = 0, ∂t ∂x ∂x ∂p j(t , x) = pf (t , x, p) dp.
(1)
Here f (t , x, p) ≥ 0 is the distribution function for bubbles in phase space, x and p are the position and momentum of the bubble, t is the time, and j is the first moment of the distribution function. It is assumed that f decreases rapidly at infinity or has compact support in p. The model is suitable for the description of a rarefied bubbly flow in the case of small pressure variations. The condition for a bubbly flow to be rarefied is given by the inequality n(t , x) =
f dp < 1.
(2)
It was shown in [5], to study the properties of kinetic equation (1), it is appropriate to make use of the Eulerian–Lagrangian coordinates x and λ by substitution the variable p = p(t , x, λ), where the function p(t , x, λ) is a solution of the Cauchy problem pt + (p − j)px = pjx ,
p(0, x, λ) = p0 (x, λ),
λ ∈ [0, 1]. (3) As a result, for the new functions p(t , x, λ) and f (t , x, λ), we obtain
f dp
(p − k − j)2
,
(6)
or, more precisely, in terms of its limit values on the real axis from upper and lower complex half-planes
χ ± (p) = 1 − n(t , x) + p2 ± π ip2
∂ f (t , x, p′ ) dp′ ∂ p′ p′ − p
∂ f (t , x, p) . ∂p
(7)
The above formula is obtained from (6) by integration by parts and application of Sokhotski–Plemelj formulae. According to [5,7], kinetic equation (1) is hyperbolic on the rapidly decreasing solution f (t , x, p) if the following conditions are valid
1 arg χ ± (p) = 0.
(8)
The argument increment is calculated when p changes from −∞ to ∞ at fixed values of variables t and x. If supp f is bounded, then conditions (8) take the form
1 arg χ + (p)/χ − (p) = 0, here the argument increment is calculated for p ∈ supp f . Hyperbolicity conditions (8) guarantee that characteristic equation χ (k + j) = 0 has no complex roots, and these conditions are necessary for the flow stability.
the integro-differential system of equations [5,7] 2.2. Example of verification of the hyperbolicity conditions
pt + (p − j)px − pjx = 0, ft + (p − j)fx = 0,
j=
ppλ f dλ.
(4)
Conditions (8) allow one to verify whether Russo–Smereka equation (1) is hyperbolic for given solution f (t , x, p). Following
52
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
satisfied at the single interior point of p = pc at which the distribution function f (p) reaches a local maximum. Since the distributions f = f (p) with one maximum obey the inequality (pc − p)f ′ (p) ≥ 0, one can conclude that the sign of the quantity Z 1 (pc ) = 1 − n + p2c
Fig. 1. Contours C + in the complex plane (Z 1 , Z 2 ) (the solid curve corresponds to a = 1/2; the dashed curve to a = 1).
f ′ (p)dp p − pc
depends on the point pc . Bubbly flow is stable if the extreme point pc is sufficiently close to zero (this guarantees that the inequality Z 1 (pc ) > 0 is fulfilled). Back to the previous example (Figs. 1 √and 2), one can conclude that the flow is stable for |a| < ac = 1/ 2 (unstable for |a| > ac ), where a = ac is the root of the equation Z 1 (a) = 0. 3. Conservative form of the model The evolution of a smooth solution of the hyperbolic system of equations can involve a gradient catastrophe. Therefore, further description of the solutions is possible only for the class of discontinuous functions. It leads to the necessity to formulate the model in a form of conservation laws. To do this, we present (4) in the form Ht + (p − j)H
x
jt + A2 − 3j /2
Fig. 2. Distribution function f (p) with one maximum (notation is the same as in Fig. 1).
[7,10], in the plane (Z 1 , Z 2 ) we construct a closed contour C consisting of the contours C + and C − . The contour C + is given parametrically by the equations Z 1 = Re{χ + (p)},
Z 2 = Im{χ + (p)},
where the complex functions χ ± (p) are defined by (7). A contour C − , which is symmetric about the Z 1 axis to the contour C + , is given by the same equation with the function χ − (p). If the point of Z 1 = 0, Z 2 = 0 lies in the domain bounded by the contour C , then the characteristic equation χ (k) = 0 has complex roots (the function χ (k) is given by (6)). Otherwise, the kinetic equation for the corresponding solution is hyperbolic. Fig. 1 shows the contours C + (the direction of circulation about the contours is positive, i.e. counter-clockwise) corresponding to the distribution function with one maximum f (p) =
1 √ exp −(p − a)2 ,
2 π
p ∈ (−∞, ∞).
This function is shown in Fig. 2 for a = 1/2 (solid curve) and a = 1 (dashed curve). As it can be seen from the graphs, hyperbolicity conditions (8) are fulfilled if the extreme point p = a is closely spaced to zero (solid curve in Fig. 2). Increasing the parameter a leads to an increment in the argument of the functions χ ± (dashed line in Fig. 1). It corresponds to existence of the complex characteristic roots and, consequently, the flow is unstable. In the theory of plasma waves the following result is known [11]: any solution of the one-dimensional linearized Vlasov equation is stable if it is defined by a distribution function with a single maximum. We show that in the kinetic theory for a bubbly flow this is not true. It is mentioned above, hyperbolicity conditions (8) are violated if the point Z 1 = 0, Z 2 = 0 is in the domain bounded 1 by the contour C . Let Z 1 → Z±∞ > 0 for p → ±∞. Then, this is possible only if the conditions Z 1 (p∗ ) < 0 and Z 2 (p∗ ) = 0 are satisfied at some point p∗ ∈ (−∞, ∞). The equality Z 2 (p∗ ) = 0 is
j=
2
= 0,
pλt + (p − j)pλ
x
= 0,
=0 pH dλ, A2 = p2 H dλ, H = pλ f . x
(9)
The first condition in (9) is the local conservation of a number of bubbles, the second condition is equivalent to the conservation of the function f along the trajectory, and the last equation is the conservation law for hydrodynamics momentum j. Note that Eqs. (9) are similar to conservation laws of the shallow water equations for shear flows [12,13]. As a rule, in the case of infinitedimensional systems to substantiate the correctness of the choice of conservation laws is difficult. In what follows, we consider solutions with small amplitude jumps. From conservation laws (9) we obtain the Hugoniot conditions at the shock front x = x(t ) moving with the velocity V = x′ (t ):
(p − j − V )H = 0, (p − j − V )pλ = 0, A2 − 3j2 /2 − Vj = 0.
(10)
Here [ϕ(t , x, λ)] = ϕ(t , x(t ) + 0, λ) − ϕ(t , x(t ) − 0, λ) is a jump of the function ϕ at the shock front. Calculating the ratio of continuous quantities (p − j − V )H and (p − j − V )pλ , we obtain the following jump relation: [f ] = 0. Let us show that systems (4) and (9) are equivalent for smooth solutions. System (9) is a consequence of Eqs. (4) obtained by simple transformations. Indeed differentiating with respect to λ the first equation in (4) gives the second one in (9). The following equality leads to the first equation in (9): 0 = pλt + ((p − j)pλ )x f + ft + (p − j)fx pλ
= Ht + (p − j)H x . Finally we note that
(pH )t + ((p − j)pH )x − pHjx = 0. Integrating this equation with respect to λ, we get the third equation in (9). Now we show that (4) is a consequence of (9). From the first two equations (9) we have ft + (p − j)fx = 0. Integrating with respect to λ the second equation in (9) gives the equality pt + (p − j)px − pjx = G(t , x),
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
53
where G(t , x) is an arbitrary function. Let us multiply this equation by H and integrate with respect to λ. Taking into account the equality
(pt + (p − j)px )H = (pH )t + ((p − j)pH )x , which holds by virtue of (9), after integration we have jt + A2 − 3j2 /2
x
= nG.
The third equation in (9) leads to the identity G ≡ 0 which proves equivalence of (4) and (9) for smooth solutions. 3.1. Approximate model Fig. 3. Graph of the function χ( ¯ k) for M = 5.
To derive differential conservation laws that approximate integro-differential model (9) we divide the segment λ ∈ [0, 1] into M intervals 0 = λ0 < λ1 < · · · < λM −1 < λM = 1 and introduce the variables pi + pi−1 , pi = p(t , x, λi ) hi = pi − pi−1 , p¯ i = 2 λi 1 H dλ. f¯i = hi λi−1 Taking into account the equality H dλ = f dp and using a piecewise constant approximation of the distribution function f (t , x, p) = f¯i (t , x), λi λi−1
A2 =
p ∈ [pi−1 , pi ]
i =1
pi
p2 f dp =
pi−1
M ¯ 3 fi h i
12
i =1
2
+
× j+
hk +
k=1
M
(11)
f¯i hi
−1
i=1 M
1 2 i =1
f¯i h2i −
M
1 qi − k
f¯i = 0
(13)
f¯i hi ,
j=
M
p¯ i f¯i hi .
i =1
i=1
f¯M < f¯M −1 < · · · < f¯l ,
that correspond to an approximation of the distribution function with a single maximum. Then the equation χ¯ (k) = 0 has a unique root k = ki on each interval k ∈ (qi−1 , qi ), i = 1, . . . , l − 1, l + 1, . . . , M. Indeed, function χ¯ (k) is monotonic and varies from −∞ to ∞ on the intervals (qi−1 , qi ) for i = 1, . . . , l − 1, and from ∞ to −∞ for i = l + 1, . . . , M. Consequently, system (11) is hyperbolic if Eq. (13) has two roots in the interval (ql−1 , ql ). Since χ (−j) = 1 − n > 0 and χ¯ (k) tends to −∞ as k → ql−1 (or k → ql ), then the inequalities
(ql−1 < −j < ql )
are the sufficient conditions for hyperbolicity of system (11). The typical form of the function χ¯ (k) for M = 5 is shown in Fig. 3. This dependence is obtained for the piecewise constant approximation of the distribution function with bounded support: f (p) = (cos p + 1)/8, p ∈ (−π , π ). Note that integro-differential model (4) has continuous characteristic spectrum k = p − j, where p ∈ supp f . 4. Local conservation laws
f¯i hi
i
hk .
Introducing infinitely many moments
k=1
To solve differential conservation laws (11) numerically, one can apply Godunov type methods. In this case, due to a large number of equations in system (11), it is convenient to use central schemes [14], which do not require exact or approximate solution of the Riemann problem. System (11) can be written in the following form ut + A(u)ux = 0,
M
pl−1 < 0 < pl
The quantities pci (t , x) included in (11) are given by the formulae p¯ i = −
qi−1 − k
−
where χ¯ (k) is a discrete analogue of the characteristic function (6), qi = pi − j and the moments n and j are given by the formulae
f¯1 < f¯2 < · · · < f¯l ,
+ f¯i hi p¯ 2i .
∂ ∂ hi + (¯pi − j)hi = 0, ∂t ∂x ∂ ¯ ∂ fi hi + p¯ i − j f¯i hi = 0, ∂t ∂x M ¯ 3 3 f i hi ∂ ∂j 2 ¯ + + fi hi p¯ i − j2 = 0. ∂t ∂ x i=1 12 2 i
1
In the hyperbolic case Eq. (13) has M + 1 real roots ki (i = 1, . . . , M + 1) in the interval (p0 − j, pM − j). Let the values f¯i satisfy the inequalities
Further, we integrate the first two equations in (9) with respect to λ over the intervals (λi−1 , λi ) and use the previous formulae. As a result we obtain a system of conservation laws consisting of 2M + 1 differential equations for unknown functions hi (t , x), f¯i (t , x), and j(t , x):
hi
M i=1
i =1
pH dλ = f¯i hi p¯ i , M
χ¯ (k) = 1 − n + (k + j)2
n=
we have the equalities
f¯it + (¯pi − j)f¯ix = 0. This means that the matrix A(u) has the following eigenvalues: k∗i = p¯ i − j (i = 1, . . . , M ). To find the other eigenvalues of A(u), one can use the equation
(12)
where u = (h1 , . . . , hM , f¯1 , . . . , f¯M , j) and A(u) is the corresponding matrix. This form is used to compute the characteristic values. As follows from (11), the functions f¯i satisfy the equations
Ak =
pk f dp (k = 1, 2, . . .)
kinetic model (1) leads to the modified Benney hydrodynamic chain
∂ Ak ∂ Ak+1 ∂ Ak ∂ A1 + − A1 − (k + 1)Ak = 0, ∂t ∂x ∂x ∂x
(14)
which was first derived by Kupershmidt (see details in [15]) in framework of the KP theory. This hydrodynamic chain is integrable. This means that (14) possesses infinitely many conservation laws (18), commuting flows and local Hamiltonian structures.
54
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
Suppose that some particular solution of (14) is found, i.e. Ak (t , x) are determined. Then the distribution function f (t , x, p) can be reconstructed. Indeed in this case one can look for solution in the form f (t , x, p) = F (p, A1 , A2 , . . .). Substituting this ansatz directly into (1) yields infinitely many equations which can be integrated (see details in [16,17]). This dependence can be presented in the Laurent series form (p → ∞) F =
1 p
+
A1 p2
+
A2 p3
+
A3 p4
+ ···.
(15)
Substitution of expansion (15) into kinetic model (1) implies Eqs. (14) again. As we have already mentioned modified Benney chain (14) possesses infinitely many local conservation laws. A generating function of these conservation laws can be obtained from kinetic model (1). The first equation in (1) divided by fp takes the form ft fp
+ (p − j)
fx
+ pjx = 0.
fp
Then taking into account (according to the theorem about differentiation of an implicit function) that pt = −ft /fp and px = −fx /fp , one can obtain
pt +
p2 2
1
= 0,
− jp
(16)
x
− P1 − P2 F − P3 F 2 − · · · ,
(17)
F and all conservation law densities Pk are the polynomials with respect to the moments Am . Explicit expressions can be obtained by substituting expansion (15) into (17) and equating factors of every degree of p. For instance, P1 = A1 , P2 = A2 − A21 , P3 = A3 − 3A1 A2 + 2A31 . Then substituting (17) into (16) and equating to zero the expressions of the same degrees of F leads to the following conservation laws
∂ Qk ∂ Pk + = 0, ∂t ∂x
(k = 1, 2, . . .)
(18)
where the first three polynomials Qk are Q1 = A2 − 3A22 /2, Q2 = A3 − 3A1 A2 + 2A31 , Q3 = A4 − 4A1 A3 − A22 + 8A21 A2 − 4A41 . The existence of an infinite number of conservation laws is a rare property for hydrodynamics models. Benney system of the long-wave theory is a well-known example of the model having an infinite number of conservation laws [18]. Both models, Russo–Smereka and Benney, are generalized hyperbolic in the sense of [6] and can be written in terms of the Riemann invariants conserving along the characteristics. 5. Method of hydrodynamic reductions The method of hydrodynamic reductions was established in [19] for Vlasov (collisionless Boltzmann) equation. In the case of Russo–Smereka kinetic model (1) the same approach is applicable. The method of hydrodynamic reductions is based on a natural observation: all moments Ak (t , x) of hydrodynamic chain (14) can be parameterized by N arbitrary functions am (t , x) only, then hydrodynamic chain (14) becomes N component hydrodynamic type system, where N is an arbitrary natural number. Example. The moment decomposition (εm are arbitrary constants) Ak =
1
N
k + 1 m=1
akt +
εm (am )k+1
(19)
(ak )2 2
− ak A1 (a) = 0,
(20)
x
where a = (a1 , . . . , aN ), the function A1 (a) is given by (19) for k = 1. It should be noted that system (19), (20) is not a sole reduction of hydrodynamic chain (14). Indeed one can compare (20) with (16). Thus reduction (20) follows directly from (16) replacing p by ak (t , x). Let us assume that all moments Ak depend on a finite number of field variables am (t , x) satisfying hydrodynamic type system (20). All equations (infinitely many) of hydrodynamic chain (14) must be consistent with each other. Without loss of generality (see details in [16,17]) we use the first equation from hydrodynamic chain (14) only (we assume summation with respect to the repeated index):
∂ A1 (a) k ∂ A1 (a) k ∂ A2 (a) k a + a − 3A1 (a) a = 0. ∂ ak t ∂ ak x ∂ ak x Taking into account (20) we obtain the following relations
where a generating function of conservation law densities p(F , A1 , A2 , . . .) can be found from (15), i.e. p=
transforms hydrodynamic chain (14) into the so-called ‘‘waterbag’’ reduction
∂ A1 (a) ∂ A1 (a) ∂ A1 (a) − 2A1 (a) m k ∂a ∂a ∂ ak ∂ A1 (a) ∂ A2 (a) k − ak + ax = 0, ∂ ak ∂ ak
am
where each bracket must vanish independently for each index k. The compatibility conditions
∂ ∂ Al ∂ ∂ Al = k i i k ∂a ∂a ∂a ∂a for each pair of distinct indices i and k lead to the overdetermined system for A1 (a(t , x))
( ai − ak )
∂ 2 A1 ∂ A1 ∂ A1 ∂ A1 ∂ A1 + k Rˆ i − i Rˆ k = 0, ∂ ai ∂ ak ∂a ∂a ∂a ∂a
i ̸= k.
(21)
Here Rˆ = Σ am ∂/∂ am is the scaling operator. A direct computation shows that all other higher equations of hydrodynamic chain (14) allow to reconstruct dependences Ak (a) in quadratures (we omit here this technical verification) if the function A1 (a) is a solution of overdetermined system (21). Since i ̸= k all second derivatives are expressible via the first derivatives except the second derivatives with coincided indices. Thus system (21) possesses a general solution parameterized by N arbitrary functions. A similar system for description of N component reductions for the remarkable Benney hydrodynamic chain was derived by Gibbons and Tsarev in [19]. By this reason we call Eqs. (21) the Gibbons–Tsarev system. In a separate paper [20] we proved that chain (14) is nothing but the so-called modified Benney hydrodynamic chain. This means that both chains have the same amount of local conservation laws. Some multi-parametric solutions of Gibbons–Tsarev system (21) were also found in [20]. Instead of substitution (19), one can choose any other particular solution of overdetermined system (21). One can see, for instance, that a family of particular solutions can be extracted from (21) by the natural ansatz A1 (a) =
N
gm (am ),
m=1
where functions gm (am ) can be found by direct substitution into (21). Indeed in this case overdetermined system becomes
ˆ i′ (ai ) = gi′ (ai )Rg ˆ k′ (ak ), gk′ (ak )Rg
i ̸= k.
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
ˆ i′ (ai ) = ai gi′′ (ai ) one can obtain Taking into account that Rg ak gk′′
where km are the roots of characteristic equation (13). This representation allows one to construct exact solutions in the form
(a ) = α gk (a ), ′
k
55
k
ri = ri0 = const (i ̸= s),
where α is an arbitrary constant. This means that gk′ (ak ) = βk (ak )α , where βk are extra arbitrary (integration) constants. Below for simplicity we restrict our consideration on two main interesting and important cases determined by the Dirac deltafunctional and Heaviside step-functional ansatz. We also recall special class of solutions [9] characterized by a linear relationship between the Riemann invariants. This class of solutions will be used for testing the numerical results in the framework of approximate model (11).
where f = f (t , x, λ) and
5.1. Generalized functions and reductions
R(t , x, λ) =
Under certain assumptions kinetic model (1) can be reduced to the system of differential equations. The following representation of the solution
rs = rs (k),
k = ks + j.
These relations define the functions p0 (k), . . . , pM (k), where k(t , x) is a solution of the equation kt + kkx = 0. The existence of Riemann invariants for system (22) follows from the fact that integro-differential equations (4) can be presented in the characteristic form [5] ft + (p − j)fx = 0,
Rt + (p − j)Rx = 0,
n(t , x) − 1 p(t , x, λ)
f (t , x, ν)pν dν
+
p(t , x, ν) − p(t , x, λ)
(23)
are the Riemann invariants (here n is given by formula (2)). 5.2. Special class of solutions
M
f =
ni (t , x)δ(p − pi (t , x))
We also consider solutions of Eq. (1) in the class of functions that are piecewise continuous in the variable p with a bounded support:
i=1
leads to the ‘‘frozen’’ bubbly flow
f = f 1 (t , x, p) θ (p − pl (t , x)) − θ (p − pr (t , x)) .
∂ ni ∂ + (pi − j)ni = 0, ∂t ∂x M ∂ pi ∂ p2i + − pi j = 0, j = ni p i . ∂t ∂x 2 i=1
Here θ is a Heaviside step function, pl and pr are the boundaries of the interval in the variable p, beyond which the distribution function is identically zero, and f 1 (t , x, p) ≥ 0 is a non-negative function which is continuously differentiable on the set {(t , x, p) : t ≥ 0, x ∈ R, p ∈ [pl , pr ]}. Substitution of this ansatz into (1) yields
Here δ(p) is a Dirac delta function. This system of equations has imaginary characteristic roots [21]. The concept of waterbag is known in plasma physics [22]. It allows one to obtain particular solutions of kinetic equation (1) described by closed system of equations. We represent the kth moment in the form [16]
∂f ∂f ∂j ∂f + (p − j) + p = 0, ∂t ∂x ∂x ∂p pr j(t , x) = pf (t , x, p) dp,
Ak =
1
M
k + 1 i=0
M
εi pik+1 ,
εi = 0,
i =0
where εi are arbitrary constants whose sum is equal to zero. Substitution of the above introduced moments into chain (14) yields a closed system of M + 1 equations for unknown functions pi (t , x):
∂ pi ∂ pi ∂j + (pi − j) − pi = 0, ∂t ∂x ∂x
j=
M 1
2 i =0
εi p2i .
(22)
This system coincides with conservation laws (11) only for the class of solutions f¯i = fi = const (in this case (11) is the exact consequence of (1)). Note that due to the following representation of the solution: M
f =
M
fi θ (p − pi−1 ) − θ (p − pi ) ≡
i =1
i=0
εi θ (p − pi ),
i =0
M
i
εi = 0,
fi+1 =
εi
k=0
system (22) can be obtained directly from kinetic model (1) (here θ is a Heaviside step function). This system of equations is hyperbolic under some conditions derived above. Waterbag reduction (22) can be written in Riemann invariants
∂ rm ∂ rm + (km + j) = 0 (m = 1, . . . , M + 1), ∂t ∂x M p −k −j n−1 m i rm = + fi ln , km + j p − k − j i − 1 m i =1
pl
(24)
∂ pl ∂ pl ∂j + (pl − j) − pl = 0, ∂t ∂x ∂x ∂ pr ∂ pr ∂j + (pr − j) − pr = 0, ∂t ∂x ∂x where f = f 1 .
The class of solutions of kinetic model (1) with a bounded support characterized by a linear relationship between the Riemann integral invariants was obtained in [9] on the basis of the following property: if the functions f (t , x, p), pl (t , x), and pr (t , x) are a solution of system (24) then the function R(t , x, p), which has been defined above by (23) in semi-Lagrangian coordinates, satisfies the equation Rt + (p − j)Rx + pjx Rp = 0. Linear relationship between the invariants R and f R + f π cot(µπ ) − b = 0
(25)
(µ and b are constants) leads to a special class of solutions [9]. In this case the distribution function f has the form f =
sin(µπ )
1 − µb(pr − pl )
π
p
+b
p − pl
µ
pr − p
,
(26)
and the functions pl and pr satisfy system (24) with the following first moment
j(pl , pr ) = µ(pr − pl ) 1 +
1+µ 2
bpl +
1−µ 2
bpr .
(27)
Description of discontinuous solutions from this special class is given by the following closed system of conservation laws nt + (1 − n)j
x
= 0,
jt + A2 (n, j) − 3j2 /2
x
= 0.
(28)
56
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
By virtue of (26) the moments n, j, and A2 are expressed in terms of pl and pr by (27) and formulae
n=1−
pl
µ
pr
A2 = µ(pr − pl )
+
1 − µb(pr − pl ) ,
1 − µ2 3
1+µ 2
pr +
1−µ 2
6. Numerical results
pl
b(pr − pl ) + bpl pr . 2
Relations n = n(pl , pr ) and j = j(pl , pr ) can be inverted if the Jacobian determinant D(n, j)
µ
pl µ = −(pr − pl ) 1 + pl − µ(pr − pl ) b D(pl , pr ) pl pr pr × 1 + pr − µ(pr − pl ) b is non-zero. In the case b = 0 this condition reads pl ̸= pr (pl and pr det
2
are not equal to zero). Below these conservation laws will be used to perform numerical simulation of discontinuous bubbly flows from special class of solution as well as to test results obtained on the basis of general model (11). Remark. In general relation (25) is not fulfilled at the shock front. This means we cannot use conservative form (28) for the description of discontinuous flows. Only the general system of conservation laws (11) can be applied for this purpose. Nevertheless, modelling of discontinuous flows with small amplitude jumps δ = [n] is possible on the basis of simplified Eqs. (28). Following [7,12], one can show that [R] = O(δ 2 ) and, consequently, Hugoniot conditions
(p − j − V )H = 0, R + f π cot(µπ ) − b = 0, A2 − 3j2 /2 − Vj = 0, which are performed for a special class of solutions, coincide with relations (10) if the second-order terms in δ are neglected. Notice that the previous jump conditions are consequence of the following conservative form Ht + (p − j)H
= 0, jt + (A2 − 3j2 /2)x = 0, x (R + f π cot(µπ ) − b)H t + (R + f π cot(µπ ) − b)(p − j)H x = 0
n+1/2
uj
= unj − Λf′j /2,
(Λ = 1t /1x)
ujn++11/2 = (unj + unj+1 )/2 + (u′j − u′j+1 )/8
− Λ f(
n+1/2 uj+1
n+1/2 uj
) − f(
(30)
) .
This scheme approximates the system of conservation laws of the form ut + (f(u))x = 0 including models (11) and (28). Here 1x is the spatial grid spacing, while 1t is the time-step satisfying the Courant condition, and unj = u(t n , xj ). The computational domain on the x axis is divided into N cells, the cell centres are denoted by xj . Values u′j /1x and f′j /1x are approximations of the first-order derivatives with respect to x, calculated according to the ‘‘ENO limiter’’ procedure [23]. At t = 0 the initial data u0j are specified. The boundary conditions u(t n , x1−j ) = u(t n , x1 ) and u(t n , xN +j ) = u(t n , xN ) are used, which allow calculations to be performed until the initial perturbations reach the boundaries of the computation domain. Scheme (30) does not require exact or approximate solution of the Riemann problem that is very convenient in our case, due to a large number of equations in system (11). It should be noted that central schemes on non-staggered grids, such as Kurganov–Tadmor scheme [24] or, in general, schemes based on the local Lax–Friedrichs flux are also appropriate. Such non-staggered approach has the advantage that it is usually simpler to assign the boundary conditions. However, in the cases under study we consider trivial boundary conditions. Therefore, the Nessyahu–Tadmor scheme is perfectly appropriated.
Let us perform a comparison between the numerical solutions obtained for differential approximation (11) and for system (28) defining solutions from this special class. At t = 0 we define the function f by formula (26), where pl (x) = 0.5 and pr (x) = 2.25 if x ∈ (−0.5, 0), otherwise pl = 0.25 and pr = 2. Initial distribution function f (x, p) and moments n(x) and j(x) are shown in Figs. 4 and 5 by dashed lines. We calculate the solution in the domain x ∈ [−1, 1] for N = 150 and M = 120. To start calculations on the base of model (11) we specify unknown variables
5.3. Fluid dynamic limit Let the distribution function have the form
(p − p¯ )2 n f = √ exp − . 2T 2π T Calculation of moments of this distribution function A1 = j = np¯ ,
In this section we present some results related to the numerical modelling of bubbly flows. We implement here the Nessyahu– Tadmor second-order central scheme [14]
6.1. Comparison with a special class of solution
of system (4).
A0 = n,
Taking T = h21 /12, p¯ = p¯ 1 , and n = f¯1 h1 it is easy to see that system (29) coincides with the differential approximation (11) for M = 1. Note that the eigenvalues k1,2 are the roots of characteristic equation (13).
A2 = (¯p2 + T )n,
A3 = (¯p2 + 3T )np¯
(h1 , . . . , hM , f¯1 h1 , . . . , f¯M hM , j)
and their substitution into the first three conservation laws (18) leads to the closed system of equations for unknown functions n(t , x), p¯ (t , x) and T (t , x):
at t = 0. Here hi = h ≡ (pr − pl )/M; f¯i is the mean value of the function f (0, x, p) in the interval (pi−1 , pi ), where pi = pl + ih (this condition is imposed at the initial time). The following parameters for the special class of solutions are chosen µ = 0.6, b = 0. A sufficiently detailed resolution in the variable p is needed due to the rapid change of the distribution function in a special class of solutions in the neighbourhood of the points pl and pr . Figs. 4 and 5 (solid lines) show the results of computation at t = 1. This test confirms that the calculations of discontinuous solutions obtained by ‘‘multilayer’’ approximation (11) and by the system for special solutions (28) give close results, at least for small amplitude jumps.
nt + (¯un)x = 0,
(np¯ )t + (nT + np¯ u¯ − n2 p¯ 2 /2)x = 0, (nT + np¯ u¯ )t + (3T + (1 − 2n)¯p2 )nu¯ x = 0,
(29)
where u¯ = (1 − n)¯p. This system of equations, obtained in [21], is hyperbolic for T > np¯ 2 /3. In fact, let us rewrite (29) in the form (12), where u = (n, p¯ , T ). The eigenvalues of A(u) are k0 = (1 − n)¯p,
k1,2 = (1 − 2n)¯p ±
(1 − n)(3T − np¯ 2 ).
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
Fig. 4. Distribution function f at t = 0 (curves 1 and 2) and t = 1 (curves 3 and 4): 1—x ∈ (−0.5, 0); 2—x ̸∈ (−0.5, 0); curves 3 and 4 refer to calculation using models (28) and (11) at x = −0.2.
57
Fig. 6. Distribution function f at t = 0 (curve 1—x < 0, curve 2—x > 0) and at t = 2 (curve 3—x = −0.5, curve 4—x = 0.5).
Fig. 7. Moments n and j at t = 0 (curves 1, 2) and at t = 2 (curves 3 and 4). Fig. 5. Density n (curves 1–3) and moment j (curves 4–6): 1 and 4—t = 0; 2 and 5—t = 1 calculation using model (28); 3 and 6—t = 1 calculation using model (11).
6.2. Kinetic roll-over In the theory of quasineutral collisionless plasma flows the following result is known [10,25]: during the evolution the kinetic roll-over of the distribution function is possible (the formation of two peaks of the distribution function which originally had a single peak). Let us establish a similar property for the considered kinetic model for a bubbly flow. Suppose that at t = 0 the bubbly flow in the half-space x < 0 is defined by a distribution function f = fl (p) and in the halfspace x > 0 by a distribution function f = fr (p). In the interval p ∈ (pl , pr ) we choose the functions fl and fr as follows: fl (p) =
exp(−p2 )
, √ 5 π erf(1)
fr (p) =
exp(−4p2 )
√
5 π erf(2)
(pr = −pl = 1)
(curves 1 and 2 in Fig. 6); outside this interval fl = fr = 0. Distribution function at t = 2 for x = −0.5 and x = 0.5 is shown in Fig. 6 (curves 3 and 4, correspondingly). It can be seen that the distribution function for fluid more saturated by bubbles (curves 1 and 3) has changed in the minimum of p. At the same time the distribution function of the less saturated fluid (curves 2 and 4) has qualitative changes due to kinetic roll-over. In the vicinity of x = 0 (point of discontinuity at t = 0) the distribution function f for t > 0 has the following form: f → fl (p) for p < 0 and f → fr (p) for p > 0. Note that the formation of two local peaks of the distribution function usually leads to loss of the hyperbolicity of the kinetic model. Fig. 7 shows plots of the hydrodynamic moments at t = 2 (solid curves) and t = 0 (dotted). In the calculations we used the following resolution N = 400 and M = 100. Solution of the
Riemann problem with other initialdata (functions fl (p) and fr (p) have one maximum and pfl dp > pfr dp) has a similar form. Increase or decrease in the number of nodes has no significant effect on the numerical results. This remark is valid for all numerical examples presented in the section. 7. Conclusion Different formulations of Russo–Smereka kinetic equation (1) are presented. Semi-Lagrangian formulation allows one to rewrite this model as infinite-dimensional system of quasilinear equations (4) and propose its conservative form (9). This formulation gives possibility to consider discontinuous solutions. However a question about the physically correct conservative form of the model remains open. For this reason we restrict our consideration for discontinuous solutions with small amplitude jumps. The procedure of averaging with respect to the Lagrangian coordinate allows one to obtain finite-dimensional approximation in the form of differential conservation laws (11). Some known fluid dynamic limits (Eqs. (29) obtained under the assumption of local thermodynamic equilibrium and ‘‘waterbag’’ reduction (22)) are special cases of this system of equations. An example of verification of the hyperbolicity conditions of the kinetic model is given. It is established that the distribution function with one peak (Figs. 1 and 2) can lead to instability for a bubbly flow. This property distinguishes the Russo–Smereka model from the Vlasov equation because any solution of the onedimensional linearized Vlasov equation is stable if it is defined by a distribution function with a single maximum. It is established that Russo–Smereka kinetic model (1) is associated with an integrable hydrodynamic chain (14). It has infinitely many conservation laws
58
A.A. Chesnokov, M.V. Pavlov / Physica D 303 (2015) 50–58
(18). The method of hydrodynamical reductions is used to construct infinitely many particular solutions for the Russo–Smereka kinetic equation. Conservation laws (11) and (28) are used to perform numerical calculations of wave propagation in a bubbly flow initiated by discontinuous Cauchy data. It is shown the correspondence of the numerical results in the framework of the proposed finitedimensional approximation and in the special class of solutions characterized by a linear relationship between the Riemann integral invariants (Figs. 4 and 5). The effect of the kinetic roll-over of the distribution function is demonstrated (Fig. 6). As a rule it leads to instability of the flow. Acknowledgements The authors thank S.L. Gavrilyuk for stimulating discussions, and an anonymous referee for useful comments and suggestions. AAC’s work was supported by RFBR grant 13-01-00249 and Integration Project of SB RAS No. 30. MVP’s work was partially supported by RF Government grant 010-220-01-077, ag. G34.31.0005, by grant of Presidium of RAS ‘‘Fundamental Problems of Nonlinear Dynamics’’, and by RFBR grant 14-01-00389. References [1] G. Russo, P. Smereka, Kinetic theory for bubble flow I: collisionless case, SIAM J. Appl. Math. 56 (1996) 327–357. [2] H. Herrero, B. Lucquin-Desreux, B. Perthame, On the motion of dispersed bubbles in a potential flow, SIAM J. Appl. Math. 60 (1999) 61–83. [3] V.M. Teshukov, S.L. Gavrilyuk, Kinetic model for the motion of compressible bubbles in perfect fluid, Eur. J. Mech. B Fluids 21 (2002) 469–491. [4] L.M. Milne-Thomson, Theoretical Hydrodynamics, Macmillan, London, 1960. [5] V.M. Teshukov, Characteristics, conservation laws, and symmetries of the kinetic equations of motion of bubbles in a fluid, J. Appl. Mech. Tech. Phys. 40 (1999) 263–275.
[6] V.M. Teshukov, Hyperbolicity of long-wave equations, Dokl. Akad. Nauk 284 (1985) 555–559. [7] V.Yu. Liapidevskii, V.M. Teshukov, Mathematical models for long-wave propagation in an inhomogeneous fluid, Izd. Sib. Otd. Ross. Akad. Nauk (2000) (in Russian) Novosibirsk. [8] A.A. Chesnokov, Exact solutions of the one-dimensional Russo–Smereka kinetic equation, J. Appl. Mech. Tech. Phys. 41 (2000) 593–603. [9] G. Russo, V.M. Teshukov, A.A. Chesnokov, Special class of solutions of the kinetic equation of a bubbly fluid, J. Appl. Mech. Tech. Phys. 46 (2005) 176–184. [10] A.K. Khe, A.A. Chesnokov, Propagation of nonlinear perturbations in a quasineutral collisionless plasma, J. Appl. Mech. Tech. Phys. 53 (2012) 657–663. [11] T. Stix, The Theory of Plasma Waves, McGarw-Hill, New York, 1962. [12] V. Teshukov, G. Russo, A. Chesnokov, Analytical and numerical solutions of the shallow water equations for 2-D rotational flows, Math. Models Methods Appl. Sci. 14 (2004) 1451–1479. [13] A.A. Chesnokov, A.K. Khe, Is Landau damping possible in a shear fluid flow? Stud. Appl. Math. 131 (2013) 343–358. [14] H. Nessyahu, E. Tadmor, Non-oscillatory central differencing schemes for hyperbolic conservation laws, J. Comput. Phys. 87 (1990) 408–463. [15] B.A. Kupershmidt, Deformations of integrable systems, Proc. R. Ir. Acad. A 83 (1983) 45–74. [16] M.V. Pavlov, Integrability of the Gibbons–Tsarev system, Amer. Math. Soc. Transl. Ser. 224 (2008) 247–253. [17] M.V. Pavlov, Algebro-geometric approach in the theory of integrable hydrodynamic type systems, Comm. Math. Phys. 272 (2007) 469–505. [18] D.J. Benney, Some properties of long nonlinear waves, Stud. Appl. Math. 52 (1973) 45–50. [19] J. Gibbons, S.P. Tsarev, Conformal maps and reductions of the Benney equations, Phys. Lett. A 258 (1999) 263–271. [20] M.V. Pavlov, Kupershmidt hydrodynamic chains and lattices, Int. Math. Res. Not. (2006). [21] G. Russo, P. Smereka, Kinetic theory for bubble flow II: fluid dynamic limit, SIAM J. Appl. Math. 56 (1996) 358–371. [22] R.C. Davidson, Methods in Nonlinear Plasma Theory, Academic Press, New York, 1972. [23] A. Harten, B. Engquist, S. Osher, S. Chakravarthy, Uniformly high order accurate essentially non-oscillatory schemes III, J. Comput. Phys. 71 (1987) 231–303. [24] A. Kurganov, E. Tadmor, New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations, J. Comput. Phys. 160 (2000) 214–282. [25] A.V. Gurevich, L.P. Pitaevskii, Nonlinear dynamics of a rarefied plasma and ionospheric aerodynamics, in: Reviews of Plasma Physics. Vol. 10, Atomizdat, Moscow, 1980, pp. 3–87 (in Russian).