M. Ismail, B. Maury & J.-F. Gerbeau, Editors
NUMERICAL SIMULATIONS OF LOW MACH COMPRESSIBLE TWO-PHASE
FLOWS: PRELIMINARY ASSESSMENT OF SOME BASIC SOLUTION
TECHNIQUES
∗, ∗∗Benjamin Braconnier
1, Jeu-Jiun Hu
2, Yang-Yao Niu
3, Boniface Nkonga
4and
Keh-Ming Shyue
5Abstract. The objective of this paper is to give a preliminary assessment of several existing ap-proaches (compressible or incompressible, preconditioned or not preconditioned) to the numerical sim-ulation of low Mach compressible two-phase flow problems in more than one space dimension. We consider a broken dam problem of Martin and Moyce [Philos. Trans. Roy. Soc. London A, 244(1952), pp. 312-324] and a wave breaking problem of Yasuda et al. [Costal Engineering, 29(1997), pp. 317-346] as our benchmark tests, where the computed results obtained by using the numerical methods can be compared with the laboratory experiments for the basic solution validations.
R´esum´e. L’objectif de ce papier est d’´evaluer la pertinence de plusieurs m´ethodes d´edi´ees `a la simula-tion num´erique d’´ecoulements `a deux phases lorsqu’elles sont utilis´ees pour approcher des ´ecoulements dont le nombre de Mach est faible. Ces m´ethodes reposent sur des mod`eles bidimensionnels gou-vernant des fluides compressibles (avec ou sans pr´econditionnement faible Mach) ou incompressibles. De mani`ere `a permettre une validation ´el´ementaire par comparaison entre r´esultats num´eriques et exp´erimentaux, deux tests nous ont servi de r´ef´erence: l’effondrement d’une colonne d’eau (Martin and Moyce [Philos. Trans. Roy. Soc. London A, 244(1952), pp. 312- 324]) et le d´eferlement d’une vague (Yasuda et al. [Costal Engineering, 29(1997), pp. 317-346]).
1.
Introduction
We are interested in solving unsteady weakly compressible two-phase flow problems where the flow speed is assumed to be much less than the fluid component speed of sound and the wavelengths of the acoustic waves are assumed to be large. Representative applications of this kind of problems are, for instance, the rising of gas bubbles in liquids, the falling of liquid drops in the air under gravitational force, and bubbly flow in liquids, see [13, 37, 48] for more examples.
∗This work was supported in part by a French/Taiwan program ORCHID: Grant of the National science council of Taiwan NSC 97-2911-1-002-015 and EGIDE (France)
∗∗Y.-Y. Niu and K.-M. Shyue were also supported in part by National Science Council of Taiwan Grants NSC 96-2623-7-216-001-D and 96-2115-M-002-008-MY3, respectively
1
GLAIZER Group innovation agency, France, e-mail: [email protected] 2
Department of Information Management, Shu-Te University, Taiwan, e-mail:[email protected] 3
Department of Mechanical and Aerospace Engineering, Chung Hua University, Taiwan, e-mail: [email protected] 4
Department of Mathematics, University of Nice, France, e-mail: [email protected] 5
Department of Mathematics, National Taiwan University, Taiwan, e-mail:[email protected], Corresponding author © EDP Sciences, SMAI 2009
One possible approach to simulate the aforementioned low speed (single- or two-phase) flow problem is to consider fully compressible models and use a standard upwind finite volume method for numerical approxima-tions (cf. [20]). When used in a explicit formulation, it is known in the literature (cf. [11, 34–36, 46, 47, 50]) that the CFL (Courant-Friedrichs-Lewy) stability condition leads to severe time step restriction, yielding difficulties as the lack of robustness of the method and also the loss of accuracy of the computed solutions.
To overcome these difficulties, a common practice in many improved compressible flow solvers is to use splitting. That is, the flow variables are separated into several parts in order to identify the “ill-conditioned” part of the flow variable (typically density or pressure) which may require some special “implicit” treatments. When this is done in a suitable manner, it is possible to bypass the stringent CFL condition dominated by the speed of sound, and use a more robust flow speed-based CFL condition for the computations. This type of methods has been first described by Harlow and Amsden [12], and then generalized by many other researchers (cf. [3, 7, 32, 42, 52] and references therein).
Besides the above approach, Colella and Pao [9] have proposed a projection-type method for which the compressible Euler equations were rewritten in terms of the combination of a Hodge decomposition of the velocity field (i.e., the velocity vector is split into a divergence-free and a curl-free part) and auxiliary pressures. With the new equations, the solution can be separated into an incompressible part which would vary on a time scale determined by the flow speed, and a compressible part that would depend on fast acoustic waves. Hence, the former part may be advanced with time step determined solely by the flow speed, and the latter part may be treated by an implicit method. In addition to that, following the low Mach number asymptotic analysis of Klainermann and Majda [17], Schneider et al. [34] have developed an asymptotic-based method by introducing elliptic constraints in the construction of numerical flux functions. This method yield a class of finite volume compressible flow solvers constituting a viable approach when the Mach number goes to zero. Lastly, in the work of Guillard and coworkers [11,28] a preconditioning technique has been used to construct the flux functions that can also lead to the correct asymptotic behavior of the solution in the zero Mach number limit; this is an extension of the approach given previously by Turkel for steady flows [46, 47].
We note that many of the low speed flow solvers mentioned above are devised for solving single phase compressible flow problems only. The successful extension of the approaches among them to multiphase flow problems include the preconditioned method of Murrone and Guillard [28] and the constrained interpolation profile method of Yabe et al. [52]. Our goal, in this work, is to explore further some of these low speed approaches in the hope to come up with a robust and accurate numerical method for general weakly compressible multiphase flow with real materials.
As a first endeavor towards the method development, we consider a popular broken dam problem of Martin and Moyce [24] as our first benchmark test, and present results obtained by using several numerical approaches based on incompressible or compressible solver for the basic solution validations. Numerical computations for a wave breaking problem of Yasuda et al. [54] is also proposed to demonstrate the accuracy as well as the performance that are attained by using a preconditioned and a not preconditioned compressible solver, see the similar work done by Helluy et al. [13].
The format of this paper is as follows. In Section 2, we review briefly a fluid-mixture type algorithm com-pressible flow solver for homogeneous two-phase flows. In Section 3, we describe a preconditioned comcom-pressible two-phase solver with surface tension effect included in the problem formulation, and in Section 4, we give an incompressible flow solver based on the artificial compressibility approach. Numerical results are presented in Section 5.
2.
A compressible homogeneous two-phase flow solver
We begin our discussion by considering a simplified compressible two-phase flow problem in two space dimen-sions. The two fluid components, consisting of gases and liquids, are non-miscible and separated by interfaces. In this problem, the flow regime of interest is assumed to be homogeneous with no jumps in the pressure and velocity across the interface, and the physical effects such as the viscosity, surface tension, and heat conduction are assumed to be small, and hence can be ignored in the problem formulation. With that, on the gas-phase
part of the domain, the fluid is governed by the compressible Euler equations ∂ ∂t ρ ρu ρv ρE + ∂ ∂x ρu ρu2 + p ρuv ρEu + pu + ∂ ∂y ρv ρuv ρv2 + p ρEv + pv = 0 0 −ρg −ρgv , (1)
where ρ is the density, u and v are the velocities in the x- and y-direction respectively, p is the pressure, E is the total energy per unit mass, and g is the gravitational acceleration. We assume that the equation of state for the gas satisfies the ideal gas law
p = (γ − 1)ρe, (2)
where e is the specific internal energy, γ is the usual ratio of specific heats (γ > 1), and so E = e + (u2 + v2
)/2. Note that the four components of Eqs. (1) express the conservation of mass, momenta in the x- and y-direction, and energy, respectively [10].
On the liquid-phase part of the domain, however, we assume the fluid to be barotropic, and use the the Tait equation of state
p(ρ) = (p0+ B) ρ ρ0
γ
− B (3)
as the constitutive law for the liquid. Here p0 and ρ0 are the reference pressure and density, respectively, and B is a pressure-like constant, see [23, 29, 53] for a typical set of data of practical importance. It is clear that because the pressure is a function of the density only, the energy equation of the Euler Eqs. (1) plays no active role to the determination of this flow behavior, and hence can be neglected. For completeness, we write down the equations of motion for the liquid as follows,
∂ ∂t ρ ρu ρv + ∂ ∂x ρu ρu2 + p(ρ) ρuv + ∂ ∂y ρv ρuv ρv2 + p(ρ) = 0 0 −ρg . (4)
To solve this two-phase flow problem numerically (irrespective of the weakly compressible flow condition), we may use a generalization of the classical shock-capturing method designed originally for single phase flows. In this case, it is known in the literature that the principal difficulty in the usual modification is the occurrence of spurious pressure oscillations when two or more fluid components are present in a grid cell (cf. [1, 15]). To overcome this difficulty, in this section, we take an approach advocated by one of the authors [40, 41] in that thermodynamically the mixture of the fluids characterized by (2) and (3) is assumed to behave like a typical non-barotropic fluid, and can be modeled by the stiffened gas equation of state,
p(ρ, e) = (γ − 1) ρe − γB. (5)
For convenience, we introduce a volume-fraction function α ∈ [0, 1] to represent the type of fluid component within a grid cell. For example, when grid cells contain only the liquid phase we may set α = 1, and so when grid cells contain only the gas phase we set α = 0. When some cells are cut by the interfaces, both of the components are present in the cells so that α ∈ (0, 1) and the volume ratio of each component are respectively given by α and 1 − α. With this definition of α, the pressure of our two-phase flow problem, in all the fluid-component scenarios within a grid cell, can be determined straightforwardly by
p = (p0+ B) ρ ρ0 γ − B if α = 1 (γ − 1) ρe − γB if α 6= 1, (6)
With the pressure law (6), the governing equations for the proposed two-phase flow problems can be written as ∂ρ ∂t + ∂ ∂x(ρu) + ∂ ∂y(ρv) = 0, ∂ ∂t(ρu) + ∂ ∂x ρu 2 + p + ∂ ∂y(ρuv) = 0, ∂ ∂t(ρv) + ∂ ∂x(ρuv) + ∂ ∂y ρv 2 + p = −ρg, ∂ ∂t(ρE) + ∂ ∂x(ρEu + pu) + ∂ ∂y(ρEv + pv) = −ρgv if α 6= 1, ∂α ∂t + u ∂α ∂x + v ∂α ∂y = 0, (7)
where the mixtures of γ and B are computed by
γ = 1 + 1 2 X ι=1 αι γι−1 ! and B= 2 X ι=1 αι γιBι γι−1 ! 1 + 2 X ι=1 αι γι−1 ! .
Note that, in above, we have set α1= α and α2= 1−α for the volume fractions occupied by the fluid-component 1 and 2, in a respective manner. With a system formulated in this way, there is no problem to compute all the state variables of interest, including the pressure from the equation of state. The initialization of the state variables in (7) for fluid-mixture cells can be made in a standard way as described in [38, 39] for numerical simulations.
We note that when the thermodynamical description of the materials of interest is limited by the stability requirement, the fluid speed of sound would belong to a set of real numbers, and so hyperbolicity of this two-phase flow model can be verified by standard characteristic analysis. For the ease of the latter reference, it is useful to write (7) into a more compact expression by
∂q ∂t + ∂f1(q) ∂x + ∂f2(q) ∂y + B1(q) ∂q ∂x + B2(q) ∂q ∂y = ψ(q), (8) with q = [ρ, ρu, ρv, ρE, α]T, f1=ρu, ρu2
+ p, ρuv, ρEu + pu, 0T, f2=ρv, ρuv, ρv2 + p, ρEv + pv, 0T , B1= diag [0, 0, 0, 0, u] , B2= diag [0, 0, 0, 0, v] , ψ = [0, 0, −ρg, −ρgv, 0]T.
If we now take the standard time-splitting version of a compressible flow solver for the numerical approx-imation of this two-phase flow model (see [19, 21] for an unsplit quasi-steady version of the approach to deal with the gravitational source terms considered here), in each time step, we would alternate between solving the homogeneous part of the equations with no source terms
∂q ∂t + ∂f1(q) ∂x + ∂f2(q) ∂y + B1(q) ∂q ∂x + B2(q) ∂q ∂y = 0, (9)
and solving the ordinary differential equations ∂q
∂t = ψ(q). (10)
In the numerical results shown in Section 5, for comparison purposes a couple of the state-of-the-art explicit methods are used to find approximate solutions for (9), while the second order Runge-Kutta method (cf. [16]) is used to find approximate solutions for (10).
To end this section, it is interesting to note that rather than using the Tait equation of state (3) as the constitutive law for the liquid, we may use the stiffened gas equation of state (5) for that instead. In this case, we may simply remove the α 6= 1 constraint from the energy equation in (7), and take the whole set of equations as a governing model for the corresponding two-phase flow problems, see [38] for the details.
3.
A preconditioned compressible two-phase solver with surface tension
Our second physical two-phase flow model assumes again the pressure and velocity equilibrium across inter-faces, and is obtained from an asymptotic limit of the non-equilibrium Bear-Nunziato system [2, 14, 27]. If we include the capillary effect, but ignore the viscous and heat conducting effects as before, the two-dimensional case of the equations can be written analogously as (8) to be of form
∂q ∂t + ∂f1(q) ∂x + ∂f2(q) ∂y + B1(q) ∂w ∂x + B2(q) ∂w ∂y = ψ(q) (11) with q = [α1ρ1, α2ρ2, ρu, ρv, ρE, α1, φ]T, w = [α1ρ1, α2ρ2, u, v, p, α1, φ]T,
f1=α1ρ1u, α2ρ2u, ρu2
+ p, ρuv, ρEu + pu, 0, 0T , f2=α1ρ1v, α2ρ2v, ρuv, ρv2 + p, ρEv + pv, 0, 0T, B1= 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 σκ 0 0 0 0 0 0 0 0 0 0 0 0 0 σκu 0 0 −β 0 0 u 0 0 0 0 0 0 0 u , B2= 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 σκ 0 0 0 0 0 0 0 0 0 0 0 0 0 σκv 0 0 0 −β 0 v 0 0 0 0 0 0 0 v , ψ = [0, 0, 0, −ρg, −ρgv, 0, 0]T.
Here ρι denotes the partial density of the fluid component ι, β is a function that takes into account of the volume fraction variation of the fluids components under pressure fluctuations (cf. [2, 14, 27]), σ is the surface tension coefficient, κ is the curvature at an interface, and φ is the level set function where φ = 0 gives the
approximate location of the interface. It is easy to see that the continuum surface force approach (cf. [4]) has been used in the equations to model the capillary effect of the problem.
To complete the model, we assume that the equation of state associated to each component ι is governed by the a linearized Mie-Gr¨uneisen equation of state [26] as
p(ρι, eι) = (γι−1) (ριeι−ριe∞ ι ) − γιp
∞
ι , (12)
where eι= eι(ρι, p) is the partial internal energy of fluid component ι. Recall that the mixture density, internal energy, and total energy are defined as
ρ = 2 X ι=1 αιρι, ρe = 2 X ι=1
αιριeι, and ρE = ρe +1 2ρ u
2 + v2
,
respectively. In this case, as before, to obtain the pressure equilibrium for multiphase grid cells, the mixtures of γ, p∞
, and e∞
can be obtained from the relations:
1 γ − 1 = 2 X ι=1 αι γι−1, γp∞ γ − 1 = 2 X ι=1 αι γιp ∞ ι γι−1, and ρe ∞ = 2 X ι=1 αιριe∞ ι . (13)
It is known in the literature that the major numerical difficulties for the resolution of this system (11) is related to the non-conservative product term for the compaction of the volume fraction. The mathematical analysis shows that, under the constitutive laws (12) and (13), this model is hyperbolic with real eigenvalues and a complete set of eigenvectors. Note that the equilibrium speed of sound would depend on the asymptotic regime of concern (cf. [14, 25, 27]), and due to the surface tension effect, the pressure would not be continuous across a material interface, but is having a jump that obeys the Laplace law locally [18].
One of the crucial requirement of a multiphase algorithm is its ability to accurately resolve a wide range of multiphase conditions, from fully compressible to the low Mach regimes. As mentioned in the introduction, the low Mach number setting is a singular asymptotic in compressible flows in which the basic feature of the flow is characterized by a large discrepancy between the magnitude of the acoustic wave-speed and the flow velocity, and the usual time-marching solvers for compressible flow is not adequate to compute the weakly compressible or nearly incompressible flows.
In this section, we are interested in a preconditioned-based method for low Mach flow and the preconditioning is obtained by the introduction of an acoustic fluctuation characteristic time τ associated to “incompressible” pressure that is smaller than the mean flow characteristic time t. Time-marching solutions are then designed from pseudo acoustic-waves obtained by altering time-derivatives of the pressure. This is similar to the procedure used in the pseudo-time (or dual-time) [6, 44, 47, 49] and pseudo-compressibility approaches [31]. The detailed description of our approach can be found in [5].
4.
An incompressible two-phase solver
Our third and final physical model for a two-phase flow problem concerns an incompressible flow, and is based on an artificial compressibility approach [8]. In two dimensions, the dimensionless form of the equations
may be written as ∂p ∂t + ∂ ∂x(βu) + ∂ ∂y(βv) = 0, ∂u ∂t + ∂ ∂x u 2 + p + ∂ ∂y(uv) = 1 Re∇ 2 u, ∂v ∂t + ∂ ∂x(uv) + ∂ ∂y v 2 + p = 1 Re∇ 2 v − 1 Frg, ∂α ∂t + ∂ ∂x(uα) + ∂ ∂y(vα) = 0, (14)
where β is the coefficient of artificial compressibility, Re is the Reynolds number, and Fr is the Froude number. To solve (14) numerically, we use a third-order MUSCL (Monotone Upstream-centered Scheme for Conserva-tion Laws) to discretize the convecConserva-tion terms, a second-order finite volume method to discrete the viscous terms, and the two-step Runge-Kutta time integration is used for time marching. For the unsteady problems, we apply the Crank-Nicolson method to the original incompressible Navier-Stokes equations, and then use pseudo-time and the artificial compressibility technique, see [22] for the details. Similar methods of this kind can be found, for example, in the work done by Wackers and Koren [51] and Qian et al. [33].
5.
Numerical results
We now present results for two low Mach test problems obtained by using various numerical approaches described in the previous sections. In order to highlight the validity and the accuracy of the methods, the numerical results are compared with experimental data.
5.1. Broken dam problem
To begin with, we consider the well-known broken dam problem of Martin and Moyce [24]. Here the initial condition consists of a water column with 0.06m wide and 0.12m high that is confined to the left in a rectangular container of size [0, 0.5] × [0, 0.15] m2
with air inside in the standard atmospheric condition, see Fig. 1 for illustration. It is clear that under the effect of the gravity, the water column would collapse subsequently, and the leading water front would move to the right of the domain. Note that in this case the Mach number is typically very small over the time period under consideration.
To solve this problem numerically, we first consider a compressible homogeneous flow model described in Section 2, and employ a standard high-resolution wave propagation method for the homogeneous part of the equations (cf. [20,41]). Numerical results for this problem with three different mesh sizes 100 × 30, 200 × 60, and 400 × 120 are presented in Figs. 1 and 2. That is, to demonstrate the basic feature of the temporal behavior of the water column, Fig. 1 shows the pseudo-color plots of the volume fraction at six different times t = 0, 0.066, 0.109, 0.164, 0.222, and 0.281s, and to validate the computed solutions, Fig. 2 shows the time history of the water column height at the left boundary and also the leading water front position at the bottom boundary, where the numerical results are in comparison with the experimental results provided by Martin and Moyce [24]. From the former figure, the existence of a kink-like structure near the upper-right corner of the water column is clearly seen when the meshes are of sizes 100 × 30 and 200 × 60, and it is not visible as the mesh is refined further. From the latter figure, it is easy to observe that the computed solution with the finest 400 × 120 grid gives the better comparison results on the water column height than the ones with the coarser grids, and somewhat surprisingly the computed solution with the intermediate 200 × 60 grid achieves the better results on the leading front position than the use of the remaining grids. Note that it remains to be an open question now and will be investigated more in the future on the cause of the aforementioned kink-like structure in the solution and also the anomalous convergence behavior of the leading front position.
To carry out the above computation we have used the Courant number 0.9, and employed solid wall boundary conditions to all the boundaries. The material constants are for the air: ρ = 1.2 kg/m3, γ = 1.4, and for the water: ρ0 = 103 kg/m3, p0 = 1 atm, γ = 7, and B = 3000 atm. The gravitational constant have been set to g = 9.81 m/s2.
Figure 1. Pseudo-color plots of volume fraction for the broken dam problem obtained by using a compressible homogeneous model described in Section 2, and a high-resolution wave propagation method. Results are shown at six different times t = 0, 0.066, 0.109, 0.164, 0.222, and 0.281s with three different mesh points 100 × 30, 200 × 60, and 400 × 120.
Before presenting numerical results obtained by using the other mathematical models, we pause to examine the effect of the equation of state for the water in this compressible two-phase flow model, and for comparison purposes, we also perform the computations by using the AUSM-based third order ENO and WENO meth-ods [30]. Figure 3 gives the validation results for the time history of the water column height at the left boundary and also the leading water front position at the bottom boundary, where the stiffened gas equation of state (5) is used for the water with γ = 4.4, B = 6000 atm, and the grid size is 400 × 120. From the figure, it is easy to observe that the equation-of-state effect is not significant for this problem, and qualitatively there is not much difference between the numerical results.
0 0.5 1 1.5 2 2.5 3 3.5 0.2 0.4 0.6 0.8 1 Experiment 100 × 30 200 × 60 400 × 120 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 1 2 3 4 5 6 Experiment 100 × 30 200 × 60 400 × 120 tpg/a y / (2 a ) t/p2g/a x / a
Water column height
Leading water front position
Figure 2. Comparison between numerical solutions shown in Fig. 1 and experimental results presented in Martin and Moyce [24]. On the top, results on the water column height are shown, and on the bottom, results on the leading water front position are shown. A mesh-refinement sequence 2i(100 × 30) for i = 1, 2, 3 are used for comparison. Here we have a = 0.06m and g = 9.81 m/s2
.
Our second attempt to solve this problem numerically is to employ a preconditioned compressible flow solver described in Section 3 which is an implicit method with Courant number 20 and Gauss-Seidel solver (50 iterations maximum) included. As before, in Fig. 4, we validate the solutions by comparing them to the experimental results on the water column height at the left boundary and also the leading water front position at the bottom boundary. Unlike the unpreconditioned case shown in Fig. 2, reasonable agreement of the results are observed on both the height and the front location as the mesh is refined. It should be mentioned that from quite a few numerical evidences (not shown here) the surface-tension effect at the air-water interface does not play an important role at the early stage of the motion of this collapse of water column. It would be a significant effect however right after the eventual breaking of the water column at a later time due to the wave collision to the container’s right and top boundaries. Note that the material constants we use here are for the air: ρ = 1 kg/m3
, γ = 1.4, p∞
= 0, e∞
= 0, and for the water: ρ = 103 kg/m3
, γ = 4.4, p∞
= 6000 atm, and e∞ = 0.
Our third and final attempt to solve this problem numerically is to employ an incompressible flow solver based on the artificial compressibility approach described briefly in Section 4. In the method, a dual-time strategy is used with a sub-iteration step included to achieve the divergence-free flow condition at every physical time step. In Fig. 5, we present the validation results on the water column height at the left boundary and also the leading water front position at the bottom boundary. with a 200 × 60 mesh points. It is interesting to see that this incompressible flow solver works rather good for this broken dam problem. Note that results obtained by
0 0.5 1 1.5 2 2.5 3 3.5 0.2 0.4 0.6 0.8 1 Experiment Wave ENO WENO 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 1 2 3 4 5 6 Experiment Wave ENO WENO tpg/a y / (2 a ) t/p2g/a x / a
Water column height
Leading water front position
Figure 3. Numerical results for the study of the equation-of-state effect on the broken dam problem with a 400×120 grid. The stiffened gas equation of state (5) is used in the compressible homogeneous model described in Section 2 when we carry out the computation by using the AUSM-based third order ENO and WENO methods, while as in Fig. 2 the Tait equation of state (3) is used with the wave propagation method.
using compressible flow solvers as shown in Figs. 2 and 4 are included in the graph also as for comparison. Here the mesh point is 200 × 60. Here the parameters we use in the run for the method are β = 1, Re = 200, Fr = 1, and the density ratio between the water and the air is 1000.
Table 1 gives the computational cost in CPU time for all the runs presented in this subsection. It is clear that the preconditioned compressible flow solver described in Section 3 gives a somewhat better performance than the other two numerical approaches. It is also interesting to see that the standard compressible flow solver (with wave propagation method) works faster than the AUSM-based method and the incompressible flow solver, when viewing on the same level of the grid size.
5.2. Wave breaking problem
We are concerned with a wave breaking problem that has been studied extensively in the literature, see Helluy et al. [13] and the references cited therein. The initial condition is a rightward-moving solitary wave which can be computed by the Tanaka method [45] traveling toward a step-like reef on the right, see Fig. 6 for the schematic setup, and see [13] for more details.
As a first attempt to solve this problem numerically, we use the same explicit compressible flow solver with the high-resolution wave propagation method as before for the broken dam problem, and we present the results in Figs. 7 and 8 where the programs are run with three different mesh sizes 400 × 50, 800 × 100, and 1600 × 200.
0 0.5 1 1.5 2 2.5 3 3.5 0.2 0.4 0.6 0.8 1 Experiment 100 × 30 200 × 60 400 × 120 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 1 2 3 4 5 6 Experiment 100 × 30 200 × 60 400 × 120 tpg/a y / (2 a ) t/p2g/a x / a
Water column height
Leading water front position
Figure 4. Comparison between numerical solutions obtained by using a preconditioned com-pressible flow solver described in Section 3 and experimental results for the broken dam problem. The graph of the solutions are displayed in the same manner as in Fig. 2.
Table 1. Computational cost for a broken dam problem measured in CPU time
Method Mesh CPU time (sec) CPU type
Compressible solver(wave-based) 100 × 30 492 AMD Dural-Core Opteron
200 × 60 3782 2220 2.8GHz
400 × 120 31783
Compressible solver(AUSM-ENO) 400 × 120 105200 IBM P595 Compressible solver(AUSM-WENO) 400 × 120 170000
Precond. compressible solver 100 × 30 352 Intel Core 2 Duo 3.0GHz 200 × 60 2453
400 × 120 21780
Incompressible solver 200 × 60 9804 Intel Pentium 4 3.4GHz
From the pseudo-color plots of the volume fraction shown in Fig. 7 at four different times t = 1.2, 1.35, 1.6, and 1.8s, the breaking of the water wave is clearly seen as the mesh is refined. To validate the computed solutions, Fig. 8 shows the time history of the elevation of the water surface at the three different gauges P2, P3, and
0 0.5 1 1.5 2 2.5 3 3.5 0.2 0.4 0.6 0.8 1 Experiment Incompressible solver Compressible solver Precond. Compressible solver
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 1 2 3 4 5 6 Experiment Incompressible solver Compressible solver Precond. Compressible solver
tpg/a y / (2 a ) t/p2g/a x / a
Water column height
Leading water front position
Figure 5. Comparison between numerical solutions obtained by using an incompressible flow solver described in Section 4 and experimental results for the broken dam problem. Results obtained by using compressible flow solvers as shown in Figs. 2 and 4 are included in the graph also. Here the mesh point is 200 × 60.
P4 which are located at x = 2m, 2.515m, and 3.02m, respectively. Note that for this problem, experimental time evolution of the water elevations are available athttp://helluy.univ-tln.fr/soliton.htm. From the graph, it is easy to notice the sensible agreement of the results as the mesh turns finer.
We continue solving this wave breaking problem by the same preconditioned compressible flow solver as done before. Figure 9 gives the Schlieren-type images of volume fraction at four different time t = 1.2, 1.35, 1.6, and 1.8s where a 1500 × 200 mesh points is used in the computation. As far as the qualitative global free-surface structure is concerned, we notice reasonable resolution of the solution as in comparison with the results shown in Fig. 7 (the 1600 × 200 grid case) and also in [13]. To validate the results, in Fig. 10, we again present the time history of the elevation of the water surface at the three different gauges P2, P3, and P4, noticing good results with the experimental data and also the results obtained by using the explicit compressible solver.
Table 2 gives the computational cost in CPU time for all the runs presented in this subsection, and also the compressible and incompressible results reported in [13] as for comparison. It is clear that the implicit preconditioned compressible flow solver now out-performs the explicit wave-propagation based compressible solver by a factor of 3.5 approximately, and is the most robust one among all the methods considered here.
6.
Conclusion
To assess the feasibility of the existing approaches for numerical resolution of low Mach two-phase flow problems, in this work, we have taken three popular solution techniques into consideration, that is, a standard
Figure 6. Initial setup for the wave breaking problem.
Figure 7. Pseudo-color plots of volume fraction for the wave breaking problem obtained by using the same explicit compressible flow solver as shown in Fig. 1 for the broken dam problem. Results are shown at four different times t = 1.2, 1.35, 1.6, and 1.8s with three different mesh points 400 × 50, 800 × 100, and 1600 × 200.
explicit compressible solver, an implicit preconditioned compressible solver, and an implicit incompressible solver based on artificial compressibility. Numerical results presented here for the two benchmark tests; the broken dam and wave breaking problems, indicate good potential to achieve highly accurate results by many of the aforementioned methods. If efficiency of the method is taken into account, the implicit preconditioned
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.05 0.1 0.15 0.2 Experiment 400 × 50 800 × 100 1600 × 200 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.05 0.1 Experiment 400 × 50 800 × 100 1600 × 200 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.05 0.1 0.15 Experiment 400 × 50 800 × 100 1600 × 200 Gauge P2 Gauge P3 Gauge P4 W a te r su rf a ce (m ) W a te r su rf a ce (m ) W a te r su rf a ce (m ) t(sec)
Figure 8. Comparison between the computed water surface for the solutions shown in Fig. 7 at three different gauges with experimental results presented in Yasuda et al. [54]. A mesh-refinement sequence 2i(400 × 50) for i = 1, 2, 3 are used for comparison.
Table 2. Computational cost for a wave breaking problem measured in CPU time
Method Mesh CPU time (sec) CPU type
Compressible solver(wave-based) 400 × 50 6756 AMD Dural-Core Opteron
800 × 100 55253 2220 2.8GHz
1600 × 200 476429
Compressible solver(MUSCL-based) [13] 1500 × 200 172800 Alpha 666 MHz 64 bits Precond. compressible solver 1500 × 200 146400 Intel Quad Core Xeon 3.0GHz Incompressible solver [13] 1200 × 200 273600 Itanium 1.4GHz 64 bits
compressible solver is a better choice of the method which is as expected. Future work is to examine more implicit strategies with preconditioned or not preconditioned in the hope to devise a robust method for a wide variety of low Mach multiphase flow problems of practical importance.
Figure 9. Schlieren-type images of volume fraction for the wave breaking problem obtained by using the same implicit preconditioned compressible flow solver as shown in Fig. 4 for the broken dam problem. Results from the top to bottom are shown at four different times t = 1.2, 1.35, 1.6, and 1.8s with mesh points 1500 × 200.
References
[1] R. Abgrall. How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach. J. Comput. Phys., 125:150–160, 1996.
[2] M. R. Baer and J. W. Nunziato. A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials. Int. J. Multiphase Flow, 12(6):861–889, 1986.
[3] H. Bijl and P. Wesseling. A unified method for computing incompressible and compressible flows in boundary-fitted coordinates. J. Comput. Phys., 141:153–173, 1998.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.05 0.1 0.15 0.2 Experiment
Precond. compressible solver Compressible solver 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.05 0.1 Experiment
Precond. compressible solver Compressible solver 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.05 0.1 0.15 Experiment
Precond. compressible solver Compressible solver Gauge P2 Gauge P3 Gauge P4 W a te r su rf a ce (m ) W a te r su rf a ce (m ) W a te r su rf a ce (m ) t(sec)
Figure 10. Comparison between the computed water surface for the solutions shown in Fig. 9 at three different gauges with experimental results presented in Yasuda et al. [54]. Results obtained by using compressible flow solvers as shown in Figs. 8 is included in the graph also. Here the mesh point is 1500 × 200.
[4] J. U. Brackbill, D. B. Kothe, and C. Zemach. A continuum method for modeling surface tension. J. Comput. Phys., 100:335– 354, 1992.
[5] B. Braconnier and B. Nkonga. An all-speed relaxation scheme for interfacial flow with surface tension. Preprint, 2008. [6] P. E. O. Buelow, S. Venkateswaran, and C. L. Merkle. Stability and convergence analysis of implicit upwind schemes. Computers
& Fluids, 30:961–988, 2001.
[7] V. Casulli and D. Greenspan. Pressure method for the numerical solution of transient, compressible fluid flows. Int. J. Numer. Meth. Fluids, 4:1001–1012, 1984.
[8] A. J. Chorin. A numerical method for solving incompressible viscous flow problem. J. Comput. Phys., 2:12–26, 1967. [9] P. Colella and K. Pao. A projection method for low speed flows. J. Comput. Phys., 149:245–269, 1999.
[10] R. Courant and K. O. Friedrichs. Supersonic Flow and Shock waves. Wiley-Interscience, New York, 1948.
[11] H. Guillard and C. Viozat. On the behavior of upwind schemes in the low Mach number limit. Computers & Fluids, 28:63–86, 1999.
[12] F. H. Harlow and A. A. Amsden. Numerical calculation of almost incompressible flow. J. Comput. Phys., 3:80–93, 1968. [13] P. Helluy, F. Golay, J.-P. Caltagirone, P. Lubin, S. Vincent, D. Drevard, R. Marcer, P. Fraunie, N. Seguin, S. Grilli, A.-C.
Lesage, A. Dervieux, and O. Allain. Numerical simulations of wave breaking. ESAIM: M2AN, 39:591–607, 2005.
[14] A. K. Kapila, R. Menikoff, J. B. Bdzil, S. F. Son, and D. S. Stewart. Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations. Phys. Fluids, 13(10):3002–3024, 2001.
[15] S. Karni. Multicomponent flow calculations by a consistent primitive algorithm. J. Comput. Phys., 112:31–43, 1994. [16] D. Kincaid and W. Cheney. Numerical Analysis. Brooks/Cole, 1990.
[17] S. Klainermann and A. Majda. Compressible and incompressible fluids. Comm. Pure Appl. Math., 35:629–651, 1982. [18] L. D. Landau and E. M. Lifshitz. Fluid Mechanics. Pergamon, New York, 1959.
[19] R. J. LeVeque. Balancing source terms and flux gradients in high-resolution godunov methods: The quasi-steady wave propa-gation algorithm. J. Comput. Phys., 146:346–365, 1998.
[20] R. J. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, 2002.
[21] R. J. LeVeque and D. S. Bale. Wave propagation method for conservation laws with source terms. In Hyperbolic Problems: Theory, Numerics, Applications, pages 609–618. Birkh¨auser-Verlag, 1999. Intl. Series of Numerical Mathematics, Vol. 130. [22] S.-Y. Lin and T.-M. Wu. An adaptive multigrid finite volume scheme for incompressible Navier-Stokes equations. Int. J.
Numer. Meth. Fluids, 17:687–710, 1993.
[23] S. P. Marsh. LASL Shock Hugoniot Data. University of California Press, Berkeley, 1980.
[24] J. C. Martin and W. J. Moyce. An experimental study of the collapse of liquid columns on a rigid horizontal plane. Philos. Trans. Roy. Soc. London A, 244:312–324, 1952.
[25] J. Massoni, R. Saurel, B. Nkonga, and R. Abgrall. Proposition de m´ethodes et modeles euleriens pour les problems a interfaces entre fluides compressibles en presence de transfert de chaleur. Intl. J. Heat and Mass Transfer, 45(6):1287–1307 (in french), 2002.
[26] R. Menikoff and B. Plohr. The Riemann problem for fluid flow of real materials. Rev. Mod. Phys., 61:75–130, 1989.
[27] A. Murrone and H. Guillard. A five reduced equation model for compressible two-phase flow problems. J. Comput. Phys., 202:664–698, 2004.
[28] A. Murrone and H. Guillard. Behavior of upwind scheme in the low Mach number limit: III. preconditioned dissipation for a five equation two phase model. Computers & Fluids, 37:1209–1224, 2008.
[29] H. Nagoya, T. Obara, and K. Takayama. Underwater shock wave propagation and focusing in inhomogeneous media. In R. Brun and L. Z. Dumitrescu, editors, Proceedings of 19th Intl. Symp. on Shock Waves, Marseille, pages 439–444. Springer-Verlag, Berlin, 1995.
[30] Y.-Y. Niu. A simple and robust advection upwind flux splitting to simulate transient cavitated water-vapor flows. Numerical Heat Transfer: Part B, 51(7):679–696, 2007.
[31] R. R. Nourgaliev, T. N. Dinh, and T. G. Theofanous. A pseudo-compressibility method for the numerical simulation of incompressible multifluid flows. Int. J. Multiphase Flow, 30:901–937, 2004.
[32] G. Patnaik, R. H. Guirguis, J. P. Boris, and E. S. Oran. A barely implicit correction for flux-correct transport. J. Comput. Phys., 71:1–20, 1987.
[33] L. Qian, D. M. Causon, D. M. Ingram, C. G. Mingham. Cartesian cut cell two-fluid solver for hydraulic flow problems. J. Hydraulic Eng., 129:688–696, 2003.
[34] T. Schneider, N. Botta, K. J. Geratz, and R. Klein. Extension of finite volume compressible flow solvers to multi-dimensional, variable density zero Mach number flows. J. Comput. Phys., 155:248–286, 1999.
[35] J. Sesterhenn, B. Mueller, and H. Thomann. Flux-vector splitting for compressible low Mach number flow. Computers and Fluids, 22:441–451, 1993.
[36] J. Sesterhenn, B. Mueller, and H. Thomann. On the cancellation problem in calculating compressible low Mach number flows. J. Comput. Phys., 151:597–615, 1999.
[37] B. R. Shin, S. Yamamoto, and X. Yuan. Application of preconditioning method to gas-liquid two-phase flow computations. J. Fluids Eng., 126:605–612, 2004.
[38] K.-M. Shyue. An efficient shock-capturing algorithm for compressible multicomponent problems. J. Comput. Phys., 142:208– 242, 1998.
[39] K.-M. Shyue. A fluid-mixture type algorithm for compressible multicomponent flow with van der Waals equation of state. J. Comput. Phys., 156:43–88, 1999.
[40] K.-M. Shyue. A volume-of-fluid type algorithm for compressible two-phase flows. In Hyperbolic Problems: Theory, Numerics, Applications, pages 895–904. Birkh¨auser-Verlag, 1999. Intl. Series of Numerical Mathematics, Vol. 130.
[41] K.-M. Shyue. A volume-fraction based algorithm for hybrid barotropic and non-barotropic two-fluid flow problems. Shock Waves, 15(6):407–423, 2006.
[42] K.-M. Shyue. Mach-uniform methods for compressible two-phase flow problems: A preliminary report. Term report: NSC94-2115-M-002-016, National Science Council, Taiwan, R.O.C., 2006 (unpublished).
[43] K.-M. Shyue. A wave propagation method for compressible multicomponent problems on quadrilateral grids. Manuscript sub-mitted for publication, 2009.
[44] R. C. Swanson and E. Turkel. Pseudo-time algorithms for Navier-Stpkes equations. Appl. Num. Math., 2:321–333, 1986. [45] M. Tanaka. The stability of solitary waves. Phys. Fluids, 29:650–655, 1986.
[46] E. Turkel. Preconditioned methods for solving the incompressible and low speed compressible equations. J. Comput. Phys., 72:277–298, 1987.
[47] E. Turkel. Review of preconditioning techniques for fluid dynamics. Appl. Num. Math., 12:257–284, 1993.
[48] D. R. van der Heul, C. Vuik, and P. Wesseling. A staggered scheme for hyperbolic conservation laws applied to unsteady sheet cavitation. Comput. Visual Sci., 2:63–68, 1999.
[49] S. Venkateswaran, J. W. Lindau and. R. F. Kunz, and C. L. Merkle. Computation of multiphase mixture with compressibility effects. J. Comput. Phys., 180:54–77, 2002.
[50] G. Volpe. Performance of compressible flow codes at low Mach numbers. AIAA J., 31(1):49–56, 1993.
[51] J. Wackers and B. Koren. A surface capturing method for the efficient computation of steady water waves. J. Comput. and Appl. Math., 215:618–625, 2008.
[52] T. Yabe, F. Xiao, and T. Utsumi. The constrained interpolation profile method for multiphase analysis. J. Comput. Phys., 169:556–593, 2001.
[53] K. Yamada, H. Nagoya, and K. Takayama. Shock wave reflection and refraction over a two-fluid interface. In R. Brun and L. Z. Dumitrescu, editors, Proceedings of 19th Intl. Symp. on Shock Waves, Marseille, pages 299–304. Springer-Verlag, Berlin, 1995.
[54] T. Yasuda, H. Mutsuda, and N. Mizutani. Kinematic of overturning solitary waves and their relations to breaker types. Coastal Engineering, 29:317–346, 1997.
![Figure 2. Comparison between numerical solutions shown in Fig. 1 and experimental results presented in Martin and Moyce [24]](https://thumb-ap.123doks.com/thumbv2/9libinfo/8684387.197561/9.918.159.694.183.627/figure-comparison-numerical-solutions-experimental-results-presented-martin.webp)
