A simple and efficient finite-volume method is presented for the computation of compressible flows of two immiscible fluids at very different densities. One novel ingredient in the method is a two-fluid Osher-type scheme, which is capable of computing the cell-face flux in case of two different fluids (e.g., water and air) left and right of the cell face. The other original property of the method is that a level-set term, for distinguishing between the two fluid8, is consistently incorporated as one of the flux components. The level-set flux is properly treated by the Osher-type scheme

Brummelen, E.H. (Harald) van

Koren, B. (Barry)

Leer, B. van

Lewis, M.R. (Mervyn)

CWI's Institutional Repository

An Osher-Type and Level-Set Scheme for Two-Fluid Flow Computations Barry Koren* - Mervyn Lewis* -Harald van Brummelen** -Bram van Leer*** *CWI, P.O. Box 94019, 1090 GB Amsterdam, The Netherlands. Barry.Koren@cwi.nl, Mervyn.Lewis@cwi.nl **Delft University of Technology, P.O. Box 5058, 2600 GB Delft, The Netherlands. e.h. vanbrummelen@lr. tudeljt. nl ***The University of Michigan, Ann Arbor, MI 48109-2140, USA. bram@engin. umich. edu ABSTRACT. A simple and efficient finite-volume method is presented for the compu-tation of compressible flows of two immiscible fluids at very different densities. One novel ingredient in the method is a two-fluid Osher-type scheme, which is capable of computing the cell-face flux in case of two different fluids (e.g., water and air) left and right of the cell face. The other original property of the method is that a level-set term, for distinguishing between the two fluids, is consistently incorporated as one of the flux components. The level-set flux is properly treated by the Osher-type scheme. KEYWORDS: compressible liquid-gas flows, interface capturing, Osher-type scheme, level-set method. 1. Introduction The present paper is directed towards an efficient, physically correct finite-volume computation of two-fluid flows consisting of, e.g., compressible water and air, at uniformly subsonic (preferably low-subsonic) speeds. Other premises are that the two fluids do not mix, that vaporization and condensation phenom-ena do not occur and that surface tension can be neglected. Last years, various 584 Finite volumes for complex applications papers have been published in which specific finite-volwne methods are pre-sented for two-fluid fl.ow computations. In most of these papers, a two-fluid flux formula is proposed, in which an approximate Riemann solver is applied. E.g., in [KAR 94], a Roe-type scheme is proposed, as is in [ABG 96]. In [SAU 99], Rusanov-type and Harten-Lax-van-Leer-type schemes are presented. Not yet seen in a two-fluid context is: an Osher-type scheme. (Note that the two--fluid method to which Osher himself contributed in [FED 99], the ghost fluid method, needs a single-fluid flux formula only.) In the present paper, we pro-pose a two-fluid Osher-type flux formula. The formula is extremely simple and computationally very efficient. At the boundaries of the computational domain, it is completely consistent with that in the interior. Also of interest is that the Osher-type flux formula incorporates a level-set tenn for accurately capturing the two-fluid interface. The 'pressure-oscillation problem', which is addressed and fixed in most of the aforementioned literature, is not considered in the present paper. Fixes that have already been proposed for it can be simply carried over to the present method. 2. Flow model 2.1. Conservation equations In lD, for a sufficiently small control volume n, conservation of mass and momentum read: f d ( P ) dx + ( pu ) ( pu ) - 0 [1] J S1 dt pu pu2 + p Bnrlght - pu2 + p &n1ert - 1 with p the bulk density [2] where a is, e.g., the volume-of-water fraction, and where Pw(p) and p0 (p) are the equations of state for water and air, respectively. To balance system [1 ]-(2], the latter equations and an equation for the location of the interface ( deter-mining a) still have to be chosen. With the bulk-density formulation, in a finite-volwne discretization, total mass of the fluid will be conserved, but not necessarily the masses of the two separate fluids. In case a(x, t) is poorly resolved, the two separate masses are poorly conserved as well. Hence, an accurate resolution of the interface location(s) is of paramount importance. For this purpose, we follow a level-set approach, to be discussed in the next section. Applications 585 2.2. Level-set equation To accurately resolve the interface location(s), a level-set approach [SET 96] is more appropriate than the classical volume-of-fluid (VOF) approach [HIR 81 ], because of its better smoothness properties at precisely the point of interest: the interface. Good smoothness of the level-set function is first taken care of in the level-set function's initialization. A common approach is to initialize the level-set function as the signed distance to the initial interface, with the distance positive in, e.g., water and negative in air. (Nonlinear initializations are possible [KOR 99].) Denoting the level-set function by ¢, it is advected by, in lD: 8¢ 8¢ -+u-=0. at ax [3] Combined with the bulk-mass conservation equation from [1], quasi-linear equa-tion [3] may be written in the conservative control-volume form { d(p¢) Jo. ---;ftdx + (pucp)ao.,iah• - (pu¢)ao.1er, = 0. [4] Conservation of pip is not physical, there is no conservation law for it. The form [4] is simply practical because it is consistent with system [1]; it can be directly embedded into it. With ¢(x, t) known, for any finite volume the VOF-function a can be computed. In Section 3.1, for an equidistant finite-volume grid, a= a(¢) is worked out in detail. So, note that the VOF-fraction is used, but not the VOF-method. In the VOF-method, a transport equation for a is used. Instead, here we apply a transport equation for¢ (in the consistent form of a conservation equation for p</J). 2.3. Equation of state In, e.g., water-air computations, for both fluids, elegant use can be made of a single equation of state, Tait's: P + BPref ( P ) 'Y (1 + B)Pref = Pref ' [5] where the subscript 'ref' indicates some reference state. For water, it holds: 1 = 7, B = 3000 and for air: 1 = i, B = 0. With [5], both the water and air density, Pw(p) and p0 (p), are convex functions of pressure. Likewise, the corresponding bulk density p according to [2] is. The physical consequences of this overall convexity are that neither locally low speeds of sound (lower than in pure water or pure air), nor entropy-condition-satisfying expansion shocks 586 Finite volumes for complex applications can occur. A slight inconvenience of [2) in combination with [5] is that the calculation of p for known p and a ( a ¥- 0 and a I 1) needs to be done iteratively. 3. Discretization 3.1. Finite volumes Summarizing, for a (sufficiently small) control volume n, the system of equations considered reads { dq Jn dt dx+(f(q))an,;gh, -(J(q))&n1e£, = 0, ·-(6), f (q) = ( pur-+ p ) , purp P = a(¢)Pw(P) + (1 - a(cp))Pa(P), ...L ..L (6) [7) ( ) ( P + BwPref ) '"' ( ) Pw p = (I+ Bw)Pref Pw ref, ( P + BaPref ) -,., Pa(P) = (l + Ba)Pref (Pa)ref, [8] with a.(¢) the fraction of the size of n over which¢~ 0. The natural space discretization for [6] is a finite-volume technique. We consider cell-centered finite volumes with, for convenience, constant mesh size. This choice directly allows us to work out the discretization of a.(¢). Consider finite volume fl; and its left and right neighbors, ni-l and ni+i, respectively, and define the level-set values at the cell faces an;_i and afli+i as 2 2 [9] Then, for r/Ji ~ 0, we propose the following expression for a.;: [10] 1 ( <Pi ) Oi = 2 cp; - <Pi-½ + 1 ' (11] ¢;_l.2 ~ 0, ¢;+1.2 < 0 : a.i = ! (1 + ¢; ) ' [12] 2 <Pi - <Pi+½ <P;-l. < 0, ¢;+1. < 0: a; = ~ ( ¢; + <Pi ) . [13) 2 2 2 ¢i - <Pi-½ <Pi - ¢;+½ So, in determining ¢;_1. and ¢;+1., as well as x( ¢ = 0), use is made of piecewise 2 2 linear interpolation of ¢. The linear interpolation is exact as long as the level-set function is the signed-distance function. Applications 587 3.2. Riemann-problem approach For the finite-volume formulation, we need a formula for the flux vector across a cell face. The formula must have built-in physics for accurately cap-turing the interface. The exact solution of the lD Riemann problem on each cell face, the well-known Godunov approach, requires the use of a numerical root finder. We avoid this by considering an approximate Riemann solver. For this, we prefer Osher's [OSH 82], particularly because of its consistent boundary-condition treatment. Denoting the left and right cell-face state by q0 and q1 , and the flux formula by F(q0 , q1 ), the Osher scheme may be written as [14] with ~ the negative eigenvalue part of i. The eigenvalues of the present Jacobian are: >-1 = u --JJi, .-\2 = u, >.3 = u+ Jfi. ( ~ does not occur in the wave speeds!) The Riemann-invariant relations describing the two intermediate states ~½ and qi along the wave path in state space are for >.1 , .,\2 and >.3 successively: (15] (16] / Pf 1 ~p !Pi l~p u ... - - -dp = U1 - - -dp, s P op P op [17] Hence, the level-set function cp can only change along the subpath corresponding with the eigenvalue >.2 , i.e., across the contact discontinuity. It is invariant along the outer subpaths; physically speaking: along the latter the distance to the two-fluid interface is constant. The integrals in [15] and [17] can be written out explicitly for the equations of state given in [8]. However, when, e.g., a water-air interface is captured, explicit calculation of u1. and Pl. is hampered 2 2 by nonlinearity; a transcendental equation needs to be solved. 3.3. Linearized two-fluid Osher scheme Since cp is constant along the two outer subpaths of the wave path, along both subpaths the bulk density can only vary due to pressure changes. For flows that are low-subsonic, large density changes will not occur there and -588 Finite volumes for complex applications consequently - the integrals rt ¼~dp and rt ¼./fidp can be linearized by good approximation around Po and P1, respectively, yielding for [15]- [17]: Co u1 = u 0 - {p1 - Po)-, 2 a Po Likewise, P½ can be linearized around Po and P1: P.1. =Po+ (p1 - Po)~, 2 8 (18] [19] Pl. -Po Elimination of P½ - Po and Pi - Pt from (18] and [19] gives uro = -poeo p½-Pl . and = P1C1, 1.e., u½-u1 Co= poco, [20] For the density and level-set function in the two intermediate points it holds: [21] Ignoring all supersonic possibilities among all possible combinations of signs of u 0 - c0, u.i.. - c.1., u1 + c.a and u 1 + c1 {note the consequent, large improvement ~ a 2 s in efficiency), the linearized two-fluid scheme reads then: We remark that the real nonlinear flux functions f(qt) and !(qi) are applied, and not F(qo,q1) = f(qo) + (q½ - qo)d1J:o) if u½ 2:: 0 and F(qo, qi)= f(q1) + (q¾ - q1) dfJ!1 ) if u½ ~ 0. There is no need for the latter linearized formulae. On the contrary, as opposed to [20]-(22], they may give rise to an erroneous, ambiguous flux at u½ = 0 (steady contact discontinuity); f (qo) + (q½ -q0 ) 4f~o) and J(q1) + (q~ - q1)dJJ:1 ) may be different for u½ = 0. 3.4. Boundary-condition treatment A very favorable property of the Osher scheme is that the fluxes across the boundary faces can be computed with the same formula as that for the interior Applications 589 faces, i.e., with [22] as well. Denoting the state at the boundary by qb, in case of a left boundary qo = qb and in case of a right q1 = qb, We work out the inflow and outflow boundary, and the non-permeable boundary as a limit case. For all three it holds, for boundary at the left and right, respectively: 3.4.1. Inflow boundary Pb-Po= -Co. Ub-Uo [23] From [23], it follows that the two boundary conditions to be imposed here cannot be Ub and Pb simultaneously; when ub is imposed, Pb follows - vice versa - when Pb is imposed, Ub follows. Hence, the second boundary condition must be one for <Pb· To compute the corresponding boundary flux f(qb) = (Pb ub, Pb u~ + Pb, Pb<Pb f, the 'OD' bulk density still needs to be defined. (In 2D and 3D, the bulk density can be computed in a normal lD and 2D way, respectively.) In lD, an appropriate 'OD' choice is [24] 3.4.2. Outflow boundary Here, in addition to [23], the equations [25] are available. So, the single boundary condition to be imposed must be Ub or Pb, or some combination of both. The bulk density Pb is defined as in the inflow case. 3.4.3. Non-permeable boundary At a non-permeable boundary ( at least) Ub = 0 must be imposed, which, given [23], already determines Pb• Considering a non-permeable boundary as the limit case of an inflow boundary, <Pb must still be imposed. Considering it as the limit of outflow, <Pb follows from the interior solution ( <Pb = ¢1 for left boundary and </Jb = ¢0 for right). The outflow-limit case is to be preferred. As opposed to in the inflow-limit case, it allows the interface to freely move along the non-permeable boundary. Also here, the bulk density may be defined according to [24]. 4. Conclusions To accurately compute compressible, immiscible two-fluid flows with very large density differences (water-air flows, e.g.), we have proposed a method that uses a level-set technique to distinguish between the two fluids. The level-set equation has been incorporated consistently into the system of hyperbolic 590 Finite volumes for complex applications conservation laws. The resulting equations have been discretized through a finite-volume method. To compute the fluxes across the finite-volume walls (the level-set flux being one of the flux-vector components), we have proposed a two-fluid, linearized Osher scheme (given by [22]). The scheme allows a physically correct capturing of the interface across a single cell face, as well as a neat boundary-condition treatment (no sticking of interfaces to solid walls, e.g.). The novel scheme combines good physical properties with great simplic-ity and efficiency. There are no principal difficulties in fixing the technique for the 'pressure-oscillation problem', and in extending it to multi-D. NOTE. - This research was supported by the Dutch Technology Foundation STW. References [ABG 96] ABGRALL, R., "How to Prevent Pressure Oscillations in Multicom-ponent Flow Calculations: a Quasi Conservative Approach", J. Comput. Phys., vol. 125, 1996, pp. 150-160. [FED 99] FEDKIW, R.P., ASLAM, T., MERRIMAN, B., OSHER, S., "A Non-Oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method)", J. Comput. Phys., vol. 152, 1999, pp. 457-492. [HIR 81] HIRT, C.W., NICHOLS, B.D., "Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries", J. Comput. Phys., vol. 39, 1981, pp. 201-225. [KAR 94] KARNI, S., "Multicomponent Flow Calculations by a Consistent Primitive Algorithm", J. Comput. Phys., vol. 112, 1994, pp. 31-43. [KOR 99] KOREN, B., VENIS, A.C.J., "A Fed Back Level-Set Method for Moving Material-Void Interfaces", J. Comput. Appl. Math., vol. 101, 1999, pp. 131-152. [OSH 82] OSHER, S., SOLOMON, F., "Upwind Difference Schemes for Hyper-bolic Systems of Conservation Laws", Math. Comput., vol. 38, 1982, pp. 339-374. [SAU 99] SAUREL, R., ABGRALL, R., "A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows", J. Comput. Phys., vol. 150, 1999, pp. 425-467. [SET 96] SETHIAN, J .A., Level Set Methods: Evolving Interfaces in Geom-etry, Fluid Mechanics, Computer Vision, and Materials Science, Cam-bridge University Press, 1996. 

An Osher-type and level-set scheme for two-fluid flow computations

Abstract

Similar works

Full text

Available Versions

CWI's Institutional Repository

Pure OAI Repository