A compressible two-fluid model for the finite volume simulation of violent aerated flows. Analytical properties and numerical results
aa r X i v : . [ phy s i c s . c l a ss - ph ] O c t A compressible two-fluid model for the finitevolume simulation of violent aerated flows.Analytical properties and numerical results
Fr´ed´eric Dias a , ∗ a LRC MESO, ENS Cachan, CEA DAM, 61 Av. du President Wilson, F-94230Cachan, FRANCE
Denys Dutykh b b LRC MESO, ENS Cachan, CEA DAM, 61 Av. du President Wilson, F-94230Cachan, FRANCE
Jean-Michel Ghidaglia c c LRC MESO, ENS Cachan, CEA DAM, 61 Av. du President Wilson, F-94230Cachan, FRANCE
Abstract
In the study of ocean wave impact on structures, one often uses Froude scalingsince the dominant force is gravity. However the presence of trapped or entrainedair in the water can significantly modify wave impacts. When air is entrained inwater in the form of small bubbles, the acoustic properties in the water changedramatically and for example the speed of sound in the mixture is much smallerthan in pure water, and even smaller than in pure air. While some work has beendone to study small-amplitude disturbances in such mixtures, little work has beendone on large disturbances in air-water mixtures. We propose a basic two-fluidmodel in which both fluids share the same velocities. It is shown that this modelcan successfully mimic water wave impacts on coastal structures. Even though thisis a model without interface, waves can occur. Their dispersion relation is discussedand the formal limit of pure phases (interfacial waves) is considered. The governingequations are discretized by a second-order finite volume method. Numerical resultsare presented. It is shown that this basic model can be used to study violent aeratedflows, especially by providing fast qualitative estimates.
Key words: wave impact, two-phase flow, compressible flow, free-surface flow,finite volumes
Preprint submitted to Elsevier 21 November 2018
Introduction
One of the challenges in Computational Fluid Dynamics (CFD) is to deter-mine efforts exerted by waves on structures, especially coastal structures. Theflows associated with wave impact can be quite complicated. In particular,wave breaking can lead to flows that cannot be described by models like e.g. the free-surface Euler or Navier–Stokes equations. In a free-surface model, theboundary between the gas (air) and the liquid (water) is a surface. The liq-uid flow is assumed to be incompressible, while the gas is represented by amedium, above the liquid, in which the pressure is constant (the atmosphericpressure in general). Such a description is known to be valid for calculatingthe propagation in the open sea of waves with moderate amplitude, which donot break. Clearly it is not satisfactory when waves either break or hit coastalstructures like offshore platforms, jetties, piers, breakwaters, etc.Our goal here is to investigate a relatively simple two-fluid model that canhandle breaking waves. It belongs to the family of averaged models, in thesense that even though the two fluids under consideration are not miscible,there exists a length scale ǫ such that each averaging volume (of size ǫ ) con-tains representative samples of each of the fluids. Once the averaging process isperformed, it is assumed that the two fluids share, locally, the same pressure,temperature and velocity. Such models are called homogeneous models in theliterature. They can be seen as limiting cases of more general two-fluid modelswhere the fluids could have different temperatures and velocities [9]. Let usexplain why it can be assumed here that both fluids share the same temper-atures and velocities. There are relaxation mechanisms that indeed tend tolocally equalize these two quantities. Concerning temperatures, these are dif-fusion processes and provided no phenomenon is about to produce very stronggradients of temperature between the two fluids like e.g a nuclear reaction inone of the two fluids, one can assume that the time scale on which diffusionacts is much smaller than the time scale on which the flow is averaged. Sim-ilarly, concerning the velocities, drag forces tend to locally equalize the twovelocities. Hence for flows in which the mean convection velocity is moder-ate (a scale of time is built with the mean convection velocity and a typicallength scale) we are in the case where this time scale is much smaller than thetime scale on which velocities are equalized through drag forces. Hence, in thepresent model, the partial differential equations, which express conservationof mass (1 per fluid), balance of momentum and total energy, read as follows: ∗ Corresponding author.
Email addresses:
[email protected] (Fr´ed´eric Dias),
[email protected] (Denys Dutykh), [email protected] (Jean-Michel Ghidaglia). α + ρ + ) t + ∇ · ( α + ρ + ~u ) = 0 , (1)( α − ρ − ) t + ∇ · ( α − ρ − ~u ) = 0 , (2)( ρ~u ) t + ∇ · ( ρ~u ⊗ ~u + p I ) = ρ~g, (3)( ρE ) t + ∇ · ( ρH~u ) = ρ~g · ~u, (4)where the superscripts ± are used to denote liquid and gas respectively. Hence α + and α − denote the volume fraction of liquid and gas, respectively, andsatisfy the condition α + + α − = 1. We denote by ρ ± , ~u , p , e respectively thedensity of each phase, the velocity, the pressure, the specific internal energy, ~g is the acceleration due to gravity, ρ := α + ρ + + α − ρ − is the total density, E = e + | ~u | is the specific total energy, H := E + p/ρ is the specific totalenthalpy. In order to close the system, we assume that the pressure p is givenas a function of three parameters, namely α ≡ α + − α − , ρ and e : p = P ( α, ρ, e ) . (5)We shall discuss in Section 2 how such a function P is determined once thetwo equations of state p = P ± ( ρ ± , e ± ) are known. Equations (1)–(5) form aclosed system that we shall use in order to simulate aerated flows.The main purpose of this paper is to promote a general point of view, whichmay be useful for various applications in ocean, offshore, coastal and arcticengineering. One can say that the late Howell Peregrine was the first to makeuse of this approach. The influence of the presence of air in wave impacts is adifficult topic. While it is usually thought that the presence of air softens theimpact pressures, recent results show that the cushioning effect due to aerationvia the increased compressibility of the air-water mixture is not necessarily adominant effect [3]. First of all, air may become trapped or entrained in thewater in different ways, for example as a single bubble trapped against a wall,or as a column or cloud of small bubbles. In addition, it is not clear whichquantity is the most appropriate to measure impacts. For example some re-searchers pay more attention to the pressure impulse than to pressure peaks.The pressure impulse is defined as the integral of pressure over the short du-ration of impact. A long time ago, Bagnold [1] noticed that the maximumpressure and impact duration differed from one identical wave impact to thenext, even in carefully controlled laboratory experiments, while the pressureimpulse appears to be more repeatable. For sure, the simple one-fluid modelswhich are commonly used for examining the peak impacts are no longer ap-propriate in the presence of air. There are few studies dealing with two-fluidmodels. An exception is the work by Peregrine and his collaborators. Woodet al. [13] used the pressure impulse approach to model a trapped air pocket.Peregrine & Thiais [10] examined the effect of entrained air on a particularkind of violent water wave impact by considering a filling flow. Bullock et al.[4] found pressure reductions when comparing wave impact between fresh andsalt water where, due to the different properties of the bubbles in the two3uids, the aeration levels being much higher in salt water than in fresh water.Bredmose [2] recently performed numerical experiments on a two-fluid systemwhich has similarities with the one we will use below.The novelty of the present paper is not the finite volume method used belowbut rather the modelling of two-fluid flows. Since the model described belowinvolves neither the tracking nor the capture of a free surface, its integrationis much less costly from the computational point of view. We have chosen toreport here on the stiffest case. Should the viscosity effects become important,they can be taken into account via e.g. a fractional step method. In fact, whenviscous effects are important, the flow is easier to capture from the numericalpoint of view.The paper is organized as follows. Section 2 provides an analytical study ofthe model. In Section 3 we show that in the limit where the function α isidentically equal to − It is shown in this section how to determine the function P ( α, ρ, e ) in Eq. (5)once the two equations of state p = P ± ( ρ ± , e ± ) are known. We call Eq. (5) anextended EOS, since P ( − , ρ, e ) = P − ( ρ, e ) and P (1 , ρ, e ) = P + ( ρ, e ), where p ± = P ± ( ρ ± , e ± ) , T ± = T ± ( ρ ± , e ± ) , (6)are the EOS of each fluid. We will use the following prototypical example inthis paper. Assume that the fluid denoted by the superscript − is an ideal gas: p − = ( γ − − ρ − e − , e − = C − V T − , (7)while the fluid denoted by the superscript + obeys to the stiffened gas law(Tait’s law) [8]: p + + π + = ( γ + − ρ + e + , e + = C + V T + + π + γ + ρ + , (8)4here γ ± , C ± V , and π + are constants. For example, pure water is well describedin the vicinity of the normal conditions by taking γ + = 7 and π + = 2 . × Pa.Let us now return to the general case. In order to find the function P , thereare three given quantities: α ∈ [ − ,
1] , ρ > e > ρ ± , e ± the following system of four nonlinear equations:(1 + α ) ρ + + (1 − α ) ρ − = 2 ρ , (9)(1 + α ) ρ + e + + (1 − α ) ρ − e − = 2 ρ e , (10) P + ( ρ + , e + ) − P − ( ρ − , e − ) = 0 , (11) T + ( ρ + , e + ) − T − ( ρ − , e − ) = 0 . (12)For given values of the pressure p > T >
0, we denoteby R ± ( p, T ) and E ± ( p, T ) the solutions ( ρ ± , e ± ) to: P ± ( ρ ± , e ± ) = p , T ± ( ρ ± , e ± ) = T , (13)and then: ρ = 1 + α R + ( p, T ) + 1 − α R − ( p, T ) , (14) ρ e = 1 + α R + ( p, T ) E + ( p, T ) + 1 − α R − ( p, T ) E − ( p, T ) . (15)Finally the inversion of this system of equations leads to p = P ( α, ρ, e ) and T = T ( α, ρ, e ).Concerning the prototypical case, the following generalization of (7) is consid-ered: p − + π − = ( γ − − ρ − e − , e − = C − V T − + π − γ − ρ − . (16)Introducing γ and π defined by2 γ ( α ) − αγ + − − αγ − − , (17)2 π ( α ) γ ( α ) − αγ + − π + + 1 − αγ − − π − , (18)Eq. (14) and (15) then lead to P ( α, ρ, e ) = ( γ ( α ) − ρ e − π ( α ) , (19) T ( α, ρ, e ) = ρ e − ( λ + ( α ) π + + λ − ( α ) π − ) ρ C V ( α ) , (20)5here αC + V ( γ + −
1) + 1 − αC − V ( γ − − ! C V ( α ) = 1 + αγ + − − αγ − − , (21) λ ± ( α ) ≡ ± α γ ± − − C V ( α ) γ ± C ± V ! . (22)One can easily check that one recovers the equations of state for each fluid inthe limits α → ± In this section, we assume that the system of equations is solved in R , havingin mind the numerical computations performed below. However the extensionto 3D is trivial. The system (1)–(4) can be written as ∂ w ∂t + ∇ · F ( w ) = S ( w ) , (23)where w = ( w i ) i =1 := ( α + ρ + , α − ρ − , ρu , ρu , ρE ) , (24)and, for every ~n ∈ R , F ( w ) · ~n = ( α + ρ + ~u · ~n, α − ρ − ~u · ~n, ρ~u · ~nu + pn , ρ~u · ~nu + pn , ρH~u · ~n ) , (25) S ( w ) = (0 , , ρg , ρg , ρ~g · ~u ) . (26)The Jacobian matrix A ( w ) · ~n is defined by A ( w ) · ~n = ∂ ( F ( w ) · ~n ) ∂ w . (27)In order to compute A ( w ) · ~n , one writes Eq. (25) for F ( w ) · ~n in terms of w and p : F ( w ) · ~n = (cid:18) w w n + w n w + w , w w n + w n w + w , w w n + w n w + w + pn ,w w n + w n w + w + pn , ( w + p ) w n + w n w + w (cid:19) . (28)The Jacobian matrix (27) then has the following expression:6 ( w ) · ~n = u n α − ρ − ρ − u n α + ρ + ρ α + ρ + ρ n α + ρ + ρ n − u n α − ρ − ρ u n α + ρ + ρ α − ρ − ρ n α − ρ − ρ n − u u n + ∂p∂w n − u u n + ∂p∂w n u n + u n + ∂p∂w n u n + ∂p∂w n ∂p∂w n − u u n + ∂p∂w n − u u n + ∂p∂w n u n + ∂p∂w n u n + u n + ∂p∂w n ∂p∂w n u n (cid:16) ∂p∂w − H (cid:17) u n (cid:16) ∂p∂w − H (cid:17) u n ∂p∂w + Hn u n ∂p∂w + Hn u n (cid:16) ∂p∂w (cid:17) , where u n = ~u · ~n. Let us now compute the five derivatives ∂p/∂w i . A systematic way of doing itis to introduce a set of five independent physical variables and here we shalltake: ϕ = α, ϕ = p, ϕ = T, ϕ = u , ϕ = u . (29)The expressions of the w ′ i s in terms of the ϕ ′ j s are algebraic and explicit.Hence the Jacobian matrix ∂w i /∂ϕ j can be easily computed. Since ∂ϕ j /∂w i is its inverse matrix, one finds easily with the help of a computer algebraprogram that ∂p∂w = Γ −
12 ( u + u ) + α − ρ − χ − , (30) ∂p∂w = Γ −
12 ( u + u ) + α + ρ + χ + , (31) ∂p∂w = − (Γ − u , ∂p∂w = − (Γ − u , ∂p∂w = Γ − , (32)where χ ∓ = 1 ρ ± ( c ∓ s ) γ ∓ − − ρ ∓ ( c ± s ) γ ± − , χ + + χ − = 0 , (33)( c ± s ) ≡ C ± V γ ± ( γ ± − T = γ ± p + π ± ρ ± , (34)Γ − ≡ ( γ ( α ) − ρc s γ ( α ) p + π ( α ) . (35)In Eq. (35), we have introduced the speed of sound of the mixture c s , definedby 1 ρc s = (1 + α ) γ + ρ + ( c + s ) + (1 − α ) γ − ρ − ( c − s ) − ρa , (36)with ρa ≡ (1 + α ) ρ + ( c + s ) γ + −
1) + (1 − α ) ρ − ( c − s ) γ − − . (37)Then one can show that the Jacobian matrix A ( w ) · ~n has three distincteigenvalues: λ = u n − c s , λ , , = u n , λ = u n + c s , (38)7nd is diagonalizable on R . The expression of a set of eigenvectors can beobtained by using a computer algebra program. Remark 1 If π + = 0 and π − = 0 , then c s = γ ( α ) pρ and a = c s γ ( α ) − . Remark 2
The left hand side of (36) is positive since ρa is bounded frombelow by (1+ α ) ρ + ( c + s ) γ + + (1 − α ) ρ − ( c − s ) γ − The system of conservation laws (1)–(4) can be transformed into a set ofevolution equations for the physical variables. Let us introduce the entropyfunction s ( ~x, t ) defined by (compare with Eq. (10))2 ρ s = (1 + α ) ρ + s + + (1 − α ) ρ − s − . Proposition 1
Continuous solutions to (1)–(4) satisfy ~u t + ~u · ∇ ~u + 1 ρ ∇ p = ~g , (39) p t + ~u · ∇ p + ρc s ∇ · ~u = 0 , (40) α t + ~u · ∇ α + (1 − α ) δ ∇ · ~u = 0 , (41) s t + ~u · ∇ s = 0 , (42) where c s is given by (36)-(37) and δ is given by δ ≡ ρc s ( γ − π + − γ + π − ) ρ + ρ − ( c + s ) ( c − s ) . (43) Remark 3
For pure fluids ( α = ± ), Eq. (41) is no longer relevant and δ is not needed. One can check that the speed of sound c s is then equal to theexpected speed of sound ( c + s or c − s ) for pure fluids. The balance of entropy (42) comes from the balance( ρs ) t + ∇ · ( ρs~u ) = 0 . (44)Adding together Eqs (1) and (2) leads to ρ t + ∇ · ( ρ~u ) = 0 . (45)Combining Eqs (44) and (45) leads to Eq. (42).8 emark 4 Subtracting Eq. (1) from (2) leads to ( ρχ ) t + ∇ · ( ρχ~u ) = 0 , with χ = α + ρ + − α − ρ − ρ . (46) In the case of smooth solutions, we obtain that χ t + ~u · ∇ χ = 0 , which is an alternative to Eq. (41). From now on, we denote the set of equations (1)–(4) by (E). In order to studysmall perturbations around basic smooth and stationary solutions, it is moreconvenient to use the general set of equations (E) rewritten in the physicalvariables α , ~u , p and s : see Eq. (39)–(42).The two-fluid model (E) describes the evolution of mixtures. It can be usedfor example to study waves along a diffuse interface between a gas and aliquid. In order to find the dispersion relation for such waves, one first looksfor rest states. There is an infinity of such rest states. Then, the governingequations are linearized around a special class of rest states. Here the situationis considerably complicated by the fact that the stationary solutions are notuniform in space. Thus, we come up with a linear system of partial differentialequations which have non-constant coefficients.The steady state will be denoted by α ± , ρ ± , p , ~u and s . The special class of so-lutions we are looking for are motionless, uniform in the horizontal coordinatesand continuously stratified in the vertical direction: α ± = α ± ( z ) , ρ ± = ρ ± ( z ) , p = p ( z ) , ~u = ~ , s = s ( z ) . In this case, Eqs (39)–(42) become dpdz = − ρ ( z ) g . (47)In the case where one makes the assumption that the mixture density is con-stant, ρ ( z ) = 12 (1 + α ) ρ + + 12 (1 − α ) ρ − ≡ ρ , ρ is some positive constant, the solution to (47) is p ( z ) = p − ρ gz, (48)where p is a constant.Equation (48) combined with (19) leads to p − ρ gz = ( γ ( α ) − ρ e − π ( α ) , (49)and there are infinitely many solutions. Since we have imposed ρ ( z ) = ρ ,it makes sense to also impose a temperature which does not depend on z : T ( z ) = T . Otherwise, there would be some motion due to convection ( ~u = 0).Thus, Eq. (20) leads to the equation e = C V ( α ) T + λ + ( α ) π + + λ − ( α ) π − ρ , (50)which can be used to obtain an equation for α when combined with Eq. (49). In this section we linearize the governing equations around the basic steadystate derived in the previous section. We write down the following perturbationof the stationary solution: α = α ( z )+2 β + . . . , p = p ( z )+ q + . . . , ~u = ~ ~v + . . . , s = s ( z )+ σ + . . . , where the vector ~v has the three components ( v , v , v ).From Eqs (40)–(42), it is straightforward to check that q , β and σ satisfy theequations ∂q∂t + ~v · ∇ p ( z ) + ρ c s ∇ · ~v = 0 , (51) ∂β∂t + ~v · ∇ α ( z ) + 12 (1 − α ) δ ∇ · ~v = 0 , (52) ∂σ∂t + ~v · ∇ s ( z ) = 0 . (53)In order to obtain the velocity perturbation equation, one needs to computecarefully the density perturbation. Using that dρ ± dz = γ ± ( c ± s ) dpdz (recall that dTdz = 0) , = 12 (1 + α ( z ) + 2 β + . . . ) R + ( p ( z ) + q + . . . , T )+ 12 (1 − α ( z ) − β + . . . ) R − ( p ( z ) + q + . . . , T )= 12 (1 + α ( z ) + 2 β + . . . ) ρ + + γ + q ( c + s ) + . . . ! + 12 (1 − α ( z ) − β + . . . ) ρ − + γ − q ( c − s ) + . . . ! = ρ + β ( ρ + − ρ − ) + 12 (1 + α ) γ + ( c + s ) + (1 − α ) γ − ( c − s ) ! q + O ( || π || + || β || ) . Having computed ρ , it is easy to write the equation for ~v from Eq. (39): ∂~v∂t + 1 ρ ∇ q = " ρ + − ρ − ρ β + 12 ρ (1 + α ) γ + ( c + s ) + (1 − α ) γ − ( c − s ) ! q ~g . (54) From now on, the following notation will be used: a ( z ) := ρ c s , b ( z ) := 12 (1 − α ) δ ,c ( z ) := ρ + − ρ − ρ , d ( z ) := 12 ρ (1 + α ) γ + ( c + s ) + (1 − α ) γ − ( c − s ) ! . (55) Consider first the simple case where the acceleration due to gravity is absent: ~g = 0. It represents a major simplification since in this case we recover from(51)–(54) a system of partial differential equations with constant coefficients: ∂q∂t + ρ c s ∇ · ~v = 0 ,∂β∂t + 12 (1 − α ) δ ∇ · ~v = 0 ,∂σ∂t = 0 ,∂~v∂t + 1 ρ ∇ q = 0 . ∂ w ∂t + X i =1 B i ∂ w ∂x i = 0 , (56)where w = ( q, β, σ, ~v ). We look for plane wave solutions: w ( ~x, t ) = w e i ( ~k · ~x − ωt ) , ~x = ( x , x , x ) with x = z, ~k = ( k , k , k ) . Substituting this ansatz into equation (56) yields B ( ~k ) w = ω w , with B ( ~k ) = X i =1 k i B i or B ( ~k ) = a k a k a k b k b k b k k ρ k ρ k ρ , a = ρ c s , b = 12 (1 − α ) δ . It means that ω is an eigenvalue and w the corresponding eigenvector of thematrix B k .The dispersion relation of the system (56) is given by its characteristic poly-nomial ω ω − a ρ | ~k | ! = 0 or ω (cid:16) ω − c s | ~k | (cid:17) = 0 . Note that here we have not considered any boundary conditions and that thevertical direction does not play any particular role. This is why we have beenlooking for perturbations which are periodic in all three directions. In factthere is no dispersion in the acoustic waves we have found.
Remark 5
The computations performed above are in full agreement with thecomputations (23)–(27). The matrix B is similar to A (0) . Note that we ob-tained here 6 eigenvalues as opposed to 5 in Section 2.2. This is only due tothe fact that we performed the computations in 3D in this section as opposedto 2D in Section 2.2. .3.2 General case with gravity Let us now consider the general situation where g is different from zero. Wedrop the equation for σ since (51)–(53) is a strictly “triangular” linear systemof PDEs. As above we look for periodic solutions of the form qβ~v = ˆ q ˆ β ˆ v ( z ) e i ( ~k · ~x − ωt ) . (57)In this case ~x and ~k have only horizontal components: ~x = ( x , x ) , ~k = ( k , k ) , and ˆ v = (ˆ v , ˆ v , ˆ w ). The main difference with the previous case is that theamplitude now depends on the vertical coordinate z . It makes the analysismore complicated.Substituting the expression (57) into Eqs (51)–(54) except for (53) yields thefollowing system of ordinary differential equations: − iω ˆ q − ρ g ˆ w + a ( z ) (cid:18) i~k · ˆ v , + d ˆ wdz (cid:19) = 0 , − iω ˆ β + dαdz ˆ w + b ( z ) (cid:18) i~k · ˆ v , + d ˆ wdz (cid:19) = 0 , − iω ˆ v , + i~k ˆ qρ = ~ , − iω ˆ w + 1 ρ d ˆ qdz + c ( z ) g ˆ β + d ( z ) g ˆ q = 0 . The third equation yields the horizontal divergence of ~v in Fourier space interms of the pressure perturbation ˆ q : i~k · ˆ v , = i | ~k | ω ˆ qρ . Another algebraic identity can be obtained if we multiply the first equationby b ( z ), the second one by a ( z ) and subtract them: − iω ˆ βa ( z ) + (cid:18) dαdz a ( z ) + ρ gb ( z ) (cid:19) ˆ w + iω ˆ qb ( z ) = 0 . β : ˆ β = b ( z ) a ( z ) ˆ q + (cid:18) dαdz + ρ g b ( z ) a ( z ) (cid:19) ˆ wiω . Thus, the system governing the behaviour of the perturbations ˆ w and ˆ q isgiven by the two equations d ˆ wdz + i (cid:18) | ~k | ρ ω − ωa ( z ) (cid:19) ˆ q − ρ ga ( z ) ˆ w = 0 , (58)1 ρ d ˆ qdz + g (cid:18) d ( z ) + c ( z ) b ( z ) a ( z ) (cid:19) ˆ q + " ω + gc ( z ) (cid:18) dαdz + ρ g b ( z ) a ( z ) (cid:19) ˆ wiω = 0 . (59)Note that this analysis is consistent with the previous one without gravity.Indeed, if one takes g = 0 in (58)-(59) and assumes a periodic dependence in z with wavenumber k , one recovers the previous dispersion relation.In order to find the dispersion relation in the general case, these equationsmust completed by boundary conditions, for exampleˆ w bottom = 0; ˆ w top = 0 , if the flow occurs between two solid walls.One expects to obtain a dispersion relation as a solvability condition for thesecond-order boundary value problem. Unfortunately, at this stage, we cannotgo much further, except for the particular case where gz/ ( c ± s ) ≪
1. Then onecan perform an asymptotic expansion in this small parameter. As a result thecoefficients are polynomials in z . Restricting to the leading order terms yieldsa system of two ordinary differential equations with coefficients that are affinein z . The solution can be obtained in terms of Airy’s wave functions Ai andBi. We show now that the two-fluid model degenerates into the classical water-wave equations in the limit of an interface separating two pure fluids. Considerthe case where α is either 1 or −
1. More precisely let α := 1 − H ( z − η ( ~x, t )) , ~x = ( x , x ) (60)where H is the Heaviside step function. Physically this substitution meansthat we consider two pure fluids separated by an interface. It follows that α + α − = 0 , − α = 0 . η t + ~u h · ∇ h η = w , where ~u h = ( u , u ) and ∇ h = ( ∂ x , ∂ x ).This equation simply states that there is no mass flux across the interface.Incidentally this is no longer true in the case of shock waves. Integratingthe conservation of momentum equation (3) inside a volume moving withthe flow and enclosing the interface and using the fact there is no mass fluxacross the interface simply leads to the fact that there is no pressure jumpacross the interface. In other words, the pressure is continuous across theinterface. Integrating the entropy equation inside the same volume enclosingthe interface and using the fact there is no mass flux across the interface doesnot lead to any new information.One can now write Eqs (2)–(4) in each fluid, either in the conservative form( ρ ± ) t + ∇ · ( ρ ± ~u ± ) = 0 , (61)( ρ ± ~u ± ) t + ∇ · ( ρ ± ~u ± ⊗ ~u ± ) + ∇ p ± = ρ ± ~g , (62)( ρ ± s ± ) t + ∇ · ( ρ ± s ± ~u ± ) = 0 , (63)or in the more classical form ρ ± t + ( ~u ± · ∇ ) ρ ± + ρ ± ∇ · ~u ± = 0 , (64) ~u ± t + ( ~u ± · ∇ ) ~u ± + ∇ p ± ρ ± = ~g , (65) s ± t + ~u ± · ∇ s ± = 0 . (66)In these two systems, the superscripts + and − are used for the heavy fluid(below the interface) and the light fluid (above the interface) respectively.The system of equations we derived is nothing else than the system of adiscontinuous two-fluid system with an interface located at z = η ( ~x, t ). Alongthe interface, one has the kinematic and dynamic boundary conditions η t + ~u ± h · ∇ h η = w ± (67) p − = p + (68)This simple computation shows an interesting property of our model: it au-tomatically degenerates into a discontinuous two-fluid system where two purecompressible phases are separated by an interface. In Appendix B, we derive15 kH/2 π ω / π n=4n=3n=2n=1 Fig. 1. Dispersion relation (71) for the acoustic mode. the dispersion relation for this limit. We choose the approximate rest state η = 0 , ρ ± = ρ ± , ~u ± = ~ , p ± = p − ρ ± gz, s ± = s ± , and assume that the fluid domain is bounded by two horizontal walls locatedat z = ∓ α ± D ( D is the total depth of the domain). Let θ = ρ − /ρ +0 andintroduce S + ω = vuut − ω | ~k | ( c + s ) , T − ω = vuut ω | ~k | ( c − s ) − , S − ω = vuut − ω | ~k | ( c − s ) . (69)There is a distinction between three cases. When c − s | ~k | < ω < c + s | ~k | , one findsthe dispersion relation (B.15), which is reproduced here: ω g | ~k | (cid:16) T − ω tan( T − ω α − | ~k | D ) − θS + ω tanh( S + ω α +0 | ~k | D ) (cid:17) == (1 − θ ) S + ω T − ω tan( T − ω α − | ~k | D ) tanh( S + ω α +0 | ~k | D ) . (70)Neglecting the effects due to gravity, the dispersion relation (70) becomes T − ω tan( T − ω α − | ~k | D ) = θS + ω tanh( S + ω α +0 | ~k | D ) . (71)There is an infinite number of solutions because of the presence of the tan-gent term tan( T − ω α − | ~k | D ). Since θ is small in coastal engineering applica-tions, we expect T − ω tan( T − ω α − | ~k | D ) to be small. Then either ( T − ω ) is small,or T − ω α − | ~k | D ≈ nπ , n ∈ Z .After some simple calculations, one obtains for ( T − ω ) small ω ≈ c − s | ~k | , and for T − ω α − | ~k | D ≈ nπω n = c − s | ~k | vuut n π ( α − ) | ~k | D , with lim | ~k | D → ω n = nπc − s α − D . ω < min( c − s | ~k | , c + s | ~k | ), one finds the dispersion relation (B.16), which isreproduced here: ω g | ~k | (cid:16) S − ω tanh( S − ω α − | ~k | D ) + θS + ω tanh( S + ω α +0 | ~k | D ) (cid:17) == (1 − θ ) S − ω S + ω tanh( S − ω α − | ~k | D ) tanh( S + ω α +0 | ~k | D ) . (72)In the incompressible limit, the speeds of sound c ± s go to infinity, S ± ω → ω g | ~k | (cid:16) tanh( α − | ~k | D ) + θ tanh( α +0 | ~k | D ) (cid:17) = (1 − θ ) tanh( α − | ~k | D ) tanh( α +0 | ~k | D ) . Finally, when ω > max( c + s | ~k | , c − s | ~k | ), one finds the dispersion relation (B.17),which is reproduced here: ω g | ~k | (cid:16) T − ω tan( T − ω α − | ~k | D ) + θT + ω tan( T + ω α +0 | ~k | D ) (cid:17) == − (1 − θ ) T + ω T − ω tan( T − ω α − | ~k | D ) tan( T + ω α +0 | ~k | D ) . (73)Neglecting the effects due to gravity, the dispersion relation (73) becomes T − ω tan( T − ω α − | ~k | D ) + θT + ω tan( T + ω α +0 | ~k | D ) = 0 . (74)There is an infinite number of solutions because of the presence of the tangentterm. Again, since θ is small in coastal engineering applications, we expect T − ω tan( T − ω α − | ~k | D ) to be small. Then we are back to the first case. In this thirdcase, if we were not prescribing boundary conditions at the bottom and at thetop, we could look for perturbations which are periodic in all three directionswith wavenumber k ± in the z − direction. Equation (B.13) then gives ω = ( c ± s ) ( | ~k | + ( k ± ) ) . In other words, there is a relationship between the two vertical wavenumbers.This relationship is 1 + k +3 | ~k | ! = c − s c + s ! k − | ~k | ! . It is reminiscent of Snell’s law, which describes the relationship between theangles of incidence and refraction, when referring to waves passing through aboundary between two different isotropic media.17
A finite-volume discretization of the model
Here we describe the discretization of the model (1)–(4) by a standard cell-centered finite volume method. The computational domain Ω ⊂ R d is trian-gulated into a set of control volumes: Ω = ∪ K ∈T K . We start by integratingequation (23) on K : ddt Z K w d Ω + X L ∈N ( K ) Z K ∩ L F ( w ) · ~n KL dσ = Z K S ( w ) d Ω , (75)where ~n KL denotes the unit normal vector on K ∩ L pointing into L and N ( K ) = { L ∈ T : area( K ∩ L ) = 0 } . Then, setting w K ( t ) := 1vol( K ) Z K w ( ~x, t ) d Ω , we approximate (75) by d w K dt + X L ∈N ( K ) area( L ∩ K )vol( K ) Φ( w K , w L ; ~n KL ) = S ( w K ) , (76)where the numerical fluxΦ( w K , w L ; ~n KL ) ≈ L ∩ K ) Z K ∩ L F ( w ) · ~n KL dσ , is explicitly computed by the FVCF formula of Ghidaglia et al. [6]:Φ( v , w ; n ) = F ( v ) · ~n + F ( w ) · ~n − sgn( A n ( µ ( v , w ))) F ( w ) · ~n − F ( v ) · ~n , (77)where the Jacobian matrix A n ( µ ) is defined in (27), µ ( v , w ) is an arbitrarymean between v and w and sgn( M ) is the matrix whose eigenvectors are thoseof M but whose eigenvalues are the signs of the eigenvalues of M .So far we have not discussed the case where a control volume K meets theboundary of Ω. Here we shall only consider the case where this boundary is awall and from the numerical point of view, we only need to find the normalflux F · ~n . Since ~u ( ~x, t ) · ~n = 0 for ~x ∈ ∂ Ω , we have( F · ~n ) | ~x ∈ ∂ Ω = (0 , , p b n , p b n , , p b := p | ~x ∈ ∂ Ω , and following Ghidaglia and Pascal [7], we can take p b = p + ρu n c s , where theright-hand side is evaluated in the control volume K . Remark 1
In order to turn (76) into a numerical algorithm, we must atleast perform time discretization and give an expression for µ ( v , w ) . Since his matter is standard, we do not give the details here but instead refer toDutykh [5]. Let us also notice that formula (76) leads to a first-order schemebut in fact we use a MUSCL technique to achieve better accuracy in space [12]. In order to check the accuracy of our second-order scheme, we first solvenumerically the scalar linear advection equation ∂v∂t + ~a · ∇ v = 0 , ~a ∈ R , with smooth initial conditions with compact support in order to reduce theinfluence of boundary conditions. It is obvious that this equation will justtranslate the initial form in the direction ~a . So, we have an analytical solutionwhich can be used to quantify the error of the numerical method. On the otherhand, to measure the convergence rate, we constructed a sequence of refinedmeshes.Fig. 2 shows the error of the numerical method in L ∞ norm as a function of themesh characteristic size. The slope of the curves represents an approximationto the theoretical convergence rate. On this plot, the curve with circles forthe data points corresponds to the first order upwind scheme while the othertwo correspond to the MUSCL scheme with least-squares and Green-Gaussgradient reconstruction procedures respectively. The slope of the curve withcircles is equal approximatively to 0 .
97, which means first-order convergence.The other two curves have almost the same slope equal to 1 .
90, thus indicatinga second-order convergence rate for the MUSCL scheme. Remark that in ourimplementation of the second-order scheme the least-squares reconstructionseems to give slightly more accurate results than the Green-Gauss procedure.The next figure represents the measured CPU time in seconds, again as afunction of the mesh size. Obviously, this kind of data is extremely computerdependent but the qualitative behaviour is the same on all systems. On Fig. 3one can see that the “fastest” curve is that of the first-order upwind scheme.Then we have two almost superimposed curves referring to the second-ordergradient reconstruction on variables. Here again one can notice that the least-squares method is slightly faster than the Green-Gauss procedure. On thisfigure we represented one more curve (the highest one) which corresponds toGreen-Gauss gradient reconstruction on fluxes (it seems to be very natural inthe context of FVCF schemes). The numerical tests show that this method is19 −2 −1 −4 −3 −2 −1 h − average edge length L ∞ e rr o r Convergence rate in L ∞ norm Upwind. Conv. rate = 0.97731LeastSq. Conv. rate = 1.9051GGO on Variables. Conv. rate = 1.9012 Fig. 2. Numerical method error in L ∞ norm. −2 −1 −3 −2 −1 CPU time in secondsh T Upwind 1st orderGreen−GaussLeast SqsGG on fluxes
Fig. 3. CPU time for different finite volume schemes.
Fig. 4. Shock tube of Sod: density plot.
The geometry and initial condition for this test case are shown on Fig. 5. Ini-tially the velocity field is taken to be zero. The values of the other parametersare given in Table A.1. The mesh used in this computation contained about108000 control volumes (in this case they were triangles). The results of thissimulation are presented on Figures 6–11. Fig. 12 shows the maximal pressureon the right wall as a function of time: t max ( x,y ) ∈ × [0 , p ( x, y, t ) . We performed another computation for a mixture with α + = 0 . α − = 0 . α + = 0 . α − = 0 . α + = 0 . α − = 0 .
90 0 . .
65 0 . .
051 10 . ~g Fig. 5. Falling water column test case. Geometry and initial condition.(a) t = 0 .
005 (b) t = 0 . The geometry and initial condition for this test case are shown on Fig. 14.Initially the velocity field is taken to be zero. The values of the other param-eters are given in Table A.1. The mesh used in this computation containedabout 92000 control volumes (again they were triangles). The results of thissimulation are presented in Figures 15–21. In Fig. 22 we plot the maximal22 a) t = 0 . t = 0 . t = 0 .
15 (b) t = 0 . t = 0 . t = 0 . pressure on the bottom as a function of time: t max ( x,y ) ∈ [0 , × p ( x, y, t ) . a) t = 0 . t = 0 . t = 0 . t = 0 . p m a x / p Maximal pressure on the right wall
Fig. 12. Maximal pressure on the right wall as a function of time. Case of a heavygas. p m a x / p Maximal pressure on the right wall
Fig. 13. Maximal pressure on the right wall as a function of time. Case of a lightgas.
It is clear that the pressure exerted on the bottom reaches 2 . p due to thedrop impact at t ≈ . Remark 2
Beginning with Fig. 20 one can see some asymmetry in the solu-tion. It is not expected since the initial condition, computational domain andforcing term are fully symmetric with respect to the line x = 0 . . This discrep-ancy is explained by the use of unstructured meshes in the computation. Thearbitrariness of the orientation of the triangles introduces small perturbationswhich are sufficient to break the symmetry at the discrete level. In this article we have presented a simple mathematical model for simulatingwater wave impacts. Associated to this model, which avoids the costly captureof free surfaces, we have built a numerical solver which is: (i) second-orderaccurate on smooth solutions, (ii) stable even for solutions with very stronggradients (and solutions with shocks) and (iii) locally exactly conservativewith respect to the mass of each fluid, momentum and total energy. This lastproperty, (iii), which is certainly the most desirable from the physical point ofview, is an immediate byproduct of our cell-centered finite volume method.We have shown here the good behavior of this framework on simple test casesand we are presently working on quantitative comparisons in the context of25Sfrag replacements α + = 0 . α − = 0 . α + = 0 . α − = 0 .
10 0 . .
71 1 ~gR = 0 . Fig. 14. Geometry and initial condition for water drop test case.(a) t = 0 .
005 (b) t = 0 . real applications. A Some technical results
The constants C ± V can be calculated after simple algebraic manipulations ofequations (7), (8) and matching with experimental values at normal condi-26 a) t = 0 . t = 0 . t = 0 .
135 (b) t = 0 . t = 0 .
175 (b) t = 0 . tions: C − V ≡ p ( γ − − ρ − T ,C + V ≡ γ + p + π + ( γ + − γ + ρ +0 T . a) t = 0 .
225 (b) t = 0 . t = 0 .
325 (b) t = 0 . t = 0 . t = 0 . For example, for an air/water mixture under normal conditions we have thevalues given in Table A.1.The sound velocities in each phase are given by the following formulas:( c − s ) = γ − p − ρ − , ( c + s ) = γ + p + + π + ρ + . (A.1)28 p m a x / p Maximal pressure on the bottom
Fig. 22. Water drop test case. Maximum bottom pressure as a function of time. parameter value p P aρ +0 kg/m ρ − . kg/m T Kγ − . γ + π + . × P aC + V . Jkg · K C − V . Jkg · K g ms Table A.1Values of the parameters for an air/water mixture under normal conditions. Therather high value of the acceleration due to gravity does not correspond to anyphysical situation. Nevertheless, this value was chosen in order to accelerate thedynamic processes in the test cases.
B Dispersion relation in the pure fluid limit
Let us provide the dispersion relation for waves propagating along the inter-face in the limit of two superposed compressible heavy fluids. Consider thelinearization of the equations around the equilibrium state29 ( ~x, t ) = 0 ,ρ ± ( ~x, z, t ) = ρ ± exp − gz ( c ± s ) ! ,~u ± ( ~x, z, t ) = ~ ,p ± ( ~x, z, t ) = ( c ± s ) ρ ± exp − gz ( c ± s ) ! − ( c ± s ) ρ ± + p ,s ± = s ± If we assume that gz/ ( c ± s ) is small, an approximation to the equilibrium stateis given by η ( ~x, t ) = 0 , (B.1) ρ ± ( ~x, z, t ) = ρ ± , (B.2) ~u ± ( ~x, z, t ) = ~ , (B.3) p ± ( ~x, z, t ) = p − ρ ± gz, (B.4) s ± = s ± (B.5)We write the following perturbations of the equilibrium state: η = η + ζ + . . . , ρ ± = ρ ± + ̺ ± + . . . , ~u ± = ~ ~v ± + . . . , p ± = p ± + q ± + . . . , with in addition s ± = s ± + σ ± .The linearized equations for both fluids (64)–(66) read ∂̺ ± ∂t + ρ ± ∇ h · ~v ± h + ρ ± ∂w ± ∂z = 0 , (B.6) ∂~v ± h ∂t + 1 ρ ± ∇ h q ± = 0 , (B.7) ∂w ± ∂t + 1 ρ ± ∂q ± ∂z = 0 , (B.8) ∂q ± ∂t = (cid:16) c ± s (cid:17) ∂̺ ± ∂t . (B.9)The kinematic and dynamic boundary conditions along the interface are ∂ζ∂t = w ± ( ~x, , t ) , (B.10) q + ( ~x, , t ) − ρ +0 gζ = q − ( ~x, , t ) − ρ − gζ . (B.11)30e perform a classical perturbation analysis: we look for solutions in the form ̺ ± ~v ± h w ± q ± = ˆ ̺ ± ˆ ~v ± h ˆ w ± ˆ q ± ( z ) e i ( ~k · ~x − ωt ) , (B.12)that is periodic perturbations with wave number ~k and angular frequency ω .One eventually obtains the following second-order ODE for ˆ q ± ( z ): d ˆ q ± dz − | ~k | − ω ( c ± s ) ! ˆ q ± = 0 . (B.13)Having q ± , one can find ~v ± , ̺ ± and ζ . Assume a geometry of total depth D bounded above and below by rigid walls located at z = α − D and z = − α +0 D respectively. The boundary conditions along the horizontal walls are w ± ( ~x, ∓ α ± D, t ) = 0 . (B.14)Satisfying the kinematic and dynamic boundary conditions along the interfaceprovides a solvability condition. In other words, the solution to the linearizedproblem provides the dispersion relation ω ( ~k ), which shows the presence oftwo kinds of modes: the gravity modes and the acoustic modes.The ODE for ˆ q ± shows that the sign of ω − | ~k | ( c ± s ) plays an importantrole. There are three cases: both signs are negative, one is positive and one isnegative, both signs are positive. B.1 Case where c − s | ~k | < ω < c + s | ~k | The assumption is based on the fact that the speed of sound c + s is much higherthan the speed of sound c − s . The solution of the linearized problem is31 − = A cos h T − ω | ~k | ( z − α − D ) i e i ( ~k · ~x − ωt ) ,w − = iA | ~k | ωρ − T − ω sin h T − ω | ~k | ( z − α − D ) i e i ( ~k · ~x − ωt ) ,~v − = A ~kωρ − cos h T − ω | ~k | ( z − α − D ) i e i ( ~k · ~x − ωt ) ,̺ − = A c − s ) cos h T − ω | ~k | ( z − α − D ) i e i ( ~k · ~x − ωt ) ,ζ = A | ~k | ω ρ − T − ω sin( T − ω | ~k | α − D )e i ( ~k · ~x − ωt ) , and q + = B cosh h S + ω | ~k | ( z + α +0 D ) i e i ( ~k · ~x − ωt ) ,w + = − iB | ~k | ωρ +0 S + ω sinh h S + ω | ~k | ( z + α +0 D ) i e i ( ~k · ~x − ωt ) ,~v + = B ~kωρ +0 cosh h S + ω | ~k | ( z + α +0 D ) i e i ( ~k · ~x − ωt ) ,̺ + = B c + s ) cosh h S + ω | ~k | ( z + α +0 D ) i e i ( ~k · ~x − ωt ) ,ζ = B | ~k | ω ρ +0 S + ω sinh( S + ω | ~k | α +0 D )e i ( ~k · ~x − ωt ) . Writing that both expressions for ζ are equal and satisfying the dynamic condi-tion along the interface yields a system of two homogeneous linear equationsin A and B . The dispersion relation for the acoustic modes is obtained bysetting the determinant equal to 0: ω g | ~k | (cid:16) T − ω tan( T − ω α − | ~k | D ) − θS + ω tanh( S + ω α +0 | ~k | D ) (cid:17) == (1 − θ ) S + ω T − ω tan( T − ω α − | ~k | D ) tanh( S + ω α +0 | ~k | D ) , (B.15)where θ := ρ − /ρ +0 . B.2 Case where ω < min( c − s | ~k | , c + s | ~k | )Assume now that ω < c − s | ~k | . The solution of the linearized problem is un-changed for the liquid. For the gas it becomes32 − = A cosh h S − ω | ~k | ( z − α − D ) i ,w − = − iA | ~k | ωρ − S ω sinh h S − ω | ~k | ( z − α − D ) i ,~v − = A ~kωρ − cosh h S − ω | ~k | ( z − α − D ) i ,̺ − = A c cosh h S − ω | ~k | ( z − α − D ) i ,ζ = − A | ~k | ω ρ − S − ω sinh( S − ω | ~k | α − D ) . Writing that both expressions for ζ are equal and satisfying the dynamic con-dition along the interface yields a system of two homogeneous linear equationsin A and B . The dispersion relation for the gravity modes is obtained by set-ting the determinant equal to 0: ω g | ~k | (cid:16) S − ω tanh( S − ω α − | ~k | D ) + θS + ω tanh( S + ω α +0 | ~k | D ) (cid:17) == (1 − θ ) S − ω S + ω tanh( S − ω α − | ~k | D ) tanh( S + ω α +0 | ~k | D ) . (B.16) B.3 Case where ω > max( c + s | ~k | , c − s | ~k | )Assume now that ω > c + s | ~k | . The solution of the linearized problem for theliquid becomes q + = B cos h T + ω | ~k | ( z + α +0 D ) i e i ( ~k · ~x − ωt ) ,w + = iB | ~k | ωρ +0 T + ω sin h T + ω | ~k | ( z + α +0 D ) i e i ( ~k · ~x − ωt ) ,~v + = B ~kωρ +0 cos h T + ω | ~k | ( z + α +0 D ) i e i ( ~k · ~x − ωt ) ,̺ + = B c + s ) cos h T + ω | ~k | ( z + α +0 D ) i e i ( ~k · ~x − ωt ) ,ζ = − B | ~k | ω ρ +0 T + ω sin( T + ω | ~k | α +0 D )e i ( ~k · ~x − ωt ) . Writing that both expressions for ζ are equal and satisfying the dynamic condi-tion along the interface yields a system of two homogeneous linear equationsin A and B . The dispersion relation for the acoustic modes is obtained by33etting the determinant equal to 0: ω g | ~k | (cid:16) T − ω tan( T − ω α − | ~k | D ) + θT + ω tan( T + ω α +0 | ~k | D ) (cid:17) == − (1 − θ ) T + ω T − ω tan( T − ω α − | ~k | D ) tan( T + ω α +0 | ~k | D ) . (B.17) C Isentropic flows
Let us consider the isentropic version of the system of equations (1)–(4). Itreads ( α + ρ + ) t + ∇ · ( α + ρ + ~u ) = 0 , (C.1)( α − ρ − ) t + ∇ · ( α − ρ − ~u ) = 0 , (C.2)( ρ~u ) t + ∇ · ( ρ~u ⊗ ~u + p I ) = ρ~g, (C.3)with the equation of state p = P I ( α, ρ ) , (C.4)where the subscript I stands for isentropic.One can determine P I as follows. First consider the two equations(1 + α ) ρ + + (1 − α ) ρ − = 2 ρ , (C.5) P + I ( ρ + ) − P − I ( ρ − ) = 0 . (C.6)Given p >
0, we denote by R ± I ( p ) the solutions ρ ± to P ± I ( ρ ± ) = p . (C.7)It follows that ρ = 1 + α R + I ( p ) + 1 − α R − I ( p ) . (C.8)The inversion of this equation leads to p = P I ( α, ρ ), which is the equation ofstate that appears in (C.4).In order to see if one can go further analytically, let us consider the particularcase of stiffened gases. The equations of state are p ± + π ± = ( γ ± − ρ ± e ± , (C.9)and e ± = C ± V T ± + π ± γ ± ρ ± . roposition 2 Possible entropies s ± for stiffened gases are given by s ± = C ± V log p ± + π ± /γ ± ( ρ ± ) γ ± ! . (C.10)This expression can be easily obtained by integrating the well-known differ-ential relation T d s = d e + p d ρ ! . Remark 6
Other possible entropies s ± for stiffened gases, which differ fromthe previous ones by a constant, are given by s ± = C ± V log e ± − π ± /γ ± ρ ± ( ρ ± ) γ ± − ! . (C.11)Thus, saying that s ± is constant boils down to saying that e ± = π ± γ ± ρ ± + A ± γ ± − ρ ± ) γ ± − , (C.12)where A ± is a constant. Substituting (C.12) into the EOS (C.9) yields p ± + π ± γ ± = A ± ( ρ ± ) γ ± , (C.13)that is P ± I ( ρ ± ) = − π ± γ ± + A ± ( ρ ± ) γ ± . (C.14)Thus R ± I ( p ) = p + π ± /γ ± A ± ! /γ ± , (C.15)and the equation which gives P I ( α, ρ ) is ρ = 1 + α p + π + /γ + A + ! /γ + + 1 − α p + π − /γ − A − ! /γ − . (C.16)Even in the special case of two perfect gases where π ± = 0, this equation isin general transcendental. This is to be contrasted with the general case (thecase with a variable entropy), where P ( α, ρ, e ) can be calculated explicitly byalgebraic equations.From now on, we denote the set of equations (C.1)–(C.4) by (E I ). In order tostudy small perturbations around basic smooth and stationary solutions, it is35ore convenient to use the set of isentropic equations (E I ) rewritten in thephysical variables α , ~u , p .These equations are given in the following proposition. Proposition 3
The equivalent system to (E I ) in variables α , ~u , p is ~u t + ~u · ∇ ~u + 1 ρ ∇ p = ~g , (C.17) p t + ~u · ∇ p + ρc Is ∇ · ~u = 0 , (C.18) α t + ~u · ∇ α + (1 − α ) δ I ∇ · ~u = 0 , (C.19) where c Is and δ I are given by c Is = Γ + Γ − Γ ρρρ + ρ − pρ , Γ = 1 + α − + 1 − α + , (C.20)Γ ± ≡ ρ ± p d P ± I ( ρ ± ) d ρ ± = γ ± + π ± p , ρ = α − ρ + + α + ρ − , (C.21) δ I = Γ + − Γ − . (C.22) Remark 7
In the one-fluid case (take for example α + = 1 , α − = 0 , α = 1 ),one finds c s = c Is = c + s , (C.23) while, in the two-fluid case, c s = c Is . The analysis for the dispersion relation is quite similar to the general case.The steady state is denoted by α ± , ρ ± , p and ~u . We look for a special classof solutions which are motionless, uniform in the horizontal coordinates andcontinuously stratified in the vertical direction: α ± = α ± ( z ) , ρ ± = ρ ± ( z ) , p = p ( z ) , ~u = ~ . Again we take ρ to be constant. One must solve(1 + α ( z )) R + I ( p ( z )) + (1 − α ( z )) R − I ( p ( z )) = 2 ρ , (C.24)in order to find α ( z ) and ρ ± . It is easy to see that α ( z ) = 2 ρ − R + I ( p ( z )) − R − I ( p ( z )) R + I ( p ( z )) − R − I ( p ( z )) , (C.25)with p ( z ) given by (48).The analysis is the same as before, except that c s and δ are replaced by thevalues for the isentropic case. 36 ontents c − s | ~k | < ω < c + s | ~k | ω < min( c − s | ~k | , c + s | ~k | ) 32B.3 Case where ω > max( c + s | ~k | , c − s | ~k | ) 33C Isentropic flows 34References 38 eferences [1] R.A. Bagnold. Interim report on wave pressure research. Proc. Inst. Civil Eng. ,12:201–26, 1939.[2] H. Bredmose. Flair: A finite volume solver for aerated flows. Technical report,2005.[3] G. N. Bullock, C. Obhrai, D. H. Peregrine, and H. Bredmose. Violent breakingwave impacts. Part 1: Results from large-scale regular wave tests on verticaland sloping walls.
Coastal Engineering , 54:602–617, 2007.[4] G.N. Bullock, A.R. Crawford, P.J. Hewson, M.J.A. Walkden, and P.A.D. Bird.The influence of air and scale on wave impact pressures.
Coastal Engineering ,42:291–312, 2001.[5] D. Dutykh.
Mathematical modelling of tsunami waves . PhD thesis, Centre deMath´ematiques et Leurs Applications, 2007.[6] J.-M. Ghidaglia, A. Kumbaro, and G. Le Coq. On the numerical solution to twofluid models via cell centered finite volume method.
Eur. J. Mech. B/Fluids ,20:841–867, 2001.[7] J.-M. Ghidaglia and F. Pascal. The normal flux method at the boundaryfor multidimensional finite volume approximations in CFD.
Eur. J. Mech.B/Fluids , 24:1–17, 2005.[8] S.K. Godunov, A. Zabrodine, M. Ivanov, A. Kraiko, and G. Prokopov.
R´esolution num´erique des probl`emes multidimensionnels de la dynamique desgaz . Editions Mir, Moscow, 1979.[9] M. Ishii.
Thermo-Fluid Dynamic Theory of Two-Phase Flow . Eyrolles, Paris,1975.[10] D.H. Peregrine and L. Thais. The effect of entrained air in violent water impacts.
J. Fluid Mech. , 325:377–97, 1996.[11] G. A. Sod. A survey of several finite difference methods for systems of nonlinearhyperbolic conservation laws.
J. Comput. Phys. , 43:1–31, 1978.[12] B. van Leer. Upwind and high-resolution methods for compressible flow: Fromdonor cell to residual-distribution schemes.
Communications in ComputationalPhysics , 1:192–206, 2006.[13] D.J. Wood, D.H. Peregrine, and T. Bruce. Wave impact on wall using pressure-impulse theory. I. Trapped air.
Journal of Waterway, Port, Coastal and OceanEngineering , 126(4):182–190, 2000., 126(4):182–190, 2000.