- Email: [email protected]
Efficient numerical approximation of compressible multi-material flow for unstructured meshes
Computers & Fluids 32 (2003) 571–605 www.elsevier.com/locate/compfluid
Efficient numerical approximation of compressible multi-material flow for unstructured meshes R. Abgrall a
a,b,* ,
B. Nkonga a, R. Saurel
c
Math ematiques Appliqu ees de Bordeaux, Universit e Bordeaux I, 351 Cours de la Liberation, 33405 Talence cedex, France b Institut Universitaire de France, France c IUSTI, Universit e Aix Marseille I, 13453 Marseille Cedex 13, France
Received 19 October 2000; received in revised form 31 January 2002; accepted 4 February 2002
Abstract We consider the numerical approximation of multi-dimensional multi-material flows. This is a difficult topic related to the numerical smearing of contact discontinuities (or material interfaces or slip lines). Any Eulerian scheme will produce an artificial mixing zone. In this artificial mixture, the computation of thermodynamical variables (pressure, sound speed, temperature, . . .) is difficult to achieve correctly. For the stiff cases considered in this paper, with solid–liquid–gas interfaces, small errors on the thermodynamical variables lead to the blow up of the computation in many cases. In this paper, we review and explain this problem, then we propose solutions for 1D and 2D flows. Contrarily to front-tracking techniques, our schemes use the same formulation everywhere on the mesh. We provide several examples with the stiffened gas equation of state (EOS): inert materials (water–air), and chemically reactive materials (solid explosive– Plexiglass–air). 2003 Elsevier Science Ltd. All rights reserved.
1. Introduction and physical model The computation of multi-material flows is more and more investigated because many physical situations involve interfaces and mixtures of several materials. For example, bubble-gas problems, interface problems, problems involving jets of different physical properties in rocket engines, etc. * Corresponding author. Address: Mathe´matiques Appliquees de Bordeaux, Universite Bordeaux I, 351 Cours de la Liberation, 33405 Talence cedex, France. Tel.: +33-5-56-84-60-68; fax: +33-5-56-84-26-26. E-mail addresses: [email protected] (R. Abgrall), [email protected] (B. Nkonga), richard@ iusti.univ-mrs.fr (R. Saurel).
0045-7930/03/$ - see front matter 2003 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 3 0 ( 0 2 ) 0 0 0 1 2 - 9
572
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
There exists a large literature about the numerical discretization of these problems: Eulerian schemes, Lagrangian schemes, front tracking techniques [5], level set methods [6], etc. We are interested in the development of Eulerian schemes where the equation of state (EOS) is described by parameters (such as the ratio of specific heats) that may vary within the computational domain, but the analytical formulation is the same everywhere. We also want to compute flows by using the same numerical formulation everywhere and by introducing few modifications of existing conservative schemes in order to make easy the evolution of existing codes. One of our goals is to construct schemes that can be modified easily to take into account creation of species by chemical processes, or able to consider viscous phenomena. This is not possible by the Eulerian schemes using the level set formalism, for example. In this paper, we consider gases modeled by the stiffened gas equation of state [11,12], but our technique could easily be adapted to more complex cases following for instance [25]. The stiffened gas EOS is a modification of the ideal gas EOS. It enables to consider complex situations such as compressible liquid–gas problems. More precisely, we consider problems where the pressure is related to the internal energy by: p ¼ ðcðx; tÞ 1Þq cðx; tÞp1 ðx; tÞ
ð1Þ
where the thermodynamical parameters c and p1 depend on time and space. This formulation permits to consider problems where the thermodynamical parameters are piece-wise constant. More complex situations, where the EOS is of the Mie Gr€ uneisen form, can also be considered by freezing in each computational cell the thermodynamical parameters. ‘‘Naive’’ extensions of the classical conservative schemes such as Roe scheme, often lead to severe numerical problems. First, the thermodynamical parameters may reach values out of their range of validity, and second, the numerical results appear to be oscillatory and non-physical at interfaces (negative pressure, negative arguments in the computation of the sound speed). This is not a drawback of a particular scheme, but a consequence of the conservative formulation of the schemes [1–4,15]. Many investigations have been done to define an efficient numerical approximation for multimaterial flows, see for example [1,3,14] for a review. Most of the time, these oscillations are barely visible, but in some cases (large pressure or density ratio at interfaces, very different EOS), they can lead to the blow up of the computation. Therefore, the major challenge is the numerical characterization of the material properties, in the mixing zone and the handling of pressure oscillations across material interfaces. The previous investigations done in this topic are always based on 1D analysis or in several dimensions but for very regular meshes [21]. In this basically 1D techniques, the system of PDEs is discretized in such a way that the mass, momentum and total energy are locally conserved, hence the correct shock speed can be computed. The flow composition is updated such that, in the case of a perfect contact discontinuity, the velocity and the pressure remain uniform. However, this method cannot be easily generalized to several dimensions, especially when the mesh is an unstructured or a very distorted structured mesh, because there are no natural direction in the mesh. In this paper, we show how to extend the method of [2,21] to general unstructured meshes. A careful analysis is conducted and enables to alleviate the remaining difficulties contained in [2] about high order extensions, and to proceed the multi-dimensional case easily. More precisely, we first discuss our model, then consider the continuous case to illustrate the general principles of the
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
573
paper in one and several dimensions. The next section is devoted to the illustration of these principles to 1D case (the Roe, Godunov and Rusanov schemes, with a particular emphasis on the Roe scheme); then, we discuss the multi-dimensional case, and show which additional modifications have to be made to take properly into account the physics of slip lines. We conclude by showing how the method can be used for flows with chemistry, and different test cases, inert and reactive, are discussed.
2. The physical model We consider a flow which dynamic is governed by the compressible Euler equations. In this context, the conservative formulation of the multi-fluid flow reads: ot q þ r ðquÞ ¼ 0 ot m þ r ðm uÞ þ rp ¼ 0
ð2Þ
ot E þ r ðqH uÞ ¼ 0 where q is the density, u the velocity, m ¼ qu the momentum, E the total energy, p the hydrodynamic pressure and H ¼ ðE þ pÞ=q the total enthalpy. This system is generally written in the following form: ot W þ divðFÞ ¼ 0 The pressure is related to the other variables by the equation of state (1). For simplicity reasons, we rewrite the ‘‘stiffened’’ gas EOS as: bp þ h ¼ E 12qu u where b :¼
1 c1
and h :¼
cp1 c1
are the thermodynamical parameters. Their evolutions are governed by the equations: ot b þ u rb ¼ 0;
and ot h þ u rh ¼ 0
ð3Þ
We note that an equivalent conservative form of the equations can be written, ot q þ r ðquÞ ¼ 0 ot ðquÞ þ r ðqu uÞ þ rp ¼ 0 ot E þ r ð qH uÞ ¼ 0
ð4Þ
ot ðqbÞ þ r ðqbuÞ ¼ 0 ot ðqhÞ þ r ðqhuÞ ¼ 0 Other equivalent systems could be obtained. We rely on this particular form to explain our technique, even though the other systems might also have been considered instead. The system (2), (3) or (4) is hyperbolic. In the 1D case, the eigenvalues are u, u c and u þ c where the sound speed is defined by:
574
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
c2 ¼
cðp þ p1 Þ q
Across a contact discontinuity, the pressure and the velocity remain constant. In the multidimensional case, the pressure and the normal velocity are continuous across a slip line. 3. General principles We begin with a description of the governing principles of our technique on system (2) in 1D. Then we recall some known problems in several dimensions; they already exist for perfect gas computations, but are more pronounced in the present one. We describe how to cure them. The application of these techniques, for 1D and 2D schemes, are given in Section 4.5. 3.1. The one-dimensional case Consider Dx > 0, xi ¼ iDx for i 2 Z, Dt > 0 and the conservative scheme dq þ mDFq ¼ 0 dqu þ mDFqu ¼ 0 dE þ mDFE ¼ 0
ð5Þ
In (5), we make use of the notation, for any function f, df : ¼ f nþ1 f n , Df : ¼ fiþ1=2 fi1=2 and m ¼ Dt=Dx. The question is to find a discretization of the non-conservative system (3) compatible with the physics of (2) and (3), namely the velocity and the pressure should stay continuous across material interfaces. It is well known that the smooth solutions of (2) are solutions of ot q þ uox q þ qox u ¼ 0
ð6Þ
qot u þ quox u þ ox p ¼ 0
ð7Þ
ot p þ uox p þ qc2 ox u ¼ 0
ð8Þ
The last equation is equivalent to the transport of physical entropy. The system (6)–(8) can be obtained from (2) and (3) by algebraic manipulations. The first equation is obtained by developing ot q þ ox qu. The second one is uðot q þ ox quÞ þ ot qu þ ox qu2 þ ox p ¼ 0 and the last one is obtained in two steps: • first u2 ot q þ uox q þ pox u ¼ ð ot E þ ox ½uðE þ pÞ Þ u ot qu þ ox qu2 þ ox p þ ðot q þ ox quÞ 2 • second, since q ¼ bp þ h, ot q þ uox q þ pox u ¼ pðot b þ uox bÞ þ ðot h þ uox hÞ þ b ot p þ uox p þ qc2 ox u
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
575
Our technique is to mimic these algebraic manipulations and to identify discrete operators that are consistent with the continuous operators between the brackets above. We have dE ¼ dðqÞ þ 12dðqu2 Þ. The idea is to write a scheme for the internal energy by combining the equations for E, qu and q. For an explicit scheme, dðqu2 Þ ¼ ðquÞ
nþ1
du þ un dqu
ð9Þ
dðquÞ ¼ qnþ1 du þ un dq and then, dðqÞ ¼ dE
un þ unþ1 un unþ1 dqu þ dq 2 2
ð10Þ
that is un þ unþ1 un unþ1 dðbp þ hÞ þ m DFE DFqu þ DFq ¼ 0 2 2
ð11Þ
We rewrite the flux variation of any Q-scheme [9], DF ¼ Fjþ1=2 Fj1=2 as DF ¼ Fjþ1=2 Fj1=2 ¼ 12 Wjþ1=2 þ Wj1=2 where Wjþ1=2 ¼ Fjþ1 Fj QðWjþ1 Wj Þ . In the sequel, we often call Wjþ1=2 a fluctuation. The fluctuation form of the scheme is DF ¼ 12 Wjþ1=2 þ Wj1=2 T
n nþ1 Denoting by ‘j the vector ‘j ¼ ððunj unþ1 j Þ=2; ðuj þ uj Þ=2; 1Þ , Eq. (11) can be rewritten as
dðbp þ hÞ ¼ dðbp þ hÞjþ1=2 þ dðbp þ hÞj1=2 where T dðbp þ hÞjþ1=2 ¼ 12‘j Wq ; Wqu ; WE jþ1=2 T dðbp þ hÞj1=2 ¼ 12‘j Wq ; Wqu ; WE j1=2 Then, we can rewrite, for an explicit scheme, dðbp þ hÞjþ1=2 as: dðbp þ hÞjþ1=2 ¼ pjn dbjþ1=2 þ dhjþ1=2 þ bnþ1 j dpjþ1=2 T The last step is to identify in 12‘j Wq ; Wqu ; WE jþ1=2 the terms that are consistent with n bnþ1 j dpjþ1=2 , pj dbjþ1=2 and dhjþ1=2 in the limit of a mesh refinement. A guide is that the evolution scheme for p should preserve the pressure in the case of a contact discontinuity, as in [2,22], so the idea is to write the fluctuations in term of variations of the primitive variables, in particular the velocity and the pressure. This can be done very naturally, because all the known schemes have, up to our knowledge, their natural formulation in term of these variables. Thus, it becomes easy to identify the sets of variations that are consistent with the density, velocity and pressure evolution equations. Since the model involves thermodynamical parameters, after having identified
576
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
the previous sets of variations, the remaining terms describe the discrete evolution of these thermodynamical parameters: this provides the fluctuations needed to update b and a. In the Sections 4.1, 4.2 and 4.3, we illustrate this technique for RoeÕs, GodunovÕs and RusanovÕs schemes and, in Section 4.4, for a second order MUSCL type scheme. We follow exactly these principle to make things clear. Remark 3.1. For an implicit scheme, we would have switched the superscript n þ 1 and n in (9) and made the modification accordingly in (10). 3.2. The multi-dimensional case In several space dimensions, the situation is a bit more complex, because slip lines are intrinsically ill represented in finite volume schemes, even though the defaults are generally weak. These problems appear first on the pressure, and then on the other variables, even in the perfect gas case, see [26]. Hence, we have to correct several problems: a multi-dimensional effect, and a problem due to the equation of state. In order to explain this, let us recall the partial differential equation (2) defining the behavior of in-viscid smooth flow, ot q þ r ðquÞ ¼ 0 ot m þ r ðm uÞ þ rp ¼ 0 ot E þ r ðqH uÞ ¼ 0 If C is any sub-cell of X, then for non-smooth solutions of Eq. (2) becomes: Z Z d W dx þ F n dr ¼ 0 dt C oC where F is the flux. This formulation suggests the finite volume scheme, where the building block is the solution, exact or approximate, of a Riemann problem, F ¼0 ot W þ on ~ ?
W ðg ¼ 0; g ; t ¼ 0Þ ¼
WL WR
if g < 0 if g > 0
ð12Þ
where ðg; g? Þ are the coordinates in the direct frame ðg; g? Þ. We have kgk ¼ 1 and g? is the unit vector normal to g at M 2 oC. We introduce U¼ug
and V ¼ u g?
This system reads: ot q þ og ðqUÞ ¼ 0 ot qu þ og ðquUÞ þ gx og p ¼ 0 ot qv þ og ðqvUÞ þ gy og p ¼ 0 ot E þ og ðqUHÞ ¼ 0
ð13Þ
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
577
and by simple manipulations, we get ot q þ og qU ¼ 0 ot qU þ og ðqU2 Þ þ og p ¼ 0
ð14Þ
ot qV þ og ðqUVÞ ¼ 0 ot E þ og ðqUH Þ ¼ 0
In practice, the numerical fluxes are constructed from (13) and not (14) because of conservation issues. However, this has some drawbacks. As pointed out in [26], the computations of slip lines can be inaccurate. If one consider for example a slip line initially aligned with the mesh, like on Fig. 1, where the normal velocity is non-zero, then even for a perfect gas, it is easy to see that for any up-wind solver, the pressure will not remain uniform. This phenomena is quite similar to what happens for mixture of perfect gases and is due to the non-correct computation of the tangential velocity. Consider for example a supersonic slip line, in the situation sketched in Fig. 1. With obvious notations (the slip line is initially located at i ¼ 0), and for a positive velocity u, we have n n n ¼ q mu q q ð15Þ qnþ1 1;j 1;j 1;j 0;j ðquÞnþ1 1;j
¼
ðquÞn1;j
mu
2
qn1;j
qn0;j
ð16Þ
nþ1 n n n ðqvÞ1;j ¼ ðqvÞ1;j mu ðqvÞ1;j ðqv Þ0;j
ð17Þ
nþ1 n n n ¼ E1;j mu E1;j E1;j E0;j
ð18Þ
From (15) and (16), we get unþ1 1;j ¼ u: the x-component of velocity remains uniform. From (15) and (17), the y-component of the velocity at ð1; jÞ lies in between of v0;j and v1;j under the CFL condition um 6 1, namely
Fig. 1. Slip line configuration.
578
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
vnþ1 1;j ¼
qn1;j ð1 muÞvn1;j þ qn0;j muv0;j qn1;j ð1 muÞ þ qn0;j mu
ð19Þ
From (18), using the perfect gas EOS and the uniformity of ui;j at times tn and tnþ1 , and the uniformity of the pressure at time tn , we get 2 nþ1 n p1;j ðvn1;j Þ p1;j ðvnþ1 u 1;j Þ þ ¼ þ m v21;j v20;j 2 c1 2 c1 2 2
In order to recover the uniformity of the pressure, we need 2 ðvn1;j Þ ðvnþ1 u 1;j Þ ¼ m v21;j v20;j 2 2 2 2
which is generally incompatible with (19). This pollutes the evaluation of the kinetic energy, and then that of the pressure. In the case of mixture of gases or the stiffened equation of state, this phenomena may lead to the blow up of the calculation, for example when vL vR is large, see [22] for example. The scheme proposed in [22] used a Cartesian mesh and the multi-dimensional scheme uses a splitting technique. A solution to this problem is to follow the same ideas as in 1D. We write the conservative scheme over one time step in fluctuation form. At one cell interface, the fluctuations depend locally on the variations of the density, the pressure, the normal velocity to the cell interface, and the variation of tangential velocity to the cell interface. In fact, the tangential velocity plays a special role. Like the normal velocity, it appears also through its square (in the kinetic energy). But contrarily to the normal velocity, the variation of the tangential kinetic energy does not disappear, by essence, through a slip line. Thus, we have two type of discrete equations involving the tangential velocity: one on the tangential velocity, and one on the tangential kinetic energy, like in the simplified example above. Thus there is an incompatibility. The idea is, like in [22], to introduce a special discrete equation on the tangential kinetic energy which is consistent with the continuous one, namely ð20Þ ot qV2 þ Uox qV2 ¼ 0 that can easily be obtained. At the beginning of every time step, we initialize the variables of this equation with the values obtained at the end of the previous time step. Then we compute the tangential kinetic energy with the discrete equation. We use this tangential kinetic energy to accurately compute the pressure (or the internal energy), and we proceed through the next time step. The discrete equation on the tangential kinetic energy is only an auxiliary equation. This method is the one developed in [22] for Cartesian meshes. In the case of a general mesh, the principles are the same, but the realization is a bit more subtle. This is developed in Section 4.5. 3.3. Additional comments We have presented here design principles only. There are many possible ways of modifying a scheme according to our technique, and this will become clear in the examples we develop in the next section.
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
579
However, our experience shows that the best modified schemes are those which satisfy additional constraints arising from the maximum principle. Eq. (3) are transport equations, thus b and h satisfy a maximum principle. This is also true for the tangential kinetic energy from (20). For these reasons, we look for schemes on b, h and the kinetic tangential energy that satisfy a discrete maximum principle.
4. Illustration on some schemes In this section, we illustrate the design principles of Section 3 to several numerical fluxes: the Roe scheme, the Godunov scheme the Rusanov scheme. We also explain how to apply this method to get high order accurate schemes in time. 4.1. Example of Roe’s scheme We first recall some basic facts about RoeÕs scheme [18]. The average Jacobian matrix JðWj ; Wjþ1 Þ satisfies the jump condition: F ðWiþ1 Þ F ðWi Þ ¼ JðWj ; Wjþ1 ÞðWiþ1 Wi Þ The Jacobian A of the flux F depends on the velocity u, the total enthalpy H and additional parameters, namely the derivative of the pressure with respect to the density and the internal energy, pq , p . To obtain an average Jacobian matrix that reduces to the standard one for a perfect gas, it is useful to note that it should depend on the Roe average of the velocity, the Roe average of the total enthalpy, and the following jump relation should be satisfied, f pm Djþ1=2 m þ f pE Djþ1=2 E þ pf Djþ1=2 p ¼ peq Djþ1=2 q þ f qb Djþ1=2 ðqbÞ þ p qh Djþ1=2 ðqhÞ
ð21Þ
There are many possible choice for peq , pe , etc. What is important is that at the end, the characteristic decomposition of Djþ1=2 W :¼ Wjþ1 Wj , in terms of the eigenvectors of JðWj ; Wjþ1 Þ, is very simple and only depend on averaged quantities such as the velocity, the enthalpy, the jþ1=2 , cjþ1=2 , etc. The calculations are given in speed of sound, etc., which are denoted by ujþ1=2 , H Appendix A, as well as all the components of the averaged matrix and the right and left eigenvectors. The numerical diffusion is controlled by the function u; for the Roe scheme without entropy correction uðuÞ ¼ juj. Setting u ðX Þ ¼ ðX uðX ÞÞ=2, Eq. (10) can be written explicitly as 1 0 D p Djþ1=2 p jþ1=2 A djþ1=2 q ¼ dbpjþ1=2 þ djþ1=2 h ¼ Ajþ1=2 u ðujþ1=2 Þ@Djþ1=2 q 2 Bjþ1=2 2 cjþ1=2 cjþ1=2 Cjþ1=2 qjþ1=2
Djþ1=2 u þ Djþ1=2 mu ðujþ1=2 Þpi Djþ1=2 b mu ðujþ1=2 ÞDjþ1=2 h cjþ1=2
ð22Þ
580
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
with Ajþ1=2 ¼ mð u2jþ1=2 2 ujþ1=2 e q jþ1=2 Þ u jþ1=2 þ 2e jþ1=2 þ ujþ1=2cjþ1=2 Þ e ujþ1=2 þ cjþ1=2 ÞÞ ðH u jþ1=2 ðujþ1=2 þ cjþ1=2 Þ þ e q jþ1=2 2Bjþ1=2 ¼ mu ðð jþ1=2 ujþ1=2cjþ1=2 Þ e þ mu ðð ujþ1=2 cjþ1=2 ÞÞ ðH u jþ1=2 ðujþ1=2 cjþ1=2 Þ þ e q jþ1=2 jþ1=2 þ ujþ1=2cjþ1=2 Þ e 2C ¼ mu ðð ujþ1=2 þ cjþ1=2 ÞÞ ðH u jþ1=2 ðujþ1=2 þ cjþ1=2 Þ þ e q jþ1=2 jþ1=2 ujþ1=2cjþ1=2 Þ e ujþ1=2 cjþ1=2 ÞÞ ðH u jþ1=2 ðujþ1=2 cjþ1=2 Þ þ e q jþ1=2 mu ðð Djþ1=2 ¼ mu ð ujþ1=2 Þðpjþ1=2 pi ÞDjþ1=2 b where e u jþ1=2 ¼
unj þ unþ1 j ; 2
pjþ1=2 ¼
pffiffiffiffiffi pffiffiffiffiffiffiffiffiffi qj pjþ1 þ qjþ1 pj pffiffiffiffiffi pffiffiffiffiffiffiffiffiffi ; qj þ qjþ1
and
e q jþ1=2 ¼
unj unþ1 j 2
qjþ1=2 ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffi qj qjþ1
The terms djþ1=2 q þ dj1=2 q with Djþ1=2 p Djþ1=2 p Djþ1=2 u þ þ Bjþ1=2 þ Cjþ1=2 q dj1=2 q ¼ Ajþ1=2 u Djþ1=2 q þ Djþ1=2 c2 c2 c are consistent with upx þ qc2 ux , i.e. in the limit of a mesh refinement, and for smooth solution, we get upx þ qc2 ux . n We also have dðbpÞjþ1=2 ¼ bnþ1 j dp þ pj djþ1=2 b, so that a natural discretization is dbjþ1=2 ¼ mu ð ujþ1=2 ÞDjþ1=2 b dhjþ1=2 ¼ mu ð ujþ1=2 ÞDjþ1=2 h dbj1=2 ¼ muþ ð uj1=2 ÞDj1=2 b dhj1=2 ¼ muþ ð uj1=2 ÞDj1=2 h
ð23Þ
and dbj ¼ dbjþ1=2 þ dbj1=2 dhj ¼ dhjþ1=2 þ dhj1=2 It clearly prevents oscillations in the pressure (the contacts are perfectly described). The scheme (23) also admits a maximum principle on b and h under the CFL-like condition: m max u ð ujþ1=2 Þ 6 1 j
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
581
Fig. 2. Solution of the Riemann problem between states Wj and Wjþ1=2 . Wjþ1=2 is the state on the line x ¼ xjþ1=2 .
Moreover, the conservation of the density, momentum and total energy is guaranteed by construction. 4.2. Example of Godunov scheme Consider FðWj ; Wjþ1 Þ ¼ F ðWjþ1=2 Þ the flux for Godunov scheme [10] for (4) between the left and right states Wj and Wjþ1 (see Fig. 2). As we have seen above, the important quantities are not the flux but the residuals Wjþ1=2 ¼ FðWjþ1 ; Wj Þ F ðWj Þ :¼ F ðWjþ1=2 Þ F ðWj Þ; Wj1=2 ¼ F ðWj Þ FðWj ; Wj1 Þ :¼ F ðWj Þ F ðWj1=2 Þ
Thanks to the particular structure of GodunovÕs scheme, we can introduce Jjþ1=2 (resp. Jj1=2 ) and Wj (resp. Wj1=2 and Wj ). Calling the averaged entries the Roe average between the states Wjþ1=2 f of Jjþ1=2 (resp. Jj1=2 ) by f jþ1=2 and f jþ1=2 (resp. f j1=2 and f j1=2 , and setting Djþ1=2 ð:Þ ¼ ð:Þjþ1=2 ð:Þj (resp. Dj1=2 ð:Þ ¼ ð:Þj ð:Þj1=2 ), we can use the formulas (23) with the correct choice of Roe average to define the scheme, namely djþ1=2 b ¼ mujþ1=2 Djþ1=2 b djþ1=2 h ¼ mujþ1=2 Djþ1=2 h dj1=2 b ¼ muj1=2 Dj1=2 b
ð24Þ
dj1=2 h ¼ muj1=2 Dj1=2 h This scheme admits a maximum principle on b and h under a CFL condition. Mass, momentum and energy are also conserved.
582
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
4.3. Example of the Rusanov scheme The RusanovÕs flux [19] reads: F ðWi ; Wjþ1 Þ ¼ 12ðF ðWj Þ þ F ðWjþ1 Þ lðWjþ1 Wj Þ
so that the residuals can be written as: wjþ1=2 ¼ 12 Jjþ1=2 ljþ1=2 Id ðWjþ1 Wj Þ wj1=2 ¼ 12 Jj1=2 þ lj1=2 Id ðWj Wj1 Þ so that once more the formula (23) can be used with the functions u ð:Þ and uþ ð:Þ redefined by: u ðujþ1=2 Þ ¼ ðujþ1=2 ljþ1=2 Þ=2 and uþ ðuj1=2 Þ ¼ ðuj1=2 þ lj1=2 Þ=2, namely ðujþ1=2 ljþ1=2 Þ Djþ1=2 b 2 ðujþ1=2 ljþ1=2 Þ Djþ1=2 h djþ1=2 h ¼ m 2 ðuj1=2 þ lj1=2 Þ dj1=2 b ¼ m Dj1=2 b 2 ðuj1=2 þ lj1=2 Þ Dj1=2 h dj1=2 h ¼ m 2 djþ1=2 b ¼ m
ð25Þ
This scheme admits a maximum principle on b and h under a CFL condition. Mass, momentum and energy are also conserved. 4.4. Space and time accuracy We consider a predictor corrector scheme, and assume for simplicity that the predictor and corrector steps are or the form þ þ ; Wjþ1=2 Þ F ðWj1=2 ; Wj1=2 Þ ð26Þ Wjnþa ¼ Wjn am F ðWjþ1=2 where a ¼ 12 or 1, Wj1=2 are evaluated via the MUSCL method, but the same type of technique can also be extended to the Harten–YeeÕs scheme [27]. The trick is to rewrite the scheme in fluctuation form, Wjnþa ¼ Wjn am wajþ1=2 þ waj1=2
where þ ; Wjþ1=2 Þ F ðWj Þ :¼ wjþ1=2 þ wjþ1=4 wajþ1=2 ¼ F ðWjþ1=2 þ waj1=2 ¼ F ðWj Þ F ðWj1=2 ; Wj1=2 Þ :¼ wj1=2 þ wj1=4 ; þ ; Wjþ1=2 Þ F ðWjþ1=2 Þ, wjþ1=4 ¼ F ðWjþ1=2 Þ F ðWj Þ, so that we can apply the where wjþ1=2 ¼ F ðWjþ1=2 same machinery as before. If the base scheme is RoeÕs, GodunovÕs or RusanovÕs, we see that the
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
583
predictor and corrector are the one given by the formula (23), (24) or (25) evaluated at the corresponding states, and modified by the following corrective term djþ1=2 b ¼ m u ðujþ1=2 ÞDjþ1=2 b þ u ðujþ1=4 ÞDjþ1=4 b djþ1=2 h ¼ m u ðujþ1=2 ÞDjþ1=2 h þ u ðujþ1=4 ÞDjþ1=4 h ð27Þ dj1=2 b ¼ m uþ ðuj1=2 ÞDj1=2 b þ uþ ðujþ1=4 ÞDj1=4 b dj1=2 h ¼ m uþ ðuj1=2 ÞDj1=2 h þ uþ ðuj1=4 ÞDj1=4 h where the functions u ð:Þ and uþ ð:Þ depends of the solver used in the predictor (resp. corrector) step. We have used the notations: þ
Djþ1=2 ð:Þ ¼ ð:Þjþ1=2 ð:Þjþ1=2 ; Dj1=2 ð:Þ ¼ ð:Þþ j1=2 ð:Þj1=2 ;
Djþ1=4 ð:Þ ¼ ð:Þ jþ1=2 ð:Þj ; Dj1=4 ð:Þ ¼ ð:Þj ð:Þþ j1=2 ð:Þj1=2 , ð:Þj1=2 the average values between ð:Þþ j1=2 and ð:Þj1=2 , ð:Þjþ1=4 and ð:Þjþ1=4 the average þ values between ð:Þ jþ1=2 and ð:Þj , ð:Þj1=4 and ð:Þj1=4 the average values between ð:Þj and ð:Þj1=2 .
4.5. An illustration for multi-dimensional finite volume approximation In this section, we develop an extension of the previous scheme to multi-dimensional problems. We start from a standard conservative numerical scheme for (2). The scheme is Dt X / ð28Þ Wi nþ1 ¼ Wi n ai j2#ðiÞ ij with the conserved variables are W ¼ ðq; m; EÞT , Dt is the time step, ai is the area of the cell (see Fig. 3), #ðiÞ is the set of the neighboring cells of the cell i and /ij is the fluctuation associated to the cell interface ðijÞ. This fluctuation is computed from the solution of a Riemann problem at the cell interface. We illustrate the approach in the case of the Roe scheme. We use the notations of Appendix B where construction of Roe scheme in 2D is sketched. As pointed out above, the method can easily be extended to Godunov and RusanovÕ schemes, following the same ideas as in Sections 4.1, 4.2 and 4.3. We start with some notations: T uu V2 ; Q? ¼ ; W ¼ ðq; qU; EÞT ; F ¼ qU; qU2 þ p; qH U ; 2 2 The system (4), in the local coordinates is given by: 8 ot W þ og F ¼ 0; > > < o ðqVÞ þ o ðqUVÞ ¼ 0; Q¼
t
g
ot ðqQ? Þ þ og ðqUQ? Þ ¼ 0; > > 2 : bp þ h ¼ E qU2 qQ?
Z ¼ ðW; qVÞT
ð29Þ
584
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Fig. 3. Dual cell Ci associated with triangulation vertex Mi in 2D space.
We also set 0
W
1
C B B qV C C Z¼B B qb C; A @
0
1
C B B qVU C C G¼B B qbU C; A @ qhU
qh W ¼ M1 ij Z
F
0
1
0
B0 g B x Mij ¼ B B0 0 @ 0 g? x
0 gy 0 g? y
0
1
0C C C; 1C A 0
Z ¼ Mij W ;
and J ¼ oZ G
The numerical dissipation of Roe scheme is
e ÞDZ ¼ uðUÞðDðqÞ nÞr1 þ uðU c~Þ n H r2 þ uðU þ c~Þ n þ H r3 uð J 2 2 Db bDh hDb þ uðUÞ qDðVÞ r4 þ uðUÞq r5 þ uðUÞ r6 b qb
ð30Þ
with, see Appendix B, Dp DU DU and H ¼ q e c~2 c The fluctuation Wij at the cell interface is given by: 0 Z1 Wij n o 1 B qb C ~ ~ Wij ¼ DG uðJÞDZ ¼ u ðJÞDZ ¼ @ Wij A 2 Wqh ij n¼
ð31Þ
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
585
Fig. 4. Area (ai ) of the cell and partial area (aij ) associated to a given segment ði; jÞ.
Therefore, the update of the conservatives variables W on the internal cells, is Wi nþ1 ¼ Wi n
X aij Dt X 1 Z nþ1 Mij Wij ¼ M1 ij Zij ai j2#ðiÞ a j2#ðiÞ i
where Znþ1 ¼ Znij ðDt=aij ÞWZ ij . Here, aij is the component of the cell area associated to the ij interface (see Fig. 4). We denote by Wij the first three components of WZ ij and obtain: Wnþ1 ¼ Wnij ij
Dt Wij aij
ð32Þ
This relation provides the updated variables associated to the sub-cell aij of the cell ai . Therefore, the projection on the cell ai of the variables is obtained by averaging the sub-cells variables weighted by the ratio aij =ai . This is one of the keys of the approach proposed in this paper. Indeed, with this decomposition, we can focus the study on the sub-cells aij : the situation is now very similar to the 1D case. Let us denote dXij ¼ Xijnþ1 Xijn and K ¼ qQ ¼ 12qU2 þ 12qV2 ¼ K þ K ? Here, the local frame is defined by the normal to aij at the interface between the cells ai and aj , and the orthogonal vector, such that the frame is direct. We have 1 nþ1 nþ1 qij Uij dUij þ Unij dðqUÞij dEij ¼ dðqÞij þ dKij? þ dKij ; dKij ¼ 2 n nþ1 and dðqUÞij ¼ qij dUij þ Uij dqij . For simplicity, we denote: Uni Unþ1 Unþ1 þ Uni ij ij T ¼ ; Uij ¼ and ‘ij ¼ ðQk ij ; Uij ; 1Þ 2 2 k ? Therefore, dðqÞij þ dKij ¼ Qij dqij Uij dðqUÞij þ dEij so that denoting rij ¼ Dt=aij Qk ij
n ? bnþ1 ij dpij þ pi dbij þ dhij þ dKij ¼ rij ‘ij Wij
ð33Þ
586
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Setting u ðX Þ ¼ ðX uðX ÞÞ=2, the components of the numerical fluctuation are given by: Wqij ¼ u ðUÞðDq nÞ þ c~ a 1 n þ a2 H p q WqU ij ¼ UWij þ Wij
WEij ¼ HWqij þ N where Wpij ¼ c~2 a 2 n þ a1 H , c2 ðDq nÞ þ qVDV þ pDb þ Dh ; N ¼ UWpij þ u ðUÞ b~ a 1 ¼
u ðU þ c~Þ þ u ðU c~Þ ; 2~ c
a 2 ¼
u ðU þ c~Þ u ðU c~Þ 2~ c
Using these expressions of the fluctuations associated to the edge ði; jÞ, we compute: 2 ~ ? Wq þ u ðUÞqVDV ‘ij Wij ¼ Sij þ b c~a Dp þ a q~ c DU þ u ðUÞðpi Db þ DhÞ þ Q 1 2 ij
ð34Þ
where k ~ k Wq þ Wpij U U þ u ðUÞðp pi Þ Sij ¼ Qij Uij U þ Q ij ij In general, Sij is consistent with zero and equal to zero if the velocity and the pressure field is a ~ c2 DU is consistent with space–time constant. rij b~c~a 1 Dp is consistent with bUog p and rij ba2 q~ 2 bqc og U. Therefore a natural discretization of the EOS parameters is obtained by a combination of Eqs. (33) and (34) and an identification of the different terms with respect to the non-conservative formulation in the direction gij : dbij ¼ rij u ðUÞDij b
and dhij ¼ rij u ðUÞDij h
ð35Þ
If we denote dbji and dhji the variations on the cell j associated to the interface, it is easy to see that: dbji ¼ rji uþ ðUÞDij b and dhij ¼ rji uþ ðUÞDij h where rji ¼ Dt=aji . These relations are obtained under the assumption that the update of the tangential kinetic e ? Wq þ u ðUÞqVDV : We set energy follows: dKij? ’ rij Q ij ? e ? Wq þ u ðUÞqVDV WKij ¼ Q ij
ð36Þ
Because the stability of the contact discontinuity is marginal in 2D and 3D [24], the previous is generally violated when the tangential velocity is not continuous, even for perfect gas EOS. This is why we independently discretize the tangential kinetic energy, namely ðKij? Þnþ1 ¼ ðKij? Þn rij WKij
?
In order to take into account this auxiliary relation in the numerical approach, we define a local update of the internal energy q by: nþ1 ðqÞij
¼
Eijnþ1
nþ1 2 ðqnþ1 ij Uij Þ
2qnþ1 ij
ðKij? Þ
nþ1
¼ Ein Kin rij WEij þ rij WKij
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
587
where qU ðUni Þ2 q rij Wij qU n q ¼ þ þ W W Wij þ 2U i ij ij 2 2qni From this discretization, we can construct an average internal energy on the cell that is locally consistent. Of course, this energy is different to the one that can be computed from the conservative variables on the cell at the time tnþ1 . Therefore we propose to update the total energy by the b nþ1 satisfying the relation: value E i nþ1 nþ1 X aij b nþ1 mi mi ðqÞnþ1 ¼ E ij i a 2qnþ1 i j2#ðiÞ i
WKij
? WKij
Uni WqU ij
This means that the modified value of the energy at the time tnþ1 is given by a non-conservative formulation: mnþ1 mnþ1 mn mn X E K i i b nþ1 ¼ En Dt E W i n i ð37Þ W ij ij þ i i ai j2#ðiÞ 2qi 2qnþ1 i The update of the material parameters in the cell is defined as follows: Dt X Dt X n nþ1 n bnþ1 ¼ b u ðUÞD b and h ¼ h u ðUÞDij h ij i i i i ai j2#ðiÞ ai j2#ðiÞ
ð38Þ
Summary: The numerical scheme is then defined by the classical conservative formulation for the mass and momentum. The total energy update is obtained by a non-conservative formulation (37). The material parameters are given by (38). The resulting non-conservative scheme is a more general formulation of the idea developed in [21] with a splitting approach. A similar procedure is used in [3]: the important variables are the primitive variables and not the conservative ones. To illustrate this point, we present a quasi-1D test case where this default affects dramatically the solution of the Euler equation in 2D (see Fig. 8). 5. The reactive model We now detail the extension of this method for materials that the physical state changes during chemical reaction: solid energetic materials. This can easily be done because the formulation is the same everywhere. Solid explosives start from the solid state before reaction, then become a solid–gas mixture during reaction and a pure gas mixture at the end of the reaction. We first give some basis regarding the thermodynamics of such materials. We then give the some elements regarding their chemistry. We finally detail the way to account of the equation of state changes during the decomposition process and the way to couple it with the numerical treatment of the interfaces. This is important for all applications involving interactions of detonations with neighboring materials. 5.1. Flow model When dealing with detonation waves in condensed explosives the pressure and velocity conditions are so high that liquids and solids become compressible and governed by the gas dynamics
588
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
equations. Thus, system (2) becomes valid even for solids. The closure of such system is achieved by Eq. (1) that needs some modifications to account for changes in chemical composition of the mixture. This will be detailed in the next subsection. However, let us mention that the system (2) is not sufficient when solving detonation waves with an Eulerian method as detailed previously. It must be complemented by an equation allowing the determination of the phases concentration, varying inside the reaction zone of the detonation wave. We define the gas volume fraction k, whose initial value is zero in the solid explosive and that tends to the unity at the end of the reaction. The explosive being moving, the volume fraction equation reads: ð39Þ ot qk þ r ðqukÞ ¼ qk_f where k_ represents the rate of creation of the gas phase. Assuming that the rate of decomposition obeys a simple pressure kinetics (Vieille law), we get: k_ ¼ ð1 kÞapn
ð40Þ
where a and n are the kinetics parameters of the explosive. The function f is an indicator that allows the determination of the explosive concentration in the mesh. Indeed, the decomposition law must be applied only at points where the explosive is present and not in inert materials. Spatial localization of the explosive is realized by the introduction of function f that equals zero initially everywhere, except in the explosive where it equals 1. Its evolution is given by: ot qf þ r ðquf Þ ¼ 0
ð41Þ
5.2. Mixture equation of state During the reaction process the explosive passes from solid state to gaseous products. The equation of state (1) must be reformulated to deal with this transformation. During the detonation propagation, there are three states in the flow: state 1 where the solid explosive is intact, state 2 where the explosive is completely reacted and composed exclusively of gases and an intermediate state composed of a mixture of solid and gases inside the reaction zone. To deal with this various states we propose a mixture equation of state able to reproduce all of them in a unique formulation. Each pure fluid is assumed to obey the stiffened gas EOS. Thus the internal energy of the solid reads: s ¼
bs þ hs þ cs qs
with
cs ¼ s0
hs qs0
where s0 , cs and qs0 are respectively the formation energy of the solid, the internal energy and density under standard conditions. For the gas phase, we assume the same formulation but instead of choosing standard conditions as reference, we retain the Chapman-Jouguet (CJ) point. Thus the internal energy reads: g ¼
bg þ hg þ cg qg
with cg ¼ cj
hcj qcj
where hcj ¼ ðpcj þ cg p1;g Þ=ðcg 1Þ and the characteristics states qcj , cj and pcj can be obtained from a thermochemical code like CHEETAH [8].
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
589
With the internal energies definition, the mixture equation of state is readily obtained under the assumption of pressure equilibrium between phases. This assumption is justified because of the small size of the solid particles into the mixture and of the high speed of sound of both two phases. Indeed, the pressure relaxation time is closed to the time required for an acoustic wave to cross a solid grain diameter. The mixture internal energy reads: q ¼ kqg g þ ð1 kÞqs s
ð42Þ
and now is a function of q, p and k. According to the previous definitions for the internal energies (s and g ) and assuming that gas and solid are at the same pressure p ¼ pg ¼ ps , we recover an EOS of the same type as the pure phases. It reads: ¼
bþh þc q
with
cðkÞ ¼ kcg þ ð1 kÞcs ; bðkÞ ¼ kbg þ ð1 kÞbs ; hðkÞ ¼ khg þ ð1 kÞhs
ð43Þ
It is now necessary to determine the decomposition parameters of the kinetics law of the explosive (parameters of the Vieille law). 5.3. Chemical kinetic The parameters of the VieilleÕs law are not the standard ones as determined in stand burner or closed bomb. Indeed, the parameters determined by such a way are those of the deflagration regime of the explosive. At the pressure levels we are interested in, these parameters cannot be determined by direct measurement. Also, the kinetic parameters we are looking for are not related to the ignition and transition regime of the explosive. Such mechanisms cannot be summarized in a simple Vieille law [16]. The kinetic parameters needed here are related only to the self-detonation propagation regime of the explosive. To determine these parameters we use an experimental observation based on the front curvature of the detonation wave. The energy responsible for the front propagation is liberated between the shock front and the sonic plane (in the frame of reference related to the shock front). The distance between these two planes being finite, multi-dimensional rare-fraction waves coming from the boundaries of the explosive product are able to propagate inside the reaction zone of the explosive. These effects result in a curvature of the shock front. Thus, it is clear that the shock curvature is a result of the coupling effects between the chemical kinetic of the explosive and the various hydrodynamics effects, especially the rare-fraction waves. Considering a cylindrical charge of given diameter, the corresponding detonation wave propagates at a given velocity with a given curvature. By measuring these variables for various charge diameters we obtain a relation between the front curvature and its velocity. This relation is intrinsic to the explosive. It is a result of its chemical kinetic and rare-fraction waves. This relation will enable us to determine the kinetic parameters in a wide range of pressures. Indeed, when the diameter varies, the detonation velocity and the pressure inside the reaction zone vary strongly. By using the equation of state for the mixture, combined with system (2) and Eq. (39) in the context of a steady 2D axi-symmetric flow in the frame of the shock and specialized onto the symmetry axis, we get an ODE system easy to solve. It corresponds to the Wood and Kirkwood
590
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
model [7]. Giving an arbitrary guess for the chemical kinetic parameter of the VieilleÕs law and the front curvature for a given charge diameter this model furnishes the corresponding stable detonation velocity. This detonation velocity accounts for the coupling between the hydrodynamics and the chemical kinetics. The kinetic parameter may be varied in order to reach the experimental detonation velocity. This type of process may be done for all front curvatures, that is for all charge diameters. When the computed detonation velocities and the experimental ones are in good agreement for any front curvature, the kinetic parameters are determined. This type of inverse procedure is detailed in [17] and [20]. 5.4. Numerical approach for the reactive model In the reactive case, the equation of state (43) is now a function of three variables. The mixture internal energy can be written as: ¼ e þ c. Thus the total energy equation ot E þ r ðqHuÞ ¼ 0 can be put under the form: ot E þ r ðqH uÞ ¼ qk_c0 0 with E ¼ qe þ qðu uÞ=2 and for any variable ðÞ ¼ dðÞ=dk. The contribution qk_c0 is related to the combustion energy release. The equations governing the material interfaces motion must also be modified. Indeed, in the reactive case, the material parameters b and h of Eq. (43) now depend on the volume fraction k. Therefore, the materials parameters are governed by the equations:
ot b þ u rb ¼ f k_b0
and ot h þ u rh ¼ f k_h0
ð44Þ
To summarize, in the reactive case, the overall system to solve reads: ot q þ rðquÞ ¼ 0 ot qu þ rðqu uÞ þ rp ¼ 0 ot E þ rðqH uÞ ¼ qk_c0 ot qf þ rðqf uÞ ¼ 0
ð45Þ
ot qk þ rðqkuÞ ¼ qk_f ot qb þ rðqbuÞ ¼ qf k_b0 ot qh þ rðqhuÞ ¼ qf k_h0 The overall system is split into a convective one and an ODE system. The convective part of this system is solved by the hyperbolic solver detailed in the preceding sections. The source terms are then accounted for by an appropriate ODE solver.
6. Numerical experiments The technique developed in this paper is now applied to three tests cases: a water–air shock tube, an interaction of a shock with an air bubble in water and a detonation wave. The com-
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
591
putations are performed with 2D unstructured meshes for the two first cases and with structured meshes for the last one. In any cases, the methodology developed in the framework of unstructured grids is applied and a second order accurate MUSCL explicit scheme is used. 6.1. Water–air shock tube We consider a square domain of 1 m 1 m length. Initially there is a material interface, located at x ¼ 0:7 m, separating two fluids: water at the left and air at the right. The two fluids are initially at rest and there is a jump on the density, the pressure and the material parameters (b and h). The solution of this Riemann problem is composed by a shock wave travelling in the air, a rarefraction travelling in the water and the material interface (contact discontinuity) [9]. The explicit second order accurate scheme used here is based on the classical RoeÕs solver with CFL ¼ 0:75. The computations of the first test case (Fig. 5) are performed on an unstructured mesh of 4506 vertices with about 750 vertices in the wave propagating direction. The numerical results are in good agreement with the analytical solution. For the sake of clarity, only the numerical solutions, at the final time T ¼ 230 ls, are plotted (Fig. 6). The 1D structure of the 1D analytical solution is preserved by this modified conservative approximation. Moreover, there are no oscillations at the material interface (Fig. 6). The solution with the coarse mesh (1782 vertices) is less accurate but the global profile of the waves are captured (Fig. 7). The solution of the second test case (Fig. 8) is, analytically, composed by the solution of the test case one and a strong contact discontinuity on the v component of the velocity. When a conservative scheme is used (Roe, Rusanov or HLL solvers), the numerical solution consist in a shock wave travelling faster than the exact one. This is a consequence of the lack of accuracy, on the tangential velocity (v), for the conservative schemes. Indeed, the pressure and the sound speed, used to numerically compute the waves speed, are directly related to the norm of the velocity. Computations performed on structured or unstructured meshes with first or
Fig. 5. Initialization of the water–air shock tube: Test 1.
592
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Fig. 6. Density, velocity pressure and EOS parameters at the final time 230 ls, RoeÕs flux and second order accurate MUSCL scheme on unstructured mesh.
second order accurate conservative schemes give rise to the same problem (Fig. 9). We have also computed some numerical solutions with the partially non-conservative scheme. Therefore, the numerical results, computed with this non-conservative scheme (Fig. 9), are very close to the analytical one.
6.2. Shock–bubble interaction In this section, we investigate the dynamics of a shock–bubble interaction. The size of the computational domain is 20 10 cm2 , with absorption boundary conditions on the right and piston conditions on the left. The initial flow is water at rest with an cylindrical air bubble, of 4 cm radius, centered at the point (5, 5) cm. The system is at the equilibrium with uniform pressure and velocity then a shock wave is produced at the left of the domain by the piston motion (Fig. 10). In this test case, an unstructured mesh of 98 489 nodes and 195 776 triangles is used. Numerical results are computed with an explicit second order accurate scheme based on the Roe solver and the CFL ¼ 0:5. In this computation, the material interface is not explicitly defined, hence we represent it by the isolines b ¼ 1:5. The speed of sound inside the air bubble is smaller than that of the surrounding water. After the first interaction with the bubble, the incident shock wave is diffracted into two waves (Fig. 11): a
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
593
Fig. 7. Comparison between the computations with the second order accurate MUSCL scheme for an unstructured coarse mesh of about 100 points in the propagation direction and an unstructured fine mesh of about 750 points in the propagation direction.
reflected rare-fraction wave moving from right to left in the water and a transmitted shock wave in the bubble. Later, when the incident wave have passed the bubble, we can observe the multiple scattering phenomena (Figs. 12 and 13). Due to the fact that air compressibility is larger than the water one, the downstream bubble interface is almost motionless during the compression (Fig. 16). We notice that the numerical scheme preserves this material interface and no spurious oscillations or non-physical states are produced during the computation. Then, the transmitted shock emerges in the water (Fig. 14 for t ’ 18:5 ms) and deforms the downstream material front. Two pairs of vortices (Figs. 15 and 16) are produced around the compressed air bubble. The upstream interface is folded and progressively a re-entrant jet is formed. This jet breaks the bubble in two parts and drives them downstream (Fig. 17). Even if the context is not the same, the behavior of the shock–bubble obtained in this investigation is similar to the experiments realized in [13].
594
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Fig. 8. Initialization of the water–air shock tube, Test case 2 (left) and numerical solution with the RusanovÕs flux, first order conservative scheme at the final time 230 ls (right): structured mesh of 705 3 nodes.
6.3. Detonation wave We consider a cylindrical Plexiglass tube filled with a solid explosive. The tube is surrounded by air. The explosive is a mixture of ammonium nitrate, fuel, water and other minor compounds. It is a non-ideal explosive: its chemical decomposition rate is very slow compared to standard solid high explosives. Consequently the consideration of the chemical kinetic effects are very important for the detonation propagation and its effects over structures. The configuration studied consists in several Plexiglass tubes of different diameter. The tube thickness is 3 mm. The initial density of Plexiglass is 1186 kg/m3 and the parameters of the equation of state are c ¼ 3:8 and p1 ¼ 2 109 Pa. The EOS parameters for the solid state of the explosive are: s0 ¼ 0, cs ¼ 2:9, p1s ¼ 4:5 108 Pa and the initial density is q0 ¼ 1000 kg/m3 . The EOS parameters of the detonation products are: cg ¼ 2:37, p1g ¼ 0, pCJ ¼ 11 109 Pa, DCJ ¼ 5875 m/s and qCJ ¼ 1477 kg/m3 . The CJ state is obtained by a thermochemical equilibrium computation using CHEETAH [8]. Finally, the kinetic parameter of Eq. (40) are a ¼ 3:34 103 Pan and n ¼ 0:8 when p > p0 and a ¼ 0 otherwise. The threshold pressure p0 is take equal to 1 GPa. The following computation is performed with a mesh of 300 75 cells and we first consider a cylindrical charge of 84 mm diameter (Fig. 18). The computational domain is 30 cm long and 6 cm high. The results are shown when the stable detonation state is obtained (Fig. 18). The detonation front is clearly visible as well as the expansion of the Plexiglass cylinder. These results are to compare with X-rays measurements as shown in Fig. 19. The qualitative agreement between the computations and experiment can be noticed. The quantitative comparison will be given after examining the effects of the charge size. We now consider another charge of 174 mm internal diameter. The mesh involves now 400 150 cells. The computed density contours are shown in Fig. 20. Again, the shock front and tube expansion are clearly visible. The numerical data are now compared with the experimental ones in Table 1. From these results an excellent agreement can be noticed. It is important to note the wide range of variation of the detonation velocity compared to the ideal one (5875 m/s). This is a result
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
595
Fig. 9. Comparison between the conservative and the non-conservative scheme. Numerical solution with the RusanovÕs flux, first order conservative scheme, at the final time 230 ls: structured mesh of 705 3 nodes.
Fig. 10. Initialization of the shock–bubble interaction test case.
596
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Fig. 11. Evolution of the reflected, the refracted and the incident shock waves. Times t ¼ 0:65, 1.30, 1.95 and 2.6 ms.
Fig. 12. Pressure lines at t ¼ 9:75 ms: refracted and diffracted shock waves.
of the coupling effects between the low rate of the chemical kinetic of this explosive, and the rarefraction waves travelling inside the reaction zone.
7. Conclusion In the case of the stiffened gas EOS, we have pointed out the different problems arising in the computation of the multi-dimensional multi-material flows: the pressure can become oscillatory across slip lines, and consequently the velocity and the other variables too. For very different EOS, and/or extreme flow configurations this leads to the breakdown of the code. The situation is even worse in several space dimensions because the kinetic velocity is not accurately computed around slip lines.
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
597
Fig. 13. Streamlines at t ¼ 9:75 ms: refracted and diffracted shock waves.
Fig. 14. Evolution of the pressure (p þ p1 ) lines. Times t ¼ 0, 6.5, 13.0, 18.5, 26.0 and 32.5 ms.
We have discussed these different problems and propose solutions. They amount to a slight modification of existing schemes; they are easily incorporated into existing codes. These
598
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Fig. 15. Evolution of the material interface (b ¼ 1:5). Times t ¼ 0, 6.5, 13.0, 18.5, 26.0 and 32.5 ms.
Fig. 16. Evolution of the material interface (b ¼ 1:5). Times t ¼ 0, 6.5, 13.0, 18.5, 26.0 and 32.5 ms.
modifications are of two types: first the identification of the numerically pertinent thermodynamical parameters that have to be updated, second more consistent calculation of the kinetic energy. We note that the second ingredient of the method arises because of the multi-dimensional nature of the problems and the class of conservative schemes (directional splitting or finite volume
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
599
Fig. 17. Vortex pairs developed around the bubble at t ¼ 26:0 ms.
Fig. 18. Detonation wave in a 84 mm diameter cylindrical charge in a Plexiglass tube. Density contours.
approach) we consider here. The drawback of our method is that the scheme is not any longer strictly conservative on the total energy. We note however that the loss of conservation occur only around slip lines and does not seem to pollute the solutions. The method is flexible, it is applicable on any type of mesh: Cartesian, unstructured, hybrid, because of a local formulation.
600
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Fig. 19. X-ray measurement of the shock front and cylinder expansion from Saurel et al. [23].
Fig. 20. Detonation wave in a 174 mm diameter cylindrical charge in a Plexiglass tube. Density contours.
Table 1 Experimental vs computational results for the propagation of detonation waves in explosive cylinders of various charge diameter. Experiments: left column. Computations: right column Charge diameter (mm)
Detonation velocity
Radius of curvature
Expansion angle
Exp. (m/s)
Comp. (m/s)
Exp. (mm)
Comp. (mm)
Exp.
Comp.
84 104 174
3060 3467 4000
3160 3369 4006
96 116 290
104 120 285
32 32 33
31 32 32
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
601
Several numerical examples, with and without the corrections, are presented. They demonstrate the efficiency of our method on triangular unstructured meshes. Acknowledgements This paper has been completed during the CEMRACS 2000 summer school, held at CIRM (Marseilles) during the summer 2000. The third author has also been partially supported by NITROCHIMIE. He addresses special thanks to S. Mencacci and T. Rousse. The experiments shown for the detonation problem have been done at CEA/DAM Moronvilliers and we particularly thank G. Rodriguez. This research has been conducted as a research project of the Laboratoire Correspondant du CEA LRC M03. Appendix A. Derivation of the Roe average for the stiffness gas EOS In the case of the stiffened gas, the equation of state is bp þ h ¼ q ¼ E 12qu2
ðA:1Þ
In the 1D case, the fluid is governed by: ot W þ ox F ¼ 0
ðA:2Þ
with initial and boundary conditions. The vector of conserved variables is W ¼ ðq; qu; E; qb; qhÞ :¼ ðq; m; E; B; HÞ and the flux is F ðW Þ ¼ ðqu; qu2 þ p; uðE þ pÞ; uqb; uqhÞ. The Jacobian of F with respect to W is 1 0 0 1 0 0 0 C B C B u2 þ pq 2u þ p p p p m E B H C B B uðH þ pÞ H þ upm uð1 þ pE Þ upB upH C ðA:3Þ C B C B ub b 0 u 0 A @ uh h 0 0 u where þ pq ¼
b
u2 2
pm ¼
u b
pE ¼
1 b
pB ¼
p qb
pH ¼
1 qb
ðA:4Þ
If we set Df ¼ fL fR , the Roe average is defined by replacing in (A.3) the Roe average of the velocity and the total enthalpy H and by defining the average derivative of the pressure such that the following jump relation is satisfied, e f e f 1 p 1 Dq þ DðqÞ DðqbÞ DðqhÞ ¼ Dp b b qb qb
ðA:5Þ
Let k 2 ½0; 1 and define for any function f, f ¼ kfL þ ð1 kÞfR and f ¼ kfR þ ð1 kÞfL . For any functions f and g, we have
602
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Dfg ¼ f Dg þ gDf pffiffiffiffiffi pffiffiffiffiffi pffiffiffiffiffiffi Last, if k ¼ qL =ð qL þ qR Þ, we have qf ¼ qf . Using the above calculation rules, we can rewrite Dp as q h 1 1 1 Db Dh Dp ¼ qD þ Dq 2 b b b b where we have set b2 ¼ bL bR . Last, by proceeding by identification, we get ! f b f e 1 e 1 1 p 1 q h 1 1 :¼ :¼ ¼ 2 :¼ :¼ 2 b b b b b qb qb qb q b The matrix of eigenvectors is easily computed, 0 1 1 B u B 2 u 2 =2Þ ðu =2Þ ðu R ¼ ð r1 r2 r3 r4 r5 Þ ¼ B B @ b 0 h 0
1 1 1 1 u uc uþc C C 2 ðu =2Þ H uc H þ uc C C 0 b b A q h h associated to the eigenvalues u, u, u, u c, u þ c where c2 ¼ b H ðu2 =2Þ . Similarly, any vector W ¼ ðA; B; C; D; EÞ in R5 can be decomposed as: W ¼
5 X
ðA:6Þ
ðA:7Þ
li ðW Þri
i¼1
We get, for W ¼ ðDq; Dm; DE; DB; DHÞ, 1 bDh hDb Dp l2 ðW Þ ¼ Dq 2 l1 l3 c h q Dh Db l3 ðW Þ ¼ b q 1 Dp l4 ðW Þ ¼ Du 2 c2 c 1 Dp q l5 ðW Þ ¼ þ Du 2 c2 c l1 ðW Þ ¼
Appendix B. 2D Roe type scheme This section is devoted to the derivation of Eq. (30) in Section 4.5. We first introduce some notations. We set
ðA:8Þ
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
1 W B qV C C B Z¼B C; @ qb A 0
1 F B qVU C C B G¼B C; @ qbU A 0
qhU
qh and
0
1 B B0 B Mij ¼ B B0 @ 0
0 gx 0 g? x
0 gy 0 g? y
1 0 C 0C C C; 1C A 0
Z ¼ Mij W ;
603
W ¼ M1 ij Z
1 S0 B US0 þ S1 C C B B H S0 þ US1 C C Jðu; H; b; h; S0 ; S1 Þ ¼ oZ G ¼ UId þ B C B VS0 C B A @ bS0 hS0 0
where S0 ¼ ðU
1 0
0
S1 ¼ ð oq p
oqU p
oE p
0
0 Þ ¼ qrZ U oqV p
oqb p
oqh p Þ ¼ rZ p
We denote by Zi and Zj the values of the conservative variables in the cells i and j. For any function f, we set f ¼ kfi þ ð1 kÞfj and f ¼ kfj þ ð1 kÞfi . These two averages have the property that DðfgÞ ¼ f Dg þ gDf Moreover, if pffiffiffiffi qi k ¼ pffiffiffiffi pffiffiffiffiffi qi þ qj simplify the following calculations. we have qg ¼ qg: These remarks e u; H ; b; h; S ; S with q~ ¼ bp þ h, We define the matrix J 0 1 S0 ¼ U 1 0 0 0 0 p 1 e U 1 V ~ þ Q S1 ¼ q b
1 q
e satisfy the jump condition DG ¼ J e DZ (DZ ¼ Dij Z ¼ e ¼ ðU2 þ V2 Þ=2. The matrix J where Q Zi Zj ). Moreover, we have the relations: S0 DZ ¼ qDU;
S1 DZ ¼ Dp
604
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
e is diagonalisable and the eigenvectors are The matrix J 1 1 0 1 0 0 1 1 1 BUC B U c~ C B U þ c~ C C C B C B B C C BeC B B BQC B H c~U C B H þ c~U C C C C B B r1 ¼ B B V C; r2 ¼ B V C; r3 ¼ B V C; C C B C B B C C B C B B A A @bA @ @ b b h 0 1 0 B0C B C B C B ~ C C r5 ¼ B B 0 C; B C B C @bA
0 1 0 B0C B C B C B1C C r6 ¼ B B0C B C B C @0A
h
0
0
1
B 0 C B C B C BVC C r4 ¼ B B 1 C; B C B C @ 0 A 0
h
q
h
e . The associated eigenvalues are where b~ c2 ¼ H Q k1 ¼ U;
k2 ¼ U c~;
k3 ¼ U þ c~;
k4 ¼ U;
k5 ¼ U;
k6 ¼ U
In the eigenvectors basis, the jump of conservative variables is given by: DZ ¼ ðDðqÞ nÞr1 þ þ
nH nþH DðqbÞ bDq r5 r2 þ r3 þ DðqVÞ VDðqÞ r4 þ 2 2 b
bDðqhÞ hDðqbÞ r6 qb
where 1 DðqbÞ bDq bDðqhÞ hDðqbÞ VDðqVÞ UDðqUÞ þ QDðqÞ þ DðEÞ ~ n¼ 2 b~ c b qb
ðA:9Þ !
Dp DðqUÞ UDðqÞ and H ¼ 2 c~ c~ After some tedious but straightforward calculations, we get ¼
n¼
Dp c~2
and H ¼ q
DU c~
References [1] Abgral R. Generalisation of the Roe scheme for the computation of a mixture of perfect gases. La Recherche Aerospatiale 1988;6:31–43. [2] Abgrall R. How to prevent pressure oscillations in multicomponent flows: a quasi conservative approach. J Comput Phys 1996;125:150–60.
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
605
[3] Abgrall R, Kami S. Compressible multifluids. J Comput Phys 2001;169:594–623. [4] Chargy D, Abgrall R, Fezoui L, Larrouturou B. Conservative numerical schemes for multicomponent inviscid flows. La Recherche Aerospatiale 1992;2:61–79. [5] Chern I-L, Glimm J, McBryan O, Plohr B, Yaniv S. Front tracking for gas dynamics. J Comput Phys 1986;62:83– 110. [6] Fedkiw RK, Aslam T, Merriman B, Osher S. A non oscillatory eulerian approach to interfaces in multimaterials flows (the ghost fluid method). J Comput Phys 1999;152:457–92. [7] Fickett W, Davis WC. Detonation. Technical report, University of California Press, Berkeley, 1979. [8] Fried E. CHEETAH 1.39 userÕs manual. Lawrence Livermore National Laboratory, UCRL-MA-117541 Rev 3, Energetic Material Center, 1996. [9] Godlewski E, Raviar PA. Numerical approximation of hyperbolic systems of conservation laws. In: Applied Mathematical Sciences, vol. 118. New York: Springer; 1995. [10] Godunov SK. A finite difference method for the numerical computation of discontinuous solutions of the equations of fluid dynamics. Mat Sb 1959;47:357–93. [11] Godunov SK, Zabrodine A, Ivanov M, Kraiko A, Prokopov G. Resolution numerique des probl emes multidimensionnels de la dynamique des gaz. Moscow: Editions Mir; 1979. [12] Hartlow F, Amsden A. Fluid dynamics. Technical Report LA-4700, Lawrence Livermore National Laboratory, 1971. [13] Hass J-F, Sturtevant B. Interaction of a weak shock waves with cylindrical and spherical gas inhomogeneities. J Fluid Mech 1987;181:41–76. [14] Karni S. Hybrid multifluid algorithms. SIAM J Sci Comput 1996;17:1019–39. [15] Karni S. Multi-component flow calculations by a consistent primitive algorithm. J Comput Phys 1996;112:31–43. [16] Massoni J, Saurel R, Baudin G, Demol G. A mechanistic model for shock initiation of solid explosives. Phys Fluids 1999;11:710–36. [17] Persson PA, Holmberg R, Lee J. Rock blasting and explosive engineering. Boca Raton: CRC Press; 1996. [18] Roe PL. Riemann solvers, parameter vectors and difference schemes. J Comput Phys 1981;43:357–72. [19] Rusanov VV. Calculation of interactions of non steady shock waves with obstacle. J Comput Math Phys 1961;1:267–79. [20] Saurel R. DPEX: Un outil numerique pour la determination des cinetiques reactionnelles des explosifs fortement non ideaux. Rapport de contrat, Centre dÕEtudes de Gramat, IUSTI, Universite de Provence, Marseilles, France, 1996. [21] Saurel R, Abgrall R. A multiphase Godunov method for compressible multifluid and multiphase flows. J Comput Phys 1999;150:425–67. [22] Saurel R, Abgrall R. A simple method for compressible multifluid flows. Siam J Sci Comp 1999;21:1115–45. [23] Saurel R, Mencacci S, Mercier S. Hydro3d: Un outil simple pour la simulation de probl emes de detonique. In: EuroPyro 99. 7 eme congr es international de pyrotechnic du GTPS, 1999. Brest. [24] Serre D. Systems of conservation laws I: Hyperbolicity, entropies, shock waves vol. XXII. Cambridge: Cambridge University Press; 1999, Translated from the French by I.N. Sneddon. [25] Shyue KM. An efficient shock capturing shock––capturing scheme for multicomponent problems. J Comput Phys 1998;142:208–42. [26] Toro EF. Some ivps for which conservative methods fail miserably. In: Proceedings of the 6th International Symposium on Computational Fluid Dynamics, Lake Tahoe, California, 1995. p. 4–8. [27] Yee HC. Construction of explicit and implicit symmetric tvd schemes and their applications. J Comput Phys 1987;68:151–79.
Efficient numerical approximation of compressible multi-material flow for unstructured meshes R. Abgrall a
a,b,* ,
B. Nkonga a, R. Saurel
c
Math ematiques Appliqu ees de Bordeaux, Universit e Bordeaux I, 351 Cours de la Liberation, 33405 Talence cedex, France b Institut Universitaire de France, France c IUSTI, Universit e Aix Marseille I, 13453 Marseille Cedex 13, France
Received 19 October 2000; received in revised form 31 January 2002; accepted 4 February 2002
Abstract We consider the numerical approximation of multi-dimensional multi-material flows. This is a difficult topic related to the numerical smearing of contact discontinuities (or material interfaces or slip lines). Any Eulerian scheme will produce an artificial mixing zone. In this artificial mixture, the computation of thermodynamical variables (pressure, sound speed, temperature, . . .) is difficult to achieve correctly. For the stiff cases considered in this paper, with solid–liquid–gas interfaces, small errors on the thermodynamical variables lead to the blow up of the computation in many cases. In this paper, we review and explain this problem, then we propose solutions for 1D and 2D flows. Contrarily to front-tracking techniques, our schemes use the same formulation everywhere on the mesh. We provide several examples with the stiffened gas equation of state (EOS): inert materials (water–air), and chemically reactive materials (solid explosive– Plexiglass–air). 2003 Elsevier Science Ltd. All rights reserved.
1. Introduction and physical model The computation of multi-material flows is more and more investigated because many physical situations involve interfaces and mixtures of several materials. For example, bubble-gas problems, interface problems, problems involving jets of different physical properties in rocket engines, etc. * Corresponding author. Address: Mathe´matiques Appliquees de Bordeaux, Universite Bordeaux I, 351 Cours de la Liberation, 33405 Talence cedex, France. Tel.: +33-5-56-84-60-68; fax: +33-5-56-84-26-26. E-mail addresses: [email protected] (R. Abgrall), [email protected] (B. Nkonga), richard@ iusti.univ-mrs.fr (R. Saurel).
0045-7930/03/$ - see front matter 2003 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 3 0 ( 0 2 ) 0 0 0 1 2 - 9
572
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
There exists a large literature about the numerical discretization of these problems: Eulerian schemes, Lagrangian schemes, front tracking techniques [5], level set methods [6], etc. We are interested in the development of Eulerian schemes where the equation of state (EOS) is described by parameters (such as the ratio of specific heats) that may vary within the computational domain, but the analytical formulation is the same everywhere. We also want to compute flows by using the same numerical formulation everywhere and by introducing few modifications of existing conservative schemes in order to make easy the evolution of existing codes. One of our goals is to construct schemes that can be modified easily to take into account creation of species by chemical processes, or able to consider viscous phenomena. This is not possible by the Eulerian schemes using the level set formalism, for example. In this paper, we consider gases modeled by the stiffened gas equation of state [11,12], but our technique could easily be adapted to more complex cases following for instance [25]. The stiffened gas EOS is a modification of the ideal gas EOS. It enables to consider complex situations such as compressible liquid–gas problems. More precisely, we consider problems where the pressure is related to the internal energy by: p ¼ ðcðx; tÞ 1Þq cðx; tÞp1 ðx; tÞ
ð1Þ
where the thermodynamical parameters c and p1 depend on time and space. This formulation permits to consider problems where the thermodynamical parameters are piece-wise constant. More complex situations, where the EOS is of the Mie Gr€ uneisen form, can also be considered by freezing in each computational cell the thermodynamical parameters. ‘‘Naive’’ extensions of the classical conservative schemes such as Roe scheme, often lead to severe numerical problems. First, the thermodynamical parameters may reach values out of their range of validity, and second, the numerical results appear to be oscillatory and non-physical at interfaces (negative pressure, negative arguments in the computation of the sound speed). This is not a drawback of a particular scheme, but a consequence of the conservative formulation of the schemes [1–4,15]. Many investigations have been done to define an efficient numerical approximation for multimaterial flows, see for example [1,3,14] for a review. Most of the time, these oscillations are barely visible, but in some cases (large pressure or density ratio at interfaces, very different EOS), they can lead to the blow up of the computation. Therefore, the major challenge is the numerical characterization of the material properties, in the mixing zone and the handling of pressure oscillations across material interfaces. The previous investigations done in this topic are always based on 1D analysis or in several dimensions but for very regular meshes [21]. In this basically 1D techniques, the system of PDEs is discretized in such a way that the mass, momentum and total energy are locally conserved, hence the correct shock speed can be computed. The flow composition is updated such that, in the case of a perfect contact discontinuity, the velocity and the pressure remain uniform. However, this method cannot be easily generalized to several dimensions, especially when the mesh is an unstructured or a very distorted structured mesh, because there are no natural direction in the mesh. In this paper, we show how to extend the method of [2,21] to general unstructured meshes. A careful analysis is conducted and enables to alleviate the remaining difficulties contained in [2] about high order extensions, and to proceed the multi-dimensional case easily. More precisely, we first discuss our model, then consider the continuous case to illustrate the general principles of the
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
573
paper in one and several dimensions. The next section is devoted to the illustration of these principles to 1D case (the Roe, Godunov and Rusanov schemes, with a particular emphasis on the Roe scheme); then, we discuss the multi-dimensional case, and show which additional modifications have to be made to take properly into account the physics of slip lines. We conclude by showing how the method can be used for flows with chemistry, and different test cases, inert and reactive, are discussed.
2. The physical model We consider a flow which dynamic is governed by the compressible Euler equations. In this context, the conservative formulation of the multi-fluid flow reads: ot q þ r ðquÞ ¼ 0 ot m þ r ðm uÞ þ rp ¼ 0
ð2Þ
ot E þ r ðqH uÞ ¼ 0 where q is the density, u the velocity, m ¼ qu the momentum, E the total energy, p the hydrodynamic pressure and H ¼ ðE þ pÞ=q the total enthalpy. This system is generally written in the following form: ot W þ divðFÞ ¼ 0 The pressure is related to the other variables by the equation of state (1). For simplicity reasons, we rewrite the ‘‘stiffened’’ gas EOS as: bp þ h ¼ E 12qu u where b :¼
1 c1
and h :¼
cp1 c1
are the thermodynamical parameters. Their evolutions are governed by the equations: ot b þ u rb ¼ 0;
and ot h þ u rh ¼ 0
ð3Þ
We note that an equivalent conservative form of the equations can be written, ot q þ r ðquÞ ¼ 0 ot ðquÞ þ r ðqu uÞ þ rp ¼ 0 ot E þ r ð qH uÞ ¼ 0
ð4Þ
ot ðqbÞ þ r ðqbuÞ ¼ 0 ot ðqhÞ þ r ðqhuÞ ¼ 0 Other equivalent systems could be obtained. We rely on this particular form to explain our technique, even though the other systems might also have been considered instead. The system (2), (3) or (4) is hyperbolic. In the 1D case, the eigenvalues are u, u c and u þ c where the sound speed is defined by:
574
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
c2 ¼
cðp þ p1 Þ q
Across a contact discontinuity, the pressure and the velocity remain constant. In the multidimensional case, the pressure and the normal velocity are continuous across a slip line. 3. General principles We begin with a description of the governing principles of our technique on system (2) in 1D. Then we recall some known problems in several dimensions; they already exist for perfect gas computations, but are more pronounced in the present one. We describe how to cure them. The application of these techniques, for 1D and 2D schemes, are given in Section 4.5. 3.1. The one-dimensional case Consider Dx > 0, xi ¼ iDx for i 2 Z, Dt > 0 and the conservative scheme dq þ mDFq ¼ 0 dqu þ mDFqu ¼ 0 dE þ mDFE ¼ 0
ð5Þ
In (5), we make use of the notation, for any function f, df : ¼ f nþ1 f n , Df : ¼ fiþ1=2 fi1=2 and m ¼ Dt=Dx. The question is to find a discretization of the non-conservative system (3) compatible with the physics of (2) and (3), namely the velocity and the pressure should stay continuous across material interfaces. It is well known that the smooth solutions of (2) are solutions of ot q þ uox q þ qox u ¼ 0
ð6Þ
qot u þ quox u þ ox p ¼ 0
ð7Þ
ot p þ uox p þ qc2 ox u ¼ 0
ð8Þ
The last equation is equivalent to the transport of physical entropy. The system (6)–(8) can be obtained from (2) and (3) by algebraic manipulations. The first equation is obtained by developing ot q þ ox qu. The second one is uðot q þ ox quÞ þ ot qu þ ox qu2 þ ox p ¼ 0 and the last one is obtained in two steps: • first u2 ot q þ uox q þ pox u ¼ ð ot E þ ox ½uðE þ pÞ Þ u ot qu þ ox qu2 þ ox p þ ðot q þ ox quÞ 2 • second, since q ¼ bp þ h, ot q þ uox q þ pox u ¼ pðot b þ uox bÞ þ ðot h þ uox hÞ þ b ot p þ uox p þ qc2 ox u
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
575
Our technique is to mimic these algebraic manipulations and to identify discrete operators that are consistent with the continuous operators between the brackets above. We have dE ¼ dðqÞ þ 12dðqu2 Þ. The idea is to write a scheme for the internal energy by combining the equations for E, qu and q. For an explicit scheme, dðqu2 Þ ¼ ðquÞ
nþ1
du þ un dqu
ð9Þ
dðquÞ ¼ qnþ1 du þ un dq and then, dðqÞ ¼ dE
un þ unþ1 un unþ1 dqu þ dq 2 2
ð10Þ
that is un þ unþ1 un unþ1 dðbp þ hÞ þ m DFE DFqu þ DFq ¼ 0 2 2
ð11Þ
We rewrite the flux variation of any Q-scheme [9], DF ¼ Fjþ1=2 Fj1=2 as DF ¼ Fjþ1=2 Fj1=2 ¼ 12 Wjþ1=2 þ Wj1=2 where Wjþ1=2 ¼ Fjþ1 Fj QðWjþ1 Wj Þ . In the sequel, we often call Wjþ1=2 a fluctuation. The fluctuation form of the scheme is DF ¼ 12 Wjþ1=2 þ Wj1=2 T
n nþ1 Denoting by ‘j the vector ‘j ¼ ððunj unþ1 j Þ=2; ðuj þ uj Þ=2; 1Þ , Eq. (11) can be rewritten as
dðbp þ hÞ ¼ dðbp þ hÞjþ1=2 þ dðbp þ hÞj1=2 where T dðbp þ hÞjþ1=2 ¼ 12‘j Wq ; Wqu ; WE jþ1=2 T dðbp þ hÞj1=2 ¼ 12‘j Wq ; Wqu ; WE j1=2 Then, we can rewrite, for an explicit scheme, dðbp þ hÞjþ1=2 as: dðbp þ hÞjþ1=2 ¼ pjn dbjþ1=2 þ dhjþ1=2 þ bnþ1 j dpjþ1=2 T The last step is to identify in 12‘j Wq ; Wqu ; WE jþ1=2 the terms that are consistent with n bnþ1 j dpjþ1=2 , pj dbjþ1=2 and dhjþ1=2 in the limit of a mesh refinement. A guide is that the evolution scheme for p should preserve the pressure in the case of a contact discontinuity, as in [2,22], so the idea is to write the fluctuations in term of variations of the primitive variables, in particular the velocity and the pressure. This can be done very naturally, because all the known schemes have, up to our knowledge, their natural formulation in term of these variables. Thus, it becomes easy to identify the sets of variations that are consistent with the density, velocity and pressure evolution equations. Since the model involves thermodynamical parameters, after having identified
576
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
the previous sets of variations, the remaining terms describe the discrete evolution of these thermodynamical parameters: this provides the fluctuations needed to update b and a. In the Sections 4.1, 4.2 and 4.3, we illustrate this technique for RoeÕs, GodunovÕs and RusanovÕs schemes and, in Section 4.4, for a second order MUSCL type scheme. We follow exactly these principle to make things clear. Remark 3.1. For an implicit scheme, we would have switched the superscript n þ 1 and n in (9) and made the modification accordingly in (10). 3.2. The multi-dimensional case In several space dimensions, the situation is a bit more complex, because slip lines are intrinsically ill represented in finite volume schemes, even though the defaults are generally weak. These problems appear first on the pressure, and then on the other variables, even in the perfect gas case, see [26]. Hence, we have to correct several problems: a multi-dimensional effect, and a problem due to the equation of state. In order to explain this, let us recall the partial differential equation (2) defining the behavior of in-viscid smooth flow, ot q þ r ðquÞ ¼ 0 ot m þ r ðm uÞ þ rp ¼ 0 ot E þ r ðqH uÞ ¼ 0 If C is any sub-cell of X, then for non-smooth solutions of Eq. (2) becomes: Z Z d W dx þ F n dr ¼ 0 dt C oC where F is the flux. This formulation suggests the finite volume scheme, where the building block is the solution, exact or approximate, of a Riemann problem, F ¼0 ot W þ on ~ ?
W ðg ¼ 0; g ; t ¼ 0Þ ¼
WL WR
if g < 0 if g > 0
ð12Þ
where ðg; g? Þ are the coordinates in the direct frame ðg; g? Þ. We have kgk ¼ 1 and g? is the unit vector normal to g at M 2 oC. We introduce U¼ug
and V ¼ u g?
This system reads: ot q þ og ðqUÞ ¼ 0 ot qu þ og ðquUÞ þ gx og p ¼ 0 ot qv þ og ðqvUÞ þ gy og p ¼ 0 ot E þ og ðqUHÞ ¼ 0
ð13Þ
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
577
and by simple manipulations, we get ot q þ og qU ¼ 0 ot qU þ og ðqU2 Þ þ og p ¼ 0
ð14Þ
ot qV þ og ðqUVÞ ¼ 0 ot E þ og ðqUH Þ ¼ 0
In practice, the numerical fluxes are constructed from (13) and not (14) because of conservation issues. However, this has some drawbacks. As pointed out in [26], the computations of slip lines can be inaccurate. If one consider for example a slip line initially aligned with the mesh, like on Fig. 1, where the normal velocity is non-zero, then even for a perfect gas, it is easy to see that for any up-wind solver, the pressure will not remain uniform. This phenomena is quite similar to what happens for mixture of perfect gases and is due to the non-correct computation of the tangential velocity. Consider for example a supersonic slip line, in the situation sketched in Fig. 1. With obvious notations (the slip line is initially located at i ¼ 0), and for a positive velocity u, we have n n n ¼ q mu q q ð15Þ qnþ1 1;j 1;j 1;j 0;j ðquÞnþ1 1;j
¼
ðquÞn1;j
mu
2
qn1;j
qn0;j
ð16Þ
nþ1 n n n ðqvÞ1;j ¼ ðqvÞ1;j mu ðqvÞ1;j ðqv Þ0;j
ð17Þ
nþ1 n n n ¼ E1;j mu E1;j E1;j E0;j
ð18Þ
From (15) and (16), we get unþ1 1;j ¼ u: the x-component of velocity remains uniform. From (15) and (17), the y-component of the velocity at ð1; jÞ lies in between of v0;j and v1;j under the CFL condition um 6 1, namely
Fig. 1. Slip line configuration.
578
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
vnþ1 1;j ¼
qn1;j ð1 muÞvn1;j þ qn0;j muv0;j qn1;j ð1 muÞ þ qn0;j mu
ð19Þ
From (18), using the perfect gas EOS and the uniformity of ui;j at times tn and tnþ1 , and the uniformity of the pressure at time tn , we get 2 nþ1 n p1;j ðvn1;j Þ p1;j ðvnþ1 u 1;j Þ þ ¼ þ m v21;j v20;j 2 c1 2 c1 2 2
In order to recover the uniformity of the pressure, we need 2 ðvn1;j Þ ðvnþ1 u 1;j Þ ¼ m v21;j v20;j 2 2 2 2
which is generally incompatible with (19). This pollutes the evaluation of the kinetic energy, and then that of the pressure. In the case of mixture of gases or the stiffened equation of state, this phenomena may lead to the blow up of the calculation, for example when vL vR is large, see [22] for example. The scheme proposed in [22] used a Cartesian mesh and the multi-dimensional scheme uses a splitting technique. A solution to this problem is to follow the same ideas as in 1D. We write the conservative scheme over one time step in fluctuation form. At one cell interface, the fluctuations depend locally on the variations of the density, the pressure, the normal velocity to the cell interface, and the variation of tangential velocity to the cell interface. In fact, the tangential velocity plays a special role. Like the normal velocity, it appears also through its square (in the kinetic energy). But contrarily to the normal velocity, the variation of the tangential kinetic energy does not disappear, by essence, through a slip line. Thus, we have two type of discrete equations involving the tangential velocity: one on the tangential velocity, and one on the tangential kinetic energy, like in the simplified example above. Thus there is an incompatibility. The idea is, like in [22], to introduce a special discrete equation on the tangential kinetic energy which is consistent with the continuous one, namely ð20Þ ot qV2 þ Uox qV2 ¼ 0 that can easily be obtained. At the beginning of every time step, we initialize the variables of this equation with the values obtained at the end of the previous time step. Then we compute the tangential kinetic energy with the discrete equation. We use this tangential kinetic energy to accurately compute the pressure (or the internal energy), and we proceed through the next time step. The discrete equation on the tangential kinetic energy is only an auxiliary equation. This method is the one developed in [22] for Cartesian meshes. In the case of a general mesh, the principles are the same, but the realization is a bit more subtle. This is developed in Section 4.5. 3.3. Additional comments We have presented here design principles only. There are many possible ways of modifying a scheme according to our technique, and this will become clear in the examples we develop in the next section.
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
579
However, our experience shows that the best modified schemes are those which satisfy additional constraints arising from the maximum principle. Eq. (3) are transport equations, thus b and h satisfy a maximum principle. This is also true for the tangential kinetic energy from (20). For these reasons, we look for schemes on b, h and the kinetic tangential energy that satisfy a discrete maximum principle.
4. Illustration on some schemes In this section, we illustrate the design principles of Section 3 to several numerical fluxes: the Roe scheme, the Godunov scheme the Rusanov scheme. We also explain how to apply this method to get high order accurate schemes in time. 4.1. Example of Roe’s scheme We first recall some basic facts about RoeÕs scheme [18]. The average Jacobian matrix JðWj ; Wjþ1 Þ satisfies the jump condition: F ðWiþ1 Þ F ðWi Þ ¼ JðWj ; Wjþ1 ÞðWiþ1 Wi Þ The Jacobian A of the flux F depends on the velocity u, the total enthalpy H and additional parameters, namely the derivative of the pressure with respect to the density and the internal energy, pq , p . To obtain an average Jacobian matrix that reduces to the standard one for a perfect gas, it is useful to note that it should depend on the Roe average of the velocity, the Roe average of the total enthalpy, and the following jump relation should be satisfied, f pm Djþ1=2 m þ f pE Djþ1=2 E þ pf Djþ1=2 p ¼ peq Djþ1=2 q þ f qb Djþ1=2 ðqbÞ þ p qh Djþ1=2 ðqhÞ
ð21Þ
There are many possible choice for peq , pe , etc. What is important is that at the end, the characteristic decomposition of Djþ1=2 W :¼ Wjþ1 Wj , in terms of the eigenvectors of JðWj ; Wjþ1 Þ, is very simple and only depend on averaged quantities such as the velocity, the enthalpy, the jþ1=2 , cjþ1=2 , etc. The calculations are given in speed of sound, etc., which are denoted by ujþ1=2 , H Appendix A, as well as all the components of the averaged matrix and the right and left eigenvectors. The numerical diffusion is controlled by the function u; for the Roe scheme without entropy correction uðuÞ ¼ juj. Setting u ðX Þ ¼ ðX uðX ÞÞ=2, Eq. (10) can be written explicitly as 1 0 D p Djþ1=2 p jþ1=2 A djþ1=2 q ¼ dbpjþ1=2 þ djþ1=2 h ¼ Ajþ1=2 u ðujþ1=2 Þ@Djþ1=2 q 2 Bjþ1=2 2 cjþ1=2 cjþ1=2 Cjþ1=2 qjþ1=2
Djþ1=2 u þ Djþ1=2 mu ðujþ1=2 Þpi Djþ1=2 b mu ðujþ1=2 ÞDjþ1=2 h cjþ1=2
ð22Þ
580
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
with Ajþ1=2 ¼ mð u2jþ1=2 2 ujþ1=2 e q jþ1=2 Þ u jþ1=2 þ 2e jþ1=2 þ ujþ1=2cjþ1=2 Þ e ujþ1=2 þ cjþ1=2 ÞÞ ðH u jþ1=2 ðujþ1=2 þ cjþ1=2 Þ þ e q jþ1=2 2Bjþ1=2 ¼ mu ðð jþ1=2 ujþ1=2cjþ1=2 Þ e þ mu ðð ujþ1=2 cjþ1=2 ÞÞ ðH u jþ1=2 ðujþ1=2 cjþ1=2 Þ þ e q jþ1=2 jþ1=2 þ ujþ1=2cjþ1=2 Þ e 2C ¼ mu ðð ujþ1=2 þ cjþ1=2 ÞÞ ðH u jþ1=2 ðujþ1=2 þ cjþ1=2 Þ þ e q jþ1=2 jþ1=2 ujþ1=2cjþ1=2 Þ e ujþ1=2 cjþ1=2 ÞÞ ðH u jþ1=2 ðujþ1=2 cjþ1=2 Þ þ e q jþ1=2 mu ðð Djþ1=2 ¼ mu ð ujþ1=2 Þðpjþ1=2 pi ÞDjþ1=2 b where e u jþ1=2 ¼
unj þ unþ1 j ; 2
pjþ1=2 ¼
pffiffiffiffiffi pffiffiffiffiffiffiffiffiffi qj pjþ1 þ qjþ1 pj pffiffiffiffiffi pffiffiffiffiffiffiffiffiffi ; qj þ qjþ1
and
e q jþ1=2 ¼
unj unþ1 j 2
qjþ1=2 ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffi qj qjþ1
The terms djþ1=2 q þ dj1=2 q with Djþ1=2 p Djþ1=2 p Djþ1=2 u þ þ Bjþ1=2 þ Cjþ1=2 q dj1=2 q ¼ Ajþ1=2 u Djþ1=2 q þ Djþ1=2 c2 c2 c are consistent with upx þ qc2 ux , i.e. in the limit of a mesh refinement, and for smooth solution, we get upx þ qc2 ux . n We also have dðbpÞjþ1=2 ¼ bnþ1 j dp þ pj djþ1=2 b, so that a natural discretization is dbjþ1=2 ¼ mu ð ujþ1=2 ÞDjþ1=2 b dhjþ1=2 ¼ mu ð ujþ1=2 ÞDjþ1=2 h dbj1=2 ¼ muþ ð uj1=2 ÞDj1=2 b dhj1=2 ¼ muþ ð uj1=2 ÞDj1=2 h
ð23Þ
and dbj ¼ dbjþ1=2 þ dbj1=2 dhj ¼ dhjþ1=2 þ dhj1=2 It clearly prevents oscillations in the pressure (the contacts are perfectly described). The scheme (23) also admits a maximum principle on b and h under the CFL-like condition: m max u ð ujþ1=2 Þ 6 1 j
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
581
Fig. 2. Solution of the Riemann problem between states Wj and Wjþ1=2 . Wjþ1=2 is the state on the line x ¼ xjþ1=2 .
Moreover, the conservation of the density, momentum and total energy is guaranteed by construction. 4.2. Example of Godunov scheme Consider FðWj ; Wjþ1 Þ ¼ F ðWjþ1=2 Þ the flux for Godunov scheme [10] for (4) between the left and right states Wj and Wjþ1 (see Fig. 2). As we have seen above, the important quantities are not the flux but the residuals Wjþ1=2 ¼ FðWjþ1 ; Wj Þ F ðWj Þ :¼ F ðWjþ1=2 Þ F ðWj Þ; Wj1=2 ¼ F ðWj Þ FðWj ; Wj1 Þ :¼ F ðWj Þ F ðWj1=2 Þ
Thanks to the particular structure of GodunovÕs scheme, we can introduce Jjþ1=2 (resp. Jj1=2 ) and Wj (resp. Wj1=2 and Wj ). Calling the averaged entries the Roe average between the states Wjþ1=2 f of Jjþ1=2 (resp. Jj1=2 ) by f jþ1=2 and f jþ1=2 (resp. f j1=2 and f j1=2 , and setting Djþ1=2 ð:Þ ¼ ð:Þjþ1=2 ð:Þj (resp. Dj1=2 ð:Þ ¼ ð:Þj ð:Þj1=2 ), we can use the formulas (23) with the correct choice of Roe average to define the scheme, namely djþ1=2 b ¼ mujþ1=2 Djþ1=2 b djþ1=2 h ¼ mujþ1=2 Djþ1=2 h dj1=2 b ¼ muj1=2 Dj1=2 b
ð24Þ
dj1=2 h ¼ muj1=2 Dj1=2 h This scheme admits a maximum principle on b and h under a CFL condition. Mass, momentum and energy are also conserved.
582
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
4.3. Example of the Rusanov scheme The RusanovÕs flux [19] reads: F ðWi ; Wjþ1 Þ ¼ 12ðF ðWj Þ þ F ðWjþ1 Þ lðWjþ1 Wj Þ
so that the residuals can be written as: wjþ1=2 ¼ 12 Jjþ1=2 ljþ1=2 Id ðWjþ1 Wj Þ wj1=2 ¼ 12 Jj1=2 þ lj1=2 Id ðWj Wj1 Þ so that once more the formula (23) can be used with the functions u ð:Þ and uþ ð:Þ redefined by: u ðujþ1=2 Þ ¼ ðujþ1=2 ljþ1=2 Þ=2 and uþ ðuj1=2 Þ ¼ ðuj1=2 þ lj1=2 Þ=2, namely ðujþ1=2 ljþ1=2 Þ Djþ1=2 b 2 ðujþ1=2 ljþ1=2 Þ Djþ1=2 h djþ1=2 h ¼ m 2 ðuj1=2 þ lj1=2 Þ dj1=2 b ¼ m Dj1=2 b 2 ðuj1=2 þ lj1=2 Þ Dj1=2 h dj1=2 h ¼ m 2 djþ1=2 b ¼ m
ð25Þ
This scheme admits a maximum principle on b and h under a CFL condition. Mass, momentum and energy are also conserved. 4.4. Space and time accuracy We consider a predictor corrector scheme, and assume for simplicity that the predictor and corrector steps are or the form þ þ ; Wjþ1=2 Þ F ðWj1=2 ; Wj1=2 Þ ð26Þ Wjnþa ¼ Wjn am F ðWjþ1=2 where a ¼ 12 or 1, Wj1=2 are evaluated via the MUSCL method, but the same type of technique can also be extended to the Harten–YeeÕs scheme [27]. The trick is to rewrite the scheme in fluctuation form, Wjnþa ¼ Wjn am wajþ1=2 þ waj1=2
where þ ; Wjþ1=2 Þ F ðWj Þ :¼ wjþ1=2 þ wjþ1=4 wajþ1=2 ¼ F ðWjþ1=2 þ waj1=2 ¼ F ðWj Þ F ðWj1=2 ; Wj1=2 Þ :¼ wj1=2 þ wj1=4 ; þ ; Wjþ1=2 Þ F ðWjþ1=2 Þ, wjþ1=4 ¼ F ðWjþ1=2 Þ F ðWj Þ, so that we can apply the where wjþ1=2 ¼ F ðWjþ1=2 same machinery as before. If the base scheme is RoeÕs, GodunovÕs or RusanovÕs, we see that the
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
583
predictor and corrector are the one given by the formula (23), (24) or (25) evaluated at the corresponding states, and modified by the following corrective term djþ1=2 b ¼ m u ðujþ1=2 ÞDjþ1=2 b þ u ðujþ1=4 ÞDjþ1=4 b djþ1=2 h ¼ m u ðujþ1=2 ÞDjþ1=2 h þ u ðujþ1=4 ÞDjþ1=4 h ð27Þ dj1=2 b ¼ m uþ ðuj1=2 ÞDj1=2 b þ uþ ðujþ1=4 ÞDj1=4 b dj1=2 h ¼ m uþ ðuj1=2 ÞDj1=2 h þ uþ ðuj1=4 ÞDj1=4 h where the functions u ð:Þ and uþ ð:Þ depends of the solver used in the predictor (resp. corrector) step. We have used the notations: þ
Djþ1=2 ð:Þ ¼ ð:Þjþ1=2 ð:Þjþ1=2 ; Dj1=2 ð:Þ ¼ ð:Þþ j1=2 ð:Þj1=2 ;
Djþ1=4 ð:Þ ¼ ð:Þ jþ1=2 ð:Þj ; Dj1=4 ð:Þ ¼ ð:Þj ð:Þþ j1=2 ð:Þj1=2 , ð:Þj1=2 the average values between ð:Þþ j1=2 and ð:Þj1=2 , ð:Þjþ1=4 and ð:Þjþ1=4 the average þ values between ð:Þ jþ1=2 and ð:Þj , ð:Þj1=4 and ð:Þj1=4 the average values between ð:Þj and ð:Þj1=2 .
4.5. An illustration for multi-dimensional finite volume approximation In this section, we develop an extension of the previous scheme to multi-dimensional problems. We start from a standard conservative numerical scheme for (2). The scheme is Dt X / ð28Þ Wi nþ1 ¼ Wi n ai j2#ðiÞ ij with the conserved variables are W ¼ ðq; m; EÞT , Dt is the time step, ai is the area of the cell (see Fig. 3), #ðiÞ is the set of the neighboring cells of the cell i and /ij is the fluctuation associated to the cell interface ðijÞ. This fluctuation is computed from the solution of a Riemann problem at the cell interface. We illustrate the approach in the case of the Roe scheme. We use the notations of Appendix B where construction of Roe scheme in 2D is sketched. As pointed out above, the method can easily be extended to Godunov and RusanovÕ schemes, following the same ideas as in Sections 4.1, 4.2 and 4.3. We start with some notations: T uu V2 ; Q? ¼ ; W ¼ ðq; qU; EÞT ; F ¼ qU; qU2 þ p; qH U ; 2 2 The system (4), in the local coordinates is given by: 8 ot W þ og F ¼ 0; > > < o ðqVÞ þ o ðqUVÞ ¼ 0; Q¼
t
g
ot ðqQ? Þ þ og ðqUQ? Þ ¼ 0; > > 2 : bp þ h ¼ E qU2 qQ?
Z ¼ ðW; qVÞT
ð29Þ
584
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Fig. 3. Dual cell Ci associated with triangulation vertex Mi in 2D space.
We also set 0
W
1
C B B qV C C Z¼B B qb C; A @
0
1
C B B qVU C C G¼B B qbU C; A @ qhU
qh W ¼ M1 ij Z
F
0
1
0
B0 g B x Mij ¼ B B0 0 @ 0 g? x
0 gy 0 g? y
0
1
0C C C; 1C A 0
Z ¼ Mij W ;
and J ¼ oZ G
The numerical dissipation of Roe scheme is
e ÞDZ ¼ uðUÞðDðqÞ nÞr1 þ uðU c~Þ n H r2 þ uðU þ c~Þ n þ H r3 uð J 2 2 Db bDh hDb þ uðUÞ qDðVÞ r4 þ uðUÞq r5 þ uðUÞ r6 b qb
ð30Þ
with, see Appendix B, Dp DU DU and H ¼ q e c~2 c The fluctuation Wij at the cell interface is given by: 0 Z1 Wij n o 1 B qb C ~ ~ Wij ¼ DG uðJÞDZ ¼ u ðJÞDZ ¼ @ Wij A 2 Wqh ij n¼
ð31Þ
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
585
Fig. 4. Area (ai ) of the cell and partial area (aij ) associated to a given segment ði; jÞ.
Therefore, the update of the conservatives variables W on the internal cells, is Wi nþ1 ¼ Wi n
X aij Dt X 1 Z nþ1 Mij Wij ¼ M1 ij Zij ai j2#ðiÞ a j2#ðiÞ i
where Znþ1 ¼ Znij ðDt=aij ÞWZ ij . Here, aij is the component of the cell area associated to the ij interface (see Fig. 4). We denote by Wij the first three components of WZ ij and obtain: Wnþ1 ¼ Wnij ij
Dt Wij aij
ð32Þ
This relation provides the updated variables associated to the sub-cell aij of the cell ai . Therefore, the projection on the cell ai of the variables is obtained by averaging the sub-cells variables weighted by the ratio aij =ai . This is one of the keys of the approach proposed in this paper. Indeed, with this decomposition, we can focus the study on the sub-cells aij : the situation is now very similar to the 1D case. Let us denote dXij ¼ Xijnþ1 Xijn and K ¼ qQ ¼ 12qU2 þ 12qV2 ¼ K þ K ? Here, the local frame is defined by the normal to aij at the interface between the cells ai and aj , and the orthogonal vector, such that the frame is direct. We have 1 nþ1 nþ1 qij Uij dUij þ Unij dðqUÞij dEij ¼ dðqÞij þ dKij? þ dKij ; dKij ¼ 2 n nþ1 and dðqUÞij ¼ qij dUij þ Uij dqij . For simplicity, we denote: Uni Unþ1 Unþ1 þ Uni ij ij T ¼ ; Uij ¼ and ‘ij ¼ ðQk ij ; Uij ; 1Þ 2 2 k ? Therefore, dðqÞij þ dKij ¼ Qij dqij Uij dðqUÞij þ dEij so that denoting rij ¼ Dt=aij Qk ij
n ? bnþ1 ij dpij þ pi dbij þ dhij þ dKij ¼ rij ‘ij Wij
ð33Þ
586
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Setting u ðX Þ ¼ ðX uðX ÞÞ=2, the components of the numerical fluctuation are given by: Wqij ¼ u ðUÞðDq nÞ þ c~ a 1 n þ a2 H p q WqU ij ¼ UWij þ Wij
WEij ¼ HWqij þ N where Wpij ¼ c~2 a 2 n þ a1 H , c2 ðDq nÞ þ qVDV þ pDb þ Dh ; N ¼ UWpij þ u ðUÞ b~ a 1 ¼
u ðU þ c~Þ þ u ðU c~Þ ; 2~ c
a 2 ¼
u ðU þ c~Þ u ðU c~Þ 2~ c
Using these expressions of the fluctuations associated to the edge ði; jÞ, we compute: 2 ~ ? Wq þ u ðUÞqVDV ‘ij Wij ¼ Sij þ b c~a Dp þ a q~ c DU þ u ðUÞðpi Db þ DhÞ þ Q 1 2 ij
ð34Þ
where k ~ k Wq þ Wpij U U þ u ðUÞðp pi Þ Sij ¼ Qij Uij U þ Q ij ij In general, Sij is consistent with zero and equal to zero if the velocity and the pressure field is a ~ c2 DU is consistent with space–time constant. rij b~c~a 1 Dp is consistent with bUog p and rij ba2 q~ 2 bqc og U. Therefore a natural discretization of the EOS parameters is obtained by a combination of Eqs. (33) and (34) and an identification of the different terms with respect to the non-conservative formulation in the direction gij : dbij ¼ rij u ðUÞDij b
and dhij ¼ rij u ðUÞDij h
ð35Þ
If we denote dbji and dhji the variations on the cell j associated to the interface, it is easy to see that: dbji ¼ rji uþ ðUÞDij b and dhij ¼ rji uþ ðUÞDij h where rji ¼ Dt=aji . These relations are obtained under the assumption that the update of the tangential kinetic e ? Wq þ u ðUÞqVDV : We set energy follows: dKij? ’ rij Q ij ? e ? Wq þ u ðUÞqVDV WKij ¼ Q ij
ð36Þ
Because the stability of the contact discontinuity is marginal in 2D and 3D [24], the previous is generally violated when the tangential velocity is not continuous, even for perfect gas EOS. This is why we independently discretize the tangential kinetic energy, namely ðKij? Þnþ1 ¼ ðKij? Þn rij WKij
?
In order to take into account this auxiliary relation in the numerical approach, we define a local update of the internal energy q by: nþ1 ðqÞij
¼
Eijnþ1
nþ1 2 ðqnþ1 ij Uij Þ
2qnþ1 ij
ðKij? Þ
nþ1
¼ Ein Kin rij WEij þ rij WKij
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
587
where qU ðUni Þ2 q rij Wij qU n q ¼ þ þ W W Wij þ 2U i ij ij 2 2qni From this discretization, we can construct an average internal energy on the cell that is locally consistent. Of course, this energy is different to the one that can be computed from the conservative variables on the cell at the time tnþ1 . Therefore we propose to update the total energy by the b nþ1 satisfying the relation: value E i nþ1 nþ1 X aij b nþ1 mi mi ðqÞnþ1 ¼ E ij i a 2qnþ1 i j2#ðiÞ i
WKij
? WKij
Uni WqU ij
This means that the modified value of the energy at the time tnþ1 is given by a non-conservative formulation: mnþ1 mnþ1 mn mn X E K i i b nþ1 ¼ En Dt E W i n i ð37Þ W ij ij þ i i ai j2#ðiÞ 2qi 2qnþ1 i The update of the material parameters in the cell is defined as follows: Dt X Dt X n nþ1 n bnþ1 ¼ b u ðUÞD b and h ¼ h u ðUÞDij h ij i i i i ai j2#ðiÞ ai j2#ðiÞ
ð38Þ
Summary: The numerical scheme is then defined by the classical conservative formulation for the mass and momentum. The total energy update is obtained by a non-conservative formulation (37). The material parameters are given by (38). The resulting non-conservative scheme is a more general formulation of the idea developed in [21] with a splitting approach. A similar procedure is used in [3]: the important variables are the primitive variables and not the conservative ones. To illustrate this point, we present a quasi-1D test case where this default affects dramatically the solution of the Euler equation in 2D (see Fig. 8). 5. The reactive model We now detail the extension of this method for materials that the physical state changes during chemical reaction: solid energetic materials. This can easily be done because the formulation is the same everywhere. Solid explosives start from the solid state before reaction, then become a solid–gas mixture during reaction and a pure gas mixture at the end of the reaction. We first give some basis regarding the thermodynamics of such materials. We then give the some elements regarding their chemistry. We finally detail the way to account of the equation of state changes during the decomposition process and the way to couple it with the numerical treatment of the interfaces. This is important for all applications involving interactions of detonations with neighboring materials. 5.1. Flow model When dealing with detonation waves in condensed explosives the pressure and velocity conditions are so high that liquids and solids become compressible and governed by the gas dynamics
588
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
equations. Thus, system (2) becomes valid even for solids. The closure of such system is achieved by Eq. (1) that needs some modifications to account for changes in chemical composition of the mixture. This will be detailed in the next subsection. However, let us mention that the system (2) is not sufficient when solving detonation waves with an Eulerian method as detailed previously. It must be complemented by an equation allowing the determination of the phases concentration, varying inside the reaction zone of the detonation wave. We define the gas volume fraction k, whose initial value is zero in the solid explosive and that tends to the unity at the end of the reaction. The explosive being moving, the volume fraction equation reads: ð39Þ ot qk þ r ðqukÞ ¼ qk_f where k_ represents the rate of creation of the gas phase. Assuming that the rate of decomposition obeys a simple pressure kinetics (Vieille law), we get: k_ ¼ ð1 kÞapn
ð40Þ
where a and n are the kinetics parameters of the explosive. The function f is an indicator that allows the determination of the explosive concentration in the mesh. Indeed, the decomposition law must be applied only at points where the explosive is present and not in inert materials. Spatial localization of the explosive is realized by the introduction of function f that equals zero initially everywhere, except in the explosive where it equals 1. Its evolution is given by: ot qf þ r ðquf Þ ¼ 0
ð41Þ
5.2. Mixture equation of state During the reaction process the explosive passes from solid state to gaseous products. The equation of state (1) must be reformulated to deal with this transformation. During the detonation propagation, there are three states in the flow: state 1 where the solid explosive is intact, state 2 where the explosive is completely reacted and composed exclusively of gases and an intermediate state composed of a mixture of solid and gases inside the reaction zone. To deal with this various states we propose a mixture equation of state able to reproduce all of them in a unique formulation. Each pure fluid is assumed to obey the stiffened gas EOS. Thus the internal energy of the solid reads: s ¼
bs þ hs þ cs qs
with
cs ¼ s0
hs qs0
where s0 , cs and qs0 are respectively the formation energy of the solid, the internal energy and density under standard conditions. For the gas phase, we assume the same formulation but instead of choosing standard conditions as reference, we retain the Chapman-Jouguet (CJ) point. Thus the internal energy reads: g ¼
bg þ hg þ cg qg
with cg ¼ cj
hcj qcj
where hcj ¼ ðpcj þ cg p1;g Þ=ðcg 1Þ and the characteristics states qcj , cj and pcj can be obtained from a thermochemical code like CHEETAH [8].
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
589
With the internal energies definition, the mixture equation of state is readily obtained under the assumption of pressure equilibrium between phases. This assumption is justified because of the small size of the solid particles into the mixture and of the high speed of sound of both two phases. Indeed, the pressure relaxation time is closed to the time required for an acoustic wave to cross a solid grain diameter. The mixture internal energy reads: q ¼ kqg g þ ð1 kÞqs s
ð42Þ
and now is a function of q, p and k. According to the previous definitions for the internal energies (s and g ) and assuming that gas and solid are at the same pressure p ¼ pg ¼ ps , we recover an EOS of the same type as the pure phases. It reads: ¼
bþh þc q
with
cðkÞ ¼ kcg þ ð1 kÞcs ; bðkÞ ¼ kbg þ ð1 kÞbs ; hðkÞ ¼ khg þ ð1 kÞhs
ð43Þ
It is now necessary to determine the decomposition parameters of the kinetics law of the explosive (parameters of the Vieille law). 5.3. Chemical kinetic The parameters of the VieilleÕs law are not the standard ones as determined in stand burner or closed bomb. Indeed, the parameters determined by such a way are those of the deflagration regime of the explosive. At the pressure levels we are interested in, these parameters cannot be determined by direct measurement. Also, the kinetic parameters we are looking for are not related to the ignition and transition regime of the explosive. Such mechanisms cannot be summarized in a simple Vieille law [16]. The kinetic parameters needed here are related only to the self-detonation propagation regime of the explosive. To determine these parameters we use an experimental observation based on the front curvature of the detonation wave. The energy responsible for the front propagation is liberated between the shock front and the sonic plane (in the frame of reference related to the shock front). The distance between these two planes being finite, multi-dimensional rare-fraction waves coming from the boundaries of the explosive product are able to propagate inside the reaction zone of the explosive. These effects result in a curvature of the shock front. Thus, it is clear that the shock curvature is a result of the coupling effects between the chemical kinetic of the explosive and the various hydrodynamics effects, especially the rare-fraction waves. Considering a cylindrical charge of given diameter, the corresponding detonation wave propagates at a given velocity with a given curvature. By measuring these variables for various charge diameters we obtain a relation between the front curvature and its velocity. This relation is intrinsic to the explosive. It is a result of its chemical kinetic and rare-fraction waves. This relation will enable us to determine the kinetic parameters in a wide range of pressures. Indeed, when the diameter varies, the detonation velocity and the pressure inside the reaction zone vary strongly. By using the equation of state for the mixture, combined with system (2) and Eq. (39) in the context of a steady 2D axi-symmetric flow in the frame of the shock and specialized onto the symmetry axis, we get an ODE system easy to solve. It corresponds to the Wood and Kirkwood
590
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
model [7]. Giving an arbitrary guess for the chemical kinetic parameter of the VieilleÕs law and the front curvature for a given charge diameter this model furnishes the corresponding stable detonation velocity. This detonation velocity accounts for the coupling between the hydrodynamics and the chemical kinetics. The kinetic parameter may be varied in order to reach the experimental detonation velocity. This type of process may be done for all front curvatures, that is for all charge diameters. When the computed detonation velocities and the experimental ones are in good agreement for any front curvature, the kinetic parameters are determined. This type of inverse procedure is detailed in [17] and [20]. 5.4. Numerical approach for the reactive model In the reactive case, the equation of state (43) is now a function of three variables. The mixture internal energy can be written as: ¼ e þ c. Thus the total energy equation ot E þ r ðqHuÞ ¼ 0 can be put under the form: ot E þ r ðqH uÞ ¼ qk_c0 0 with E ¼ qe þ qðu uÞ=2 and for any variable ðÞ ¼ dðÞ=dk. The contribution qk_c0 is related to the combustion energy release. The equations governing the material interfaces motion must also be modified. Indeed, in the reactive case, the material parameters b and h of Eq. (43) now depend on the volume fraction k. Therefore, the materials parameters are governed by the equations:
ot b þ u rb ¼ f k_b0
and ot h þ u rh ¼ f k_h0
ð44Þ
To summarize, in the reactive case, the overall system to solve reads: ot q þ rðquÞ ¼ 0 ot qu þ rðqu uÞ þ rp ¼ 0 ot E þ rðqH uÞ ¼ qk_c0 ot qf þ rðqf uÞ ¼ 0
ð45Þ
ot qk þ rðqkuÞ ¼ qk_f ot qb þ rðqbuÞ ¼ qf k_b0 ot qh þ rðqhuÞ ¼ qf k_h0 The overall system is split into a convective one and an ODE system. The convective part of this system is solved by the hyperbolic solver detailed in the preceding sections. The source terms are then accounted for by an appropriate ODE solver.
6. Numerical experiments The technique developed in this paper is now applied to three tests cases: a water–air shock tube, an interaction of a shock with an air bubble in water and a detonation wave. The com-
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
591
putations are performed with 2D unstructured meshes for the two first cases and with structured meshes for the last one. In any cases, the methodology developed in the framework of unstructured grids is applied and a second order accurate MUSCL explicit scheme is used. 6.1. Water–air shock tube We consider a square domain of 1 m 1 m length. Initially there is a material interface, located at x ¼ 0:7 m, separating two fluids: water at the left and air at the right. The two fluids are initially at rest and there is a jump on the density, the pressure and the material parameters (b and h). The solution of this Riemann problem is composed by a shock wave travelling in the air, a rarefraction travelling in the water and the material interface (contact discontinuity) [9]. The explicit second order accurate scheme used here is based on the classical RoeÕs solver with CFL ¼ 0:75. The computations of the first test case (Fig. 5) are performed on an unstructured mesh of 4506 vertices with about 750 vertices in the wave propagating direction. The numerical results are in good agreement with the analytical solution. For the sake of clarity, only the numerical solutions, at the final time T ¼ 230 ls, are plotted (Fig. 6). The 1D structure of the 1D analytical solution is preserved by this modified conservative approximation. Moreover, there are no oscillations at the material interface (Fig. 6). The solution with the coarse mesh (1782 vertices) is less accurate but the global profile of the waves are captured (Fig. 7). The solution of the second test case (Fig. 8) is, analytically, composed by the solution of the test case one and a strong contact discontinuity on the v component of the velocity. When a conservative scheme is used (Roe, Rusanov or HLL solvers), the numerical solution consist in a shock wave travelling faster than the exact one. This is a consequence of the lack of accuracy, on the tangential velocity (v), for the conservative schemes. Indeed, the pressure and the sound speed, used to numerically compute the waves speed, are directly related to the norm of the velocity. Computations performed on structured or unstructured meshes with first or
Fig. 5. Initialization of the water–air shock tube: Test 1.
592
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Fig. 6. Density, velocity pressure and EOS parameters at the final time 230 ls, RoeÕs flux and second order accurate MUSCL scheme on unstructured mesh.
second order accurate conservative schemes give rise to the same problem (Fig. 9). We have also computed some numerical solutions with the partially non-conservative scheme. Therefore, the numerical results, computed with this non-conservative scheme (Fig. 9), are very close to the analytical one.
6.2. Shock–bubble interaction In this section, we investigate the dynamics of a shock–bubble interaction. The size of the computational domain is 20 10 cm2 , with absorption boundary conditions on the right and piston conditions on the left. The initial flow is water at rest with an cylindrical air bubble, of 4 cm radius, centered at the point (5, 5) cm. The system is at the equilibrium with uniform pressure and velocity then a shock wave is produced at the left of the domain by the piston motion (Fig. 10). In this test case, an unstructured mesh of 98 489 nodes and 195 776 triangles is used. Numerical results are computed with an explicit second order accurate scheme based on the Roe solver and the CFL ¼ 0:5. In this computation, the material interface is not explicitly defined, hence we represent it by the isolines b ¼ 1:5. The speed of sound inside the air bubble is smaller than that of the surrounding water. After the first interaction with the bubble, the incident shock wave is diffracted into two waves (Fig. 11): a
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
593
Fig. 7. Comparison between the computations with the second order accurate MUSCL scheme for an unstructured coarse mesh of about 100 points in the propagation direction and an unstructured fine mesh of about 750 points in the propagation direction.
reflected rare-fraction wave moving from right to left in the water and a transmitted shock wave in the bubble. Later, when the incident wave have passed the bubble, we can observe the multiple scattering phenomena (Figs. 12 and 13). Due to the fact that air compressibility is larger than the water one, the downstream bubble interface is almost motionless during the compression (Fig. 16). We notice that the numerical scheme preserves this material interface and no spurious oscillations or non-physical states are produced during the computation. Then, the transmitted shock emerges in the water (Fig. 14 for t ’ 18:5 ms) and deforms the downstream material front. Two pairs of vortices (Figs. 15 and 16) are produced around the compressed air bubble. The upstream interface is folded and progressively a re-entrant jet is formed. This jet breaks the bubble in two parts and drives them downstream (Fig. 17). Even if the context is not the same, the behavior of the shock–bubble obtained in this investigation is similar to the experiments realized in [13].
594
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Fig. 8. Initialization of the water–air shock tube, Test case 2 (left) and numerical solution with the RusanovÕs flux, first order conservative scheme at the final time 230 ls (right): structured mesh of 705 3 nodes.
6.3. Detonation wave We consider a cylindrical Plexiglass tube filled with a solid explosive. The tube is surrounded by air. The explosive is a mixture of ammonium nitrate, fuel, water and other minor compounds. It is a non-ideal explosive: its chemical decomposition rate is very slow compared to standard solid high explosives. Consequently the consideration of the chemical kinetic effects are very important for the detonation propagation and its effects over structures. The configuration studied consists in several Plexiglass tubes of different diameter. The tube thickness is 3 mm. The initial density of Plexiglass is 1186 kg/m3 and the parameters of the equation of state are c ¼ 3:8 and p1 ¼ 2 109 Pa. The EOS parameters for the solid state of the explosive are: s0 ¼ 0, cs ¼ 2:9, p1s ¼ 4:5 108 Pa and the initial density is q0 ¼ 1000 kg/m3 . The EOS parameters of the detonation products are: cg ¼ 2:37, p1g ¼ 0, pCJ ¼ 11 109 Pa, DCJ ¼ 5875 m/s and qCJ ¼ 1477 kg/m3 . The CJ state is obtained by a thermochemical equilibrium computation using CHEETAH [8]. Finally, the kinetic parameter of Eq. (40) are a ¼ 3:34 103 Pan and n ¼ 0:8 when p > p0 and a ¼ 0 otherwise. The threshold pressure p0 is take equal to 1 GPa. The following computation is performed with a mesh of 300 75 cells and we first consider a cylindrical charge of 84 mm diameter (Fig. 18). The computational domain is 30 cm long and 6 cm high. The results are shown when the stable detonation state is obtained (Fig. 18). The detonation front is clearly visible as well as the expansion of the Plexiglass cylinder. These results are to compare with X-rays measurements as shown in Fig. 19. The qualitative agreement between the computations and experiment can be noticed. The quantitative comparison will be given after examining the effects of the charge size. We now consider another charge of 174 mm internal diameter. The mesh involves now 400 150 cells. The computed density contours are shown in Fig. 20. Again, the shock front and tube expansion are clearly visible. The numerical data are now compared with the experimental ones in Table 1. From these results an excellent agreement can be noticed. It is important to note the wide range of variation of the detonation velocity compared to the ideal one (5875 m/s). This is a result
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
595
Fig. 9. Comparison between the conservative and the non-conservative scheme. Numerical solution with the RusanovÕs flux, first order conservative scheme, at the final time 230 ls: structured mesh of 705 3 nodes.
Fig. 10. Initialization of the shock–bubble interaction test case.
596
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Fig. 11. Evolution of the reflected, the refracted and the incident shock waves. Times t ¼ 0:65, 1.30, 1.95 and 2.6 ms.
Fig. 12. Pressure lines at t ¼ 9:75 ms: refracted and diffracted shock waves.
of the coupling effects between the low rate of the chemical kinetic of this explosive, and the rarefraction waves travelling inside the reaction zone.
7. Conclusion In the case of the stiffened gas EOS, we have pointed out the different problems arising in the computation of the multi-dimensional multi-material flows: the pressure can become oscillatory across slip lines, and consequently the velocity and the other variables too. For very different EOS, and/or extreme flow configurations this leads to the breakdown of the code. The situation is even worse in several space dimensions because the kinetic velocity is not accurately computed around slip lines.
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
597
Fig. 13. Streamlines at t ¼ 9:75 ms: refracted and diffracted shock waves.
Fig. 14. Evolution of the pressure (p þ p1 ) lines. Times t ¼ 0, 6.5, 13.0, 18.5, 26.0 and 32.5 ms.
We have discussed these different problems and propose solutions. They amount to a slight modification of existing schemes; they are easily incorporated into existing codes. These
598
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Fig. 15. Evolution of the material interface (b ¼ 1:5). Times t ¼ 0, 6.5, 13.0, 18.5, 26.0 and 32.5 ms.
Fig. 16. Evolution of the material interface (b ¼ 1:5). Times t ¼ 0, 6.5, 13.0, 18.5, 26.0 and 32.5 ms.
modifications are of two types: first the identification of the numerically pertinent thermodynamical parameters that have to be updated, second more consistent calculation of the kinetic energy. We note that the second ingredient of the method arises because of the multi-dimensional nature of the problems and the class of conservative schemes (directional splitting or finite volume
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
599
Fig. 17. Vortex pairs developed around the bubble at t ¼ 26:0 ms.
Fig. 18. Detonation wave in a 84 mm diameter cylindrical charge in a Plexiglass tube. Density contours.
approach) we consider here. The drawback of our method is that the scheme is not any longer strictly conservative on the total energy. We note however that the loss of conservation occur only around slip lines and does not seem to pollute the solutions. The method is flexible, it is applicable on any type of mesh: Cartesian, unstructured, hybrid, because of a local formulation.
600
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Fig. 19. X-ray measurement of the shock front and cylinder expansion from Saurel et al. [23].
Fig. 20. Detonation wave in a 174 mm diameter cylindrical charge in a Plexiglass tube. Density contours.
Table 1 Experimental vs computational results for the propagation of detonation waves in explosive cylinders of various charge diameter. Experiments: left column. Computations: right column Charge diameter (mm)
Detonation velocity
Radius of curvature
Expansion angle
Exp. (m/s)
Comp. (m/s)
Exp. (mm)
Comp. (mm)
Exp.
Comp.
84 104 174
3060 3467 4000
3160 3369 4006
96 116 290
104 120 285
32 32 33
31 32 32
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
601
Several numerical examples, with and without the corrections, are presented. They demonstrate the efficiency of our method on triangular unstructured meshes. Acknowledgements This paper has been completed during the CEMRACS 2000 summer school, held at CIRM (Marseilles) during the summer 2000. The third author has also been partially supported by NITROCHIMIE. He addresses special thanks to S. Mencacci and T. Rousse. The experiments shown for the detonation problem have been done at CEA/DAM Moronvilliers and we particularly thank G. Rodriguez. This research has been conducted as a research project of the Laboratoire Correspondant du CEA LRC M03. Appendix A. Derivation of the Roe average for the stiffness gas EOS In the case of the stiffened gas, the equation of state is bp þ h ¼ q ¼ E 12qu2
ðA:1Þ
In the 1D case, the fluid is governed by: ot W þ ox F ¼ 0
ðA:2Þ
with initial and boundary conditions. The vector of conserved variables is W ¼ ðq; qu; E; qb; qhÞ :¼ ðq; m; E; B; HÞ and the flux is F ðW Þ ¼ ðqu; qu2 þ p; uðE þ pÞ; uqb; uqhÞ. The Jacobian of F with respect to W is 1 0 0 1 0 0 0 C B C B u2 þ pq 2u þ p p p p m E B H C B B uðH þ pÞ H þ upm uð1 þ pE Þ upB upH C ðA:3Þ C B C B ub b 0 u 0 A @ uh h 0 0 u where þ pq ¼
b
u2 2
pm ¼
u b
pE ¼
1 b
pB ¼
p qb
pH ¼
1 qb
ðA:4Þ
If we set Df ¼ fL fR , the Roe average is defined by replacing in (A.3) the Roe average of the velocity and the total enthalpy H and by defining the average derivative of the pressure such that the following jump relation is satisfied, e f e f 1 p 1 Dq þ DðqÞ DðqbÞ DðqhÞ ¼ Dp b b qb qb
ðA:5Þ
Let k 2 ½0; 1 and define for any function f, f ¼ kfL þ ð1 kÞfR and f ¼ kfR þ ð1 kÞfL . For any functions f and g, we have
602
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
Dfg ¼ f Dg þ gDf pffiffiffiffiffi pffiffiffiffiffi pffiffiffiffiffiffi Last, if k ¼ qL =ð qL þ qR Þ, we have qf ¼ qf . Using the above calculation rules, we can rewrite Dp as q h 1 1 1 Db Dh Dp ¼ qD þ Dq 2 b b b b where we have set b2 ¼ bL bR . Last, by proceeding by identification, we get ! f b f e 1 e 1 1 p 1 q h 1 1 :¼ :¼ ¼ 2 :¼ :¼ 2 b b b b b qb qb qb q b The matrix of eigenvectors is easily computed, 0 1 1 B u B 2 u 2 =2Þ ðu =2Þ ðu R ¼ ð r1 r2 r3 r4 r5 Þ ¼ B B @ b 0 h 0
1 1 1 1 u uc uþc C C 2 ðu =2Þ H uc H þ uc C C 0 b b A q h h associated to the eigenvalues u, u, u, u c, u þ c where c2 ¼ b H ðu2 =2Þ . Similarly, any vector W ¼ ðA; B; C; D; EÞ in R5 can be decomposed as: W ¼
5 X
ðA:6Þ
ðA:7Þ
li ðW Þri
i¼1
We get, for W ¼ ðDq; Dm; DE; DB; DHÞ, 1 bDh hDb Dp l2 ðW Þ ¼ Dq 2 l1 l3 c h q Dh Db l3 ðW Þ ¼ b q 1 Dp l4 ðW Þ ¼ Du 2 c2 c 1 Dp q l5 ðW Þ ¼ þ Du 2 c2 c l1 ðW Þ ¼
Appendix B. 2D Roe type scheme This section is devoted to the derivation of Eq. (30) in Section 4.5. We first introduce some notations. We set
ðA:8Þ
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
1 W B qV C C B Z¼B C; @ qb A 0
1 F B qVU C C B G¼B C; @ qbU A 0
qhU
qh and
0
1 B B0 B Mij ¼ B B0 @ 0
0 gx 0 g? x
0 gy 0 g? y
1 0 C 0C C C; 1C A 0
Z ¼ Mij W ;
603
W ¼ M1 ij Z
1 S0 B US0 þ S1 C C B B H S0 þ US1 C C Jðu; H; b; h; S0 ; S1 Þ ¼ oZ G ¼ UId þ B C B VS0 C B A @ bS0 hS0 0
where S0 ¼ ðU
1 0
0
S1 ¼ ð oq p
oqU p
oE p
0
0 Þ ¼ qrZ U oqV p
oqb p
oqh p Þ ¼ rZ p
We denote by Zi and Zj the values of the conservative variables in the cells i and j. For any function f, we set f ¼ kfi þ ð1 kÞfj and f ¼ kfj þ ð1 kÞfi . These two averages have the property that DðfgÞ ¼ f Dg þ gDf Moreover, if pffiffiffiffi qi k ¼ pffiffiffiffi pffiffiffiffiffi qi þ qj simplify the following calculations. we have qg ¼ qg: These remarks e u; H ; b; h; S ; S with q~ ¼ bp þ h, We define the matrix J 0 1 S0 ¼ U 1 0 0 0 0 p 1 e U 1 V ~ þ Q S1 ¼ q b
1 q
e satisfy the jump condition DG ¼ J e DZ (DZ ¼ Dij Z ¼ e ¼ ðU2 þ V2 Þ=2. The matrix J where Q Zi Zj ). Moreover, we have the relations: S0 DZ ¼ qDU;
S1 DZ ¼ Dp
604
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
e is diagonalisable and the eigenvectors are The matrix J 1 1 0 1 0 0 1 1 1 BUC B U c~ C B U þ c~ C C C B C B B C C BeC B B BQC B H c~U C B H þ c~U C C C C B B r1 ¼ B B V C; r2 ¼ B V C; r3 ¼ B V C; C C B C B B C C B C B B A A @bA @ @ b b h 0 1 0 B0C B C B C B ~ C C r5 ¼ B B 0 C; B C B C @bA
0 1 0 B0C B C B C B1C C r6 ¼ B B0C B C B C @0A
h
0
0
1
B 0 C B C B C BVC C r4 ¼ B B 1 C; B C B C @ 0 A 0
h
q
h
e . The associated eigenvalues are where b~ c2 ¼ H Q k1 ¼ U;
k2 ¼ U c~;
k3 ¼ U þ c~;
k4 ¼ U;
k5 ¼ U;
k6 ¼ U
In the eigenvectors basis, the jump of conservative variables is given by: DZ ¼ ðDðqÞ nÞr1 þ þ
nH nþH DðqbÞ bDq r5 r2 þ r3 þ DðqVÞ VDðqÞ r4 þ 2 2 b
bDðqhÞ hDðqbÞ r6 qb
where 1 DðqbÞ bDq bDðqhÞ hDðqbÞ VDðqVÞ UDðqUÞ þ QDðqÞ þ DðEÞ ~ n¼ 2 b~ c b qb
ðA:9Þ !
Dp DðqUÞ UDðqÞ and H ¼ 2 c~ c~ After some tedious but straightforward calculations, we get ¼
n¼
Dp c~2
and H ¼ q
DU c~
References [1] Abgral R. Generalisation of the Roe scheme for the computation of a mixture of perfect gases. La Recherche Aerospatiale 1988;6:31–43. [2] Abgrall R. How to prevent pressure oscillations in multicomponent flows: a quasi conservative approach. J Comput Phys 1996;125:150–60.
R. Abgrall et al. / Computers & Fluids 32 (2003) 571–605
605
[3] Abgrall R, Kami S. Compressible multifluids. J Comput Phys 2001;169:594–623. [4] Chargy D, Abgrall R, Fezoui L, Larrouturou B. Conservative numerical schemes for multicomponent inviscid flows. La Recherche Aerospatiale 1992;2:61–79. [5] Chern I-L, Glimm J, McBryan O, Plohr B, Yaniv S. Front tracking for gas dynamics. J Comput Phys 1986;62:83– 110. [6] Fedkiw RK, Aslam T, Merriman B, Osher S. A non oscillatory eulerian approach to interfaces in multimaterials flows (the ghost fluid method). J Comput Phys 1999;152:457–92. [7] Fickett W, Davis WC. Detonation. Technical report, University of California Press, Berkeley, 1979. [8] Fried E. CHEETAH 1.39 userÕs manual. Lawrence Livermore National Laboratory, UCRL-MA-117541 Rev 3, Energetic Material Center, 1996. [9] Godlewski E, Raviar PA. Numerical approximation of hyperbolic systems of conservation laws. In: Applied Mathematical Sciences, vol. 118. New York: Springer; 1995. [10] Godunov SK. A finite difference method for the numerical computation of discontinuous solutions of the equations of fluid dynamics. Mat Sb 1959;47:357–93. [11] Godunov SK, Zabrodine A, Ivanov M, Kraiko A, Prokopov G. Resolution numerique des probl emes multidimensionnels de la dynamique des gaz. Moscow: Editions Mir; 1979. [12] Hartlow F, Amsden A. Fluid dynamics. Technical Report LA-4700, Lawrence Livermore National Laboratory, 1971. [13] Hass J-F, Sturtevant B. Interaction of a weak shock waves with cylindrical and spherical gas inhomogeneities. J Fluid Mech 1987;181:41–76. [14] Karni S. Hybrid multifluid algorithms. SIAM J Sci Comput 1996;17:1019–39. [15] Karni S. Multi-component flow calculations by a consistent primitive algorithm. J Comput Phys 1996;112:31–43. [16] Massoni J, Saurel R, Baudin G, Demol G. A mechanistic model for shock initiation of solid explosives. Phys Fluids 1999;11:710–36. [17] Persson PA, Holmberg R, Lee J. Rock blasting and explosive engineering. Boca Raton: CRC Press; 1996. [18] Roe PL. Riemann solvers, parameter vectors and difference schemes. J Comput Phys 1981;43:357–72. [19] Rusanov VV. Calculation of interactions of non steady shock waves with obstacle. J Comput Math Phys 1961;1:267–79. [20] Saurel R. DPEX: Un outil numerique pour la determination des cinetiques reactionnelles des explosifs fortement non ideaux. Rapport de contrat, Centre dÕEtudes de Gramat, IUSTI, Universite de Provence, Marseilles, France, 1996. [21] Saurel R, Abgrall R. A multiphase Godunov method for compressible multifluid and multiphase flows. J Comput Phys 1999;150:425–67. [22] Saurel R, Abgrall R. A simple method for compressible multifluid flows. Siam J Sci Comp 1999;21:1115–45. [23] Saurel R, Mencacci S, Mercier S. Hydro3d: Un outil simple pour la simulation de probl emes de detonique. In: EuroPyro 99. 7 eme congr es international de pyrotechnic du GTPS, 1999. Brest. [24] Serre D. Systems of conservation laws I: Hyperbolicity, entropies, shock waves vol. XXII. Cambridge: Cambridge University Press; 1999, Translated from the French by I.N. Sneddon. [25] Shyue KM. An efficient shock capturing shock––capturing scheme for multicomponent problems. J Comput Phys 1998;142:208–42. [26] Toro EF. Some ivps for which conservative methods fail miserably. In: Proceedings of the 6th International Symposium on Computational Fluid Dynamics, Lake Tahoe, California, 1995. p. 4–8. [27] Yee HC. Construction of explicit and implicit symmetric tvd schemes and their applications. J Comput Phys 1987;68:151–79.