A hybrid WENO method with modified ghost fluid method for compressible two-medium flow problems
AA hybrid WENO method with modified ghost fluidmethod for compressible two-medium flow problems , Yibing Chen and Jianxian Qiu Abstract
In this paper, we develop a simplified hybrid weighted essentially non-oscillatory (WENO)method combined with the modified ghost fluid method (MGFM) [28] to simulate the com-pressible two-medium flow problems. The MGFM can turn the two-medium flow problemsinto two single-medium cases by defining the ghost fluids status in terms of the predicted theinterface status, which makes the material interface “invisible”. For the single medium flowcase, we adapt between the linear upwind scheme and the WENO scheme automatically byidentifying the regions of the extreme points for the reconstruction polynomial as same asthe hybrid WENO scheme [50]. Instead of calculating their exact locations, we only needto know the regions of the extreme points based on the zero point existence theorem, whichis simpler for implementation and saves computation time. Meanwhile, it still keeps therobustness and has high efficiency. Extensive numerical results for both one and two dimen-sional two-medium flow problems are performed to demonstrate the good performances ofthe proposed method.
Key Words: hybrid WENO scheme, two-medium flow problems, modified ghost fluidmethod, zero point existence theorem
AMS(MOS) subject classification: The research is partly supported by Science Challenge Project, No. TZ2016002 and NSAF grantU1630247. School of Mathematical Sciences, Xiamen University, Xiamen, Fujian 361005, P.R. China. E-mail:[email protected]. Institute of Applied Physics and Computational Mathematics, Beijing 100094, China. E-mail:chen [email protected]. School of Mathematical Sciences and Fujian Provincial Key Laboratory of Mathematical Modeling andHigh-Performance Scientific Computing, Xiamen University, Xiamen, Fujian 361005, P.R. China. E-mail:[email protected]. a r X i v : . [ phy s i c s . c o m p - ph ] S e p Introduction
In this paper, we propose a simplified hybrid weighted essentially non-oscillatory (WENO)method with modified ghost fluid method (MGFM) [28] for simulating compressible two-medium flow problems. For two-medium flow problems, the equation of state (EOS) wouldswitch between the different medium, which may cause numerical oscillations or inaccura-cies near the material interface. Hence, many researchers have used various additional worksand techniques to overcome this difficulty, and there are two major options to simulate thecompressible two-medium flow problems.One is the front capturing method, where the high resolution methods are applied tosuppress the non-physical oscillations near discontinuities by bringing the numerical diffusionor viscosity, which inherently exists in the method itself or is given artificially. The frontcapturing method can deal with large deformation problems and relatively easy to extendto high dimension. However, the numerical inaccuracies and oscillations are inevitable nearthe interface, therefore, various techniques were introduced by Larrouturou [20], Karni [19],Abgrall et al. [1, 2], Shyue et al. [40], Saurel et al. [36] and Chen et al. [5] to resolvethis difficulty. The other one is the front tracking method, which terms the discontinuitiesbetween the two-medium flows as internal moving interfaces. It works well at multi-materialinterfaces, but it would have difficulties about the entanglement of the Lagrangian meshesand the extension to high dimension, and there are some typical methods about the fronttracking method, such as volume of fluid (VOF) method [13], level set method [41] and otherfront tracking methods [42, 7].To combine the best properties of the front capturing and tracking methods, Fedikw etal. [6] proposed a new numerical method for treating interfaces using a level set functionin Eulerian schemes named as the ghost fluid method (GFM), which makes the interface“invisible”. In the framework of the GFM [6], the pressure and velocity at the ghost fluidnodes near the interface are defined as the local real pressure and velocity, while the den-sity is obtained by isobaric fixing. It can easily turn the two-medium flow problems into2wo single-medium flow cases, and for the single-medium flow problems, many classical andmature schemes can be applied. Hence, it provides an alternative and flexible way to sim-ulate the two-medium flow problems, and the extension to high dimension becomes fairlystraightforward. However, it may cause numerical inaccuracies in the case of a strong shockimpacting on the interface, and the reason may be that the statuses near the interface areaffected by the wave interaction and the material properties on both sides. Therefore, Liuet al. [28] developed a modified ghost fluid method (MGFM), in which a multi-materialRiemann problem is defined and solved approximately or exactly to predict the interfacialstatus, then, the predicted interfacial status is applied to define the fluid values at the ghostpoints. The MGFM combines the advantages of the GFM [6] and the implicit characteristicmethods [26, 27], and it takes the interaction of shock with the interface into consideration.Later, the interface interaction GFM (IGFM) [14], the real GFM (RGFM) [43] and the prac-tical GFM (PGFM) [44] were developed following the idea of the Riemann problem-basedtechnique in the MGFM [28]. The MGFM is robust and less problem related, and it hasbeen applied in various situations as in [25, 30, 33, 46, 23], and the accuracy analysis anderrors estimation can be seen in [24, 45]. The GFM [6] and its relevant ghost fluid methods[28, 14, 43, 44] are non-conservative near the interface, and the conservative scheme can beseen in [32].Here, we would use the MGFM to define the ghost fluid status for the two-medium flowproblems considering its great performances, and for the single-medium flow problems, manysuccessful numerical schemes can be applied for it. Among them, the finite difference or finitevolume weighted essentially non-oscillatory (WENO) schemes have been widely applied forthe single-medium flow problems which usually contain shock, contact discontinuities andsophisticated smooth structures simultaneously. And in 1994, the first WENO scheme wasconstructed by Liu, Osher and Chan [29] on the basis of ENO schemes [11, 9, 10], where allcandidate stencils were used with a nonlinear convex methodology to obtain higher orderaccuracy in the smooth regions, then, Jiang and Shu [17] proposed the third and fifth-order3nite difference WENO schemes in multi-space dimension, in which a general definition forsmoothness indicators and nonlinear weights was presented. After this, WENO schemeshave been further developed in [15, 21, 37, 47, 3, 49, 51]. However, the cost of computingthe nonlinear weights and local characteristic decompositions is still very high. Hence, Hilland Pullin [12] combined the tuned center-difference schemes with WENO schemes to expectthat the nonlinear weights would be achieved automatically in the smooth regions away fromshocks, but a switching principle was still significant. Later, Li and Qiu [22] studied the hy-brid WENO scheme using different switching principles [34], which shows different principleswould have different influences for the hybrid WENO scheme [22], and the majorities of thetroubled-cell indicators need to adjust parameters for different problems to balance betternon-oscillations near discontinuities and less computation cost, simultaneously. Hence, Qiuet al. [50, 48] used a new simple switching principle, which employed different reconstructionmethod automatically based on the locations of all extreme points of the big reconstructionpolynomial for numerical flux. Then, we develop this methodology in this paper, in whichwe only need to know the regions of the extreme points, rather than calculating their exactlocations as in [50, 48].In this paper, to keep the robustness of the MGFM [28] and high efficiency of the hy-brid WENO schemes [50, 48], we first use the methodology introduced in [28] to predictthe interfacial status based on a multi-material Riemann problem, then, the predicted in-terfacial status is applied to define the ghost fluid values, by which it turns a two-mediumflow problems into two single-medium cases. For the single-medium problems, we wouldsolve it by the hybrid WENO method, where we reconstruct the numerical flux by upwindlinear approximation directly if all extreme points of the big reconstruction polynomial fornumerical flux are located outside of the big stencil, otherwise we use the classical WENOprocedure [17]. But we only need to know the regions of the extreme points in terms of thezero point existence theorem, instead of calculating their exact locations as in [50, 48], andit is more easy for implementation and saves computation time. Meanwhile, it still keeps the4obustness of the WENO scheme [17] and the MGFM [28] to simulate the two-medium flowproblems. In addition, it has higher efficiency with less computation costs than the WENOscheme [17] for employing linear approximation straightforwardly in the smooth regions.The organization of the paper is as follows: in Section 2, the detailed implementationprocedures of the finite difference hybrid WENO scheme combined with the MGFM arepresented for two-medium flow problems. In Section 3, Extensive numerical results for gas-gas and gas-water interaction problems in one and two dimensions are presented to illustrategood performances of the proposed scheme. Concluding remarks are given in Section 4.
We first introduce the governing equations for the compressible two-medium flow prob-lems, then, we give a brief review about the finite difference hybrid WENO method [50] forsingle-medium flow problems, but we have an improvement about the identification tech-nique for the regions of the extreme points of the big reconstruction polynomial. Next, weintroduce the level set technique to track the moving interface, then, we briefly introducethe modified ghost fluid method (MGFM) [28] to define the status of ghost fluids. Finally,we give the summary of the implementation procedures.
We consider the hyperbolic conservations laws given as follows (cid:26) U t + ∇ · F ( U ) = 0 ,U ( x,
0) = U ( x ) , (2.1)where U is ( ρ, ρµ, E ) T and F ( U ) is ( ρµ, ρµ + p, µ ( E + p )) T for one dimensional Eulerequations, while for two dimensional Euler equations, U is ( ρ, ρµ, ρν, E ) T , and F ( U ) is( F ( U ) , F ( U )) with F ( U ) = ( ρµ, ρµ + p, ρµν, µ ( E + p )) T , F ( U ) = ( ρν, ρµν, ρν + p, ν ( E + p )) T . In order to close the systems, the equations of state (EOS) is still required. The γ -lawfor gas is ρe = p/ ( γ − , ρe = ( p + N ¯ B ) / ( N − , in which ¯ B = B − A , N = 7 . A = 1 . × Pa, A = 3 . × Pa, and ρ = 1000 . kg/m . We first consider one dimensional scalar hyperbolic conservation laws (cid:26) u t + f x ( u ) = 0 ,u ( x,
0) = u ( x ) . (2.1)We divide the computing domain by uniform grid points { x i } , and h is denoted as x i +1 − x i .The cell I i is defined as [ x i − / , x i +1 / ], where x i +1 / is set as x i +1 / = x i + h/
2, then, thesemi-discrete finite difference scheme of (2.1) is written as du i ( t ) dt = − h (cid:16) ˆ f i +1 / − ˆ f i − / (cid:17) , (2.2)in which u i ( t ) is represented as u ( x i , t ), and the numerical flux ˆ f i +1 / is a fifth order approx-imation of v i +1 / = v ( x i +1 / ), in which v ( x ) is defined implicitly as in [17]: f ( u ( x )) = 1 h (cid:90) x + h/ x − h/ v ( x ) dx, then, the right hand item of (2.2) is the fifth order approximation for − f x ( u ) at x i . Forthe stability of the finite difference scheme, we split the flux f ( u ) into two parts: f ( u ) = f + ( u ) + f − ( u ), in which df + ( u ) du ≥ df − ( u ) du ≤
0, and the Lax-Friedrichs flux splittingmethod is applied here as f ± ( u ) = 12 ( f ( u ) ± αu ) , where α = max u | f (cid:48) ( u ) | .Next, we introduce the detailed procedures for the reconstruction of the numerical fluxˆ f + i +1 / , which is the fifth order approximation of f + ( u ( x i +1 / )), and the reconstruction formu-las for ˆ f − i +1 / are mirror symmetric with respect to x i +1 / of that for ˆ f + i +1 / . ˆ f i +1 / is finally6aken as ˆ f + i +1 / + ˆ f − i +1 / . Now, we first give a big stencil: S = { I i − , ..., I i +2 } , then we caneasily obtain the fourth degree polynomial p ( x ) in terms of the following requirements as1 h (cid:90) I j p ( x ) dx = f + ( u j ) , j = i − , ..., i + 2 . For simplicity, ( x − x i ) h is set as ξ , then we have p ( x ) = 11920 [( − f + i − + 9 f + i − + 2134 f + i − f + i +1 + 9 f + i +2 ) − f + i − − f + i − − f + i +1 + 5 f + i +2 ) ξ + 120(12 f + i − − f + i − − f + i + 12 f + i +1 − f + i +2 ) ξ + 160(2 f + i − − f + i − − f + i +1 + f + i +2 ) ξ − f + i − − f + i − − f + i + 4 f + i +1 − f + i +2 ) ξ )] . To increase the efficiency, we also use the thought of the hybrid WENO schemes [50, 48],in which the linear upwind approximation or WENO reconstruction is applied automaticallybased on the locations of the extreme points of the big polynomial p ( x ). More explicitly, ifall extreme points are located outside of the big spatial stencil S , ˆ f + i +1 / is taken as p ( x i + )directly, otherwise the classical WENO procedures [17] would be used to reconstruct it.Unlike the hybrid WENO schemes [50, 48] solving the real zero points of p (cid:48) ( x ) exactly,we identify the regions of the extreme points of p ( x ) in terms of the zero point existencetheorem as that if the endpoint values of p (cid:48) ( x ) and the extreme values of p (cid:48) ( x ) have samesigns on the big stencil S , it means there is no single zero points of p (cid:48) ( x ) on S , that is tosay, there is no extreme points of p ( x ) located in S . Also, we present their performancesin the numerical examples, which shows the new identification skill can catch the regions forthe extreme points of the big polynomial p ( x ) as the old one in [50, 48], but it has higherefficiency. In addition, the new one has simpler implementation procedure as it only needsto solve the zero points of the quadratic polynomial p (cid:48)(cid:48) ( x ), while the old one has to calculatethe roots of the cubic polynomial p (cid:48) ( x ).Then, we review the classical WENO procedure [17] for the reconstruction of ˆ f + i +1 / .Firstly, the big stencil S is divided into three smaller stencils: S = { I i − , I i − , I i } , S = { I i − , I i , I i +1 } and S = { I i , I i +1 , I i +2 } , then, the polynomials p l ( x ) are obtained by the7ollowing requirements as1 h (cid:90) I j p l ( x ) dx = f + ( u j ) , j = i − l, ..., i − l, l = 1 , , . The explicit values of p l ( x ) at the point x i +1 / can be seen in [17], and the linear weightsare computed by p ( x i +1 / ) = (cid:80) l =1 γ l p l ( x i +1 / ), in which γ = , γ = and γ = . Tomeasure how smooth these small polynomials p l ( x ) are in the target cell I i , we use the samedefinition of smoothness indicators β l seen in [17, 38] as β l = (cid:88) α =1 (cid:90) I i h α − ( d α p l ( x ) dx α ) dx, then, the nonlinear weights are ω l = ω l (cid:80) rk =0 ω k , ω l = γ l ( β l + ε ) , l = 1 , , , where ε = 10 − , and the explicit values of ω l also can be seen in [17]. Finally, the WENOreconstruction of ˆ f + i +1 / is ˆ f + i +1 / = (cid:88) l =1 ω l p l ( x i +1 / ) . After the spatial discretization, the semi-discrete scheme (2.2) is discretized by the thirdorder TVD Runge-Kutta method [39] in time as u (1) = u n + ∆ tL ( u n ) ,u (2) = u n + u (1) + ∆ tL ( u (1) ) ,u ( n +1) = u n + u (2) + ∆ tL ( u (2) ) . (2.3) Remark 1:
For the systems, such as the one dimensional compressible Euler equations,WENO procedure is performed in the local characteristic directions to overcome the oscilla-tions nearby discontinuities as in [17], while the linear approximation is directly computedin each component. For two dimensional problems, the spatial reconstruction is performedby dimension by dimension.
We choose the next level set technique to track the moving interface. For one dimensionalproblems, the level set equation is φ t + uφ x = 0 , (2.4)8hile for two dimensional case it is written as φ t + uφ x + vφ x = 0 , (2.5)where φ is a signed distance function. u and v are the velocity of the flow in the x and y directions, respectively. We would solve the equations (2.4) and (2.5) by the fifth orderfinite difference hybrid WENO method introduced in Appendix A. However, if the velocityfield has a large gradient in the vicinity of the interface, the level set method may causeseriously distorted contours. Therefore, the re-initialization technique is needed to remedythis influence. For one dimensional problems, we can obtain the position of the interfaceexactly by Newton’s iteration method, then we re-distribute the signed distance function φ . For two dimensional case, the interface is a curve, so we need to use other ways forre-initialization, in which we solve the re-initialization equation as: φ t + S ( φ )( (cid:113) φ x + φ y −
1) = 0 , (2.6)where S is the sign function of φ , and the equation (2.6) is also calculated by the hybridWENO method shown in Appendix A. We use the modified ghost fluid method (MGFM) [28] to define the information of theghost cells as it considers the interaction of shock with the interface correctly. The mainprocedures of the MGFM are that we first predict the interface status by solving a two-medium Riemann problem exactly or approximately, then, the predicted interface status isused to define the ghost fluid status for each fluid, by which it turns a two-medium flowproblem into two single-medium flow problems.For one dimensional case, we only introduce how to define the ghost fluid status forMedium 1 in detail, and the definition of the ghost fluid status for Medium 2 is similar.Let’s assume that the interface is located between i and i + 1 seen in Figure 2.1, then, we usethe status of U i − and U i +2 to define the two-medium Riemann problem suggested in [28],9 nterface ii − i − i + 1 i + 2 i + 3 u I p I ρ LI ρ RI Medium 1 Medium 2
RealGhost u I velocity at interface p I pressure at interface ρ LI density at left-side interface ρ RI density at right-side interface Figure 2.1: Isentropic fixing for 1D two-medium flow problems.and obtain the interfacial status: u i (velocity), p i (pressure), ρ Li (density at left-side) and ρ Ri (density at right-side). Later, we take the predicted u i , p i , ρ Li as the velocity, pressureand density at the ghost point i + 1, but at these points i , i + 2 and i + 3, the pressure andvelocity are those on the real local fluid, and the density at these points is replaced by theisentropic fixing [6, 28].For two dimensional case, it would have one difficulty about the definition of the two-medium Riemann problem for there is two velocity components. However, we can know thenormal direction −→ n near the interface employing the level set function ( −→ n = ∇ φ/ |∇ φ | ),then, we obtain the normal velocity u N and tangential velocity u T , in which u N is defined as( µ, ν ) · −→ n , then, we apply the normal velocity u N , the pressure p and the density ρ to definethe two-medium Riemann problem like one dimensional case. In terms of the MGFM [28],we need to define a computation domain for each medium that includes boundary pointsand grid points in the interfacial regions by | φ | < (cid:15) , where (cid:15) is set to be 3max(∆ x, ∆ y ) forthe fifth order hybrid WENO scheme. Later, we would only introduce the definition of thestatus for Medium 1 at the points A and B (seen in Figure 2.2) in detail. To define thestatus at the point A in Medium 1, we need to find other point next to the interface ( | φ | < (cid:15) )located in the Medium 2, and let’s assume that B is the target point as the angle made bythe normal of B and A is the minimum, then we define Riemann problem in the normaldirection as U | t = t n = (cid:26) U A ,U B , in which U A = ( ρ A , u AN , p A ) and U B = ( ρ B , u BN , p B ), then we can solve it approximately or10 B Medium 1 Medium 2Interface
RealGhost
Figure 2.2: Isentropic fixing for 2D two-medium flow problems.exactly to predict the status u I (velocity), p I (pressure), ρ LI (density at left-side) and ρ RI (density at right-side). Notice that node A is located in Medium 1, then, we only need todefine the density by isentropic fixing [6, 28], but node B is located in Medium 2, therefore,we need to define the status at node B by u I , p I and ρ LI , and its tangential velocity is stillthe original one. In addition, the definition of the ghost fluid status for Medium 2 is similar. Now, we give a brief summary of the procedures for simulating two-medium flow prob-lems. Let’s assume that the flow status at t n has been obtained, then we can advance therespective quantities to t n +1 following as: Step 1.
Calculate the time step ∆ t , satisfying the stability condition over the wholerange. Step 2.
Solve the level set function φ , and obtain the locations of the interface in thenext intermediate time step introduced in Section 2.3. Step 3.
Define the Riemann problem near the interface and predict the interface status,then use it to define the fluid values at the ghost points for Mediums 1 and 2, respectively,given in Section 2.4.
Step 4.
Solve the governing equations for Mediums 1 and 2, respectively, advancing the11olution to the next intermediate time level, shown in Section 2.2.
Step 5.
Repeat Steps 2, 3 and 4 at each intermediate time step of the third order TVDRunge-Kutta method, and advance the solution from U n to U n +1 , then, re-initialize the levelset function φ . In this section, we perform the numerical results of the new simplified hybrid WENOscheme and classical WENO scheme [17] with the modified ghost fluid method for two-medium flow problems which are outlined in the previous section, and the CFL number isset as 0.6. In addition, we also make a comparison between the new identification skill forthe regions of the extreme points introduced in the previous section and the old one used inthe hybrid WENO schemes [50, 48]. The units for the density, velocity, pressure, length andtime are kg / m , m/s, Pa, m, and s, respectively.Here, we use “New/simplified hybrid WENO method” to denote the new simplified finitedifference hybrid WENO scheme with the modified ghost fluid method developed in thispaper, and use “Classical WENO method” to represent the classical finite difference WENOscheme [17] with the modified ghost fluid method. In addition, we use “Old hybrid WENOmethod” to denote the finite difference hybrid WENO scheme [50] with the modified ghostfluid method. Example 3.1.
This problem was taken from [6], and the initial conditions are( ρ, µ, p, γ ) = (cid:26) (1 , , × , . , x ∈ [0 , . , (0 . , , × , . , x ∈ [0 . , . In flow and out flow boundary conditions are applied here, and the final computed time t is 0.0007. We present the computed density ρ , velocity µ and pressure p by New/Simplifiedhybrid WENO and Classical WENO methods against the exact solution in Figure 3.1. Wecan find that the two methods both capture the location of the material interface correctly.These two schemes also have similar numerical results, and the overall results are compa-rable to analysis, but New/simplified hybrid WENO method can achieve higher efficiency12 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X D e n s i t y ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X V e l o c i t y ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X P r ess u r e Figure 3.1: Example 3.1. t=0.0007. From left to right: density; velocity; pressure. Solidline: the exact solution; plus signs: the results of Classical WENO method; squares: theresults of New/simplified hybrid WENO method. Grid points: 200.than Classical WENO method for saving 22.93% computation time. In addition, we findNew/simplified hybrid WENO method can save 9.25% CPU time than Old hybrid WENOmethod by calculation, meanwhile, there are only 15.13% and 15.14% points where theWENO procedures are computed in New/simplified hybrid WENO and Old hybrid WENOmethods, respectively, and the time history of the locations of WENO reconstruction for twomethods are given in the top of Figure 3.4. These results show that the new identificationskill in New/simplified hybrid WENO method can identify the regions of the extreme pointscorrectly as the old one in Old hybrid WENO method, but the new one has higher efficiency.The new identification technique is also simpler as it only needs to solve the zero points ofa quadratic polynomial, while the old one has to calculate the roots of a cubic polynomial.
Example 3.2.
This problem is also taken from [6], which contains a right going shock re-fracting at an air-helium interface with a reflected rarefaction wave, and the initial conditionsare given as( ρ, µ, p, γ ) = (4 . , . √ , . × , . , x ∈ [0 , . , (1 , , × , . , x ∈ [0 . , . , (0 . , , × , / , x ∈ [0 . , , with inflow/outflow boundary conditions. The initial strength of the shock is p l /p R = 15 at x = 0 .
05, and the interface of air and helium is located at x = 0 .
5. We ran the code to a finaltime of 0.0005, and the computed density ρ , velocity µ and pressure p by New/simplifiedhybrid WENO and Classical WENO methods against the exact solution are shown in Figure13 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X D e n s i t y ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X V e l o c i t y ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X P r ess u r e Figure 3.2: Example 3.2. t=0.0005. From left to right: density; velocity; pressure. Solidline: the exact solution; plus signs: the results of Classical WENO method; squares: theresults of New/simplified hybrid WENO method. Grid points: 200.3.2. We can see the contact discontinuity is located in the correct cell, and two methods havesimilar results, but New/simplified hybrid WENO method can achieve higher efficiency thanClassical WENO method for saving 16.41% CPU time. In addition, we find New/simplifiedhybrid WENO method can save 9.54% computation time than Old hybrid WENO method bycalculation, meanwhile, there are 25.07% and 24.50% points where the WENO procedures arecomputed in New/simplified hybrid WENO and Old hybrid WENO methods, respectively,and the time history of the locations of WENO reconstruction are shown in the middle ofFigure 3.4, which illustrate that the new identification skill in New/simplified hybrid WENOmethod catchs the regions of the extreme points as the old one in Old hybrid WENO method.However, the new identification technique has higher efficiency, and it is also simpler for itonly needs to calculate the roots of a quadratic polynomial.
Example 3.3.
We solve the governing equations (2.1) for one dimensional Euler equationswith the following Riemann initial conditions( ρ, µ, p, γ ) = (1 . , . √ , . × , . , x ∈ [0 , . , (1 , , × , . , x ∈ [0 . , . , (3 . , , × , . , x ∈ [0 . , , The final computed time t is up to 0.0017. This example is also taken from [6], and thecomputed density ρ , velocity µ and pressure p by New/simplified hybrid WENO and ClassicalWENO methods against the exact solution are given in Figure 3.3. The numerical resultsillustrate two schemes capture the contact discontinuity correctly, with non-oscillations and14 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X D e n s i t y ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X V e l o c i t y ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X P r ess u r e Figure 3.3: Example 3.3. T=0.0017. From left to right: density; velocity; pressure. Solidline: the exact solution; plus signs: the results of Classical WENO method; squares: theresults of New/simplified hybrid WENO method. Grid points: 200.similar comparable results, meanwhile, New/simplified hybrid WENO method saves almost25.49% computation time comparing with Classical WENO method. In addition, we findNew/simplified hybrid WENO method with the new identification skill can save 8.31% CPUtime than Old hybrid WENO method by calculation, meanwhile, there are only 8.48% and8.52% points where the WENO procedures are computed in New/simplified hybrid WENOand Old hybrid WENO methods, respectively. The time history of the locations of WENOreconstruction by the two hybrid methods are seen in the bottom of Figure 3.4. These resultsillustrate that the new identification technique in New/simplified hybrid WENO methodcatches the regions of the extreme points as the old one in Old hybrid WENO method, butNew/simplified hybrid WENO method with the new one has higher efficiency, and the newidentification technique is also simpler.
Example 3.4.
This problem is taken from [28], having a strong shock on a gas-gas interface,and the strength of the right shock wave is up to p L /p R = 100. The initial conditions aregiven as follows( ρ, µ, p, γ ) = (0 . , . √ , . × , / , x ∈ [0 , . , (0 . , , × , / , x ∈ [0 . , . , (1 , , × , . , x ∈ [0 . , . In Figure 3.5, we present the computed density ρ , velocity µ and pressure p by New/simplifiedhybrid WENO and Classical WENO methods against the exact solution at the final time0 . +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X D e n s i t y ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X V e l o c i t y ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X P r ess u r e Figure 3.5: Example 3.4. t=0.0001. From left to right: density; velocity; pressure. Solidline: the exact solution; plus signs: the results of Classical WENO method; squares: theresults of New/simplified hybrid WENO method. Grid points: 200.lem, and they capture the correct location of the interface between two gases. Comparingwith Classical WENO method, New/simplified hybrid WENO method saves almost 15.98%computation time. In addition, we find New/simplified hybrid WENO method with the newidentification skill can save 14.52% CPU time than Old hybrid WENO method by calcu-lation, meanwhile, there are 26.76% and 25.52% points where the WENO procedures areperformed in New/simplified hybrid WENO and Old hybrid WENO methods, respectively,and the time history of the locations of WENO reconstruction by two hybrid methods aregiven in the top of Figure 3.8, which illustrate the new identification skill in New/simplifiedhybrid WENO method identifies the regions of the extreme points as the old one in Oldhybrid WENO method, but New/simplified hybrid WENO method with the new one hashigher efficiency. The new identification technique in New/simplified hybrid WENO methodis also simpler as it only needs to solve the roots of a quadratic polynomial, while the oldone in Old hybrid WENO method has to calculate the zero points of a cubic polynomial.
Example 3.5.
We consider the gas-water shock tube problem taken from [33], and theinitial condition are given as( ρ, µ, p, γ ) T = (cid:26) (1270 , , × , . T , x ∈ [0 , . , (1000 , , × , . T , x ∈ [0 . , . This underwater explosion problem has extremely high pressure in the gas medium, therefore,there is a very strong shock in the water. The final computation time is 0 . +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X D e n s i t y ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X V e l o c i t y ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X P r ess u r e Figure 3.6: Example 3.5. t=00016. From left to right: density; velocity; pressure. Solid line:the exact solution; plus signs: the results of Classical WENO method; squares: the resultsof New/simplified hybrid WENO method. Grid points: 200.present the computed density ρ , velocity µ and pressure p by New/simplified hybrid WENOand Classical WENO methods against the exact solution in Figure 3.6, which illustratestwo schemes capture the location of the material interface correctly, and they have goodperformance in the smooth and discontinuous regions. In addition, New/simplified hybridWENO method saves about 16.53% CPU time. Moreover, we find New/simplified hybridWENO method with the new identification skill can save 13.93% computation time than theOld hybrid WENO method with the old identification technique by calculation, meanwhile,there are both 20.68% points where the WENO procedures are performed in New/simplifiedhybrid WENO and Old hybrid WENO methods, respectively, and the time history of thelocations of WENO reconstruction in two hybrid methods are given in the middle of Figure3.8. These results show that the new identification skill in New/simplified hybrid WENOcan catch the regions of the extreme points as the old one in Old hybrid WENO method.However, New/simplified hybrid WENO method with the new identification skill has higherefficiency, and the procedure of the new one is also simpler. Example 3.6.
This gas-water shock tube problem is taken from [6], which has higherenergy of the explosive gaseous medium than the problem given in Example 3.5, and theinitial conditions are( ρ, µ, p, γ ) T = (cid:26) (1630 , , . × , . T , x ∈ [0 , . , (1000 , , × , . T , x ∈ [0 . , . +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X D e n s i t y ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X V e l o c i t y ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X P r ess u r e Figure 3.7: Example 3.6. t=0.0001. From left to right: density; velocity; pressure. Solidline: the exact solution; plus signs: the results of Classical WENO method; squares: theresults of New/simplified hybrid WENO method. Grid points: 200.We ran the code to the final time 0.0001, then, we give the computed density ρ , velocity µ and pressure p by New/simplified hybrid WENO and Classical WENO methods againstthe exact solution in Figure 3.7, which shows two schemes work well for this tough gas-water problem with non-oscillation in the non-smooth regions, and they also capture theright interface between two mediums. In addition, New/simplified hybrid WENO methodhas higher efficiency for saving almost 15.12% CPU time. Moreover, we find New/simplifiedhybrid WENO method with the new identification skill can save 13.52% computation timethan the Old hybrid WENO method by calculation, meanwhile, there are both 32.88%points where the WENO procedures are computed in two hybrid WENO method, and thetime history of the locations of WENO reconstruction in New/simplified hybrid WENO andOld hybrid WENO methods are given in the bottom of Figure 3.8, which show that thenew identification skill New/simplified hybrid WENO method can catch the regions of theextreme points correctly as the old one, but New/simplified hybrid WENO method with thenew one has higher efficiency, and the new identification technique is also simpler. Example 3.7.
This problem is a Mach 1.22 air shock acting on a helium bubble, then,we solve the governing equations (2.1) for two dimensional Euler equations, and its physicalinitial schematic diagram is shown in the left of Figure 3.9. The reflective conditions areapplied in the upper and lower boundary, while the inflow/outflow conditions are given in the19igure 3.8: The points where the WENO procedures are performed for Examples 3.4 to 3.6.From left to right: the results of Old hybrid WENO method; the results of New/simplifiedhybrid WENO method. 20eft and right boundary, respectively. The non-dimensionalized initial conditions are givenas follows ( ρ, µ, ν, p, γ ) = (1 , , , , . , pre-shocked air , (1 . , . , , . , . , post-shocked air , (0 . , , , / , helium ,φ = (cid:112) x + y − , level set , in which φ ≤ φ > x < . .
5, 1 .
0, 2 . . × t = 4 . t = 0 .
5, we can see that the initial regular shock becomes irregular having bifurcation ofthe shock on the bubble surface for the sound speed in helium is faster than that in air,which was also illustrated in the experimental results given in [8]. At t = 1 .
0, the refractedshock inside the bubble has interacted with the rear of the bubble and enters into the air,but the incident shock just went through the top of the helium bubble, then, the wholebubble starts moving to the right. At t = 2 .
0, the incident shock has gone through the wholebubble, and the shape of the bubble begins to misshape. After this, a re-entrant jet beginsto form. At t = 4 .
0, the re-entrant jet actually has been formed, and the interface would beinstable, when the re-entrant jet becomes stronger and stronger, which would affect the rearside of the bubble and cause the bubble to collapse, and the quite fine meshes or adaptiverefinement technique might be needed seen in [18, 35]. Hence, our computation stops at thenon-dimensional time t = 4 .
0. 21 -3,-3) (4,3)x=-1.2Post-shocked Pre-shocked airair Heliumbubble (-4,-3) (3,3)x=-1.2Post-shockedwater Pre-shocked waterGas
Figure 3.9: Physical domain for Example 3.7 (left) and Example 3.8 (right)Finally, we also find that the results computed by New/simplified hybrid WENO andClassical WENO methods are similar, but New/simplified hybrid WENO method savesalmost 35.49% computation time as we use linear approximation directly for the governingequation, the level set function and its re-initialization in the smooth regions. In addition, wefind New/simplified hybrid WENO method with the new identification skill can save 9.41%computation time than Old hybrid WENO method by calculation, meanwhile, there areonly 2.84% and 2.75% points where the WENO procedures are computed in New/simplifiedhybrid WENO and Old hybrid methods at the final time step, respectively, and the locationsof WENO reconstruction computed by the two hybrid WENO methods at the final timestep are given in the top of Figure 3.12, which illustrate that the new identification skillhas similar ability as the old one in Old hybrid WENO method, but New/simplified hybridWENO method with the new one has higher efficiency, and the new one is simpler as it onlyneeds to solve the roots of a quadratic polynomial, while the old one has to calculate thezero points of a cubic polynomial.
Example 3.8.
The final example is a initial Mach 1.653 planar underwater shock inter-acting with a gas bubble in an open domain taken from [33], then we solve the governingequations (2.1) for two dimensional Euler equations with the next non-dimensionalized initial22 Y -3 -2 -1 0 1 2 3-3-2-1 X Y -3 -2 -1 0 1 2 3-3-2-10123 X Y -3 -2 -1 0 1 2 3-3-2-10123 X Y -3 -2 -1 0 1 2 3-3-2-10123 X Y -3 -2 -1 0 1 2 3-3-2-10123 X Y -3 -2 -1 0 1 2 3-3-2-10 X Y -3 -2 -1 0 1 2 3-3-2-10123 X Y -3 -2 -1 0 1 2 3-3-2-10123 Figure 3.10: Example 3.7. The results computed by Classical WENO method (left) andNew/simplified hybrid WENO method (right). 30 equally spaced density contours from 0.1to 1.6. From top to bottom are T = 0 . T = 1 . T = 2 . T = 4 .
0, respectively. Gridpoints: 280 × ρ, µ, ν, p, γ ) = (1000 , , , , . , pre-shocked water , (1176 . , . , , , . , post-shocked water , (1 , , , . , gas ,φ = (cid:112) x + y − , level set , where φ ≤ φ > x < . t = 0 . t = 0 . t = 0 .
357 and = 0 . .
481 before the form of thestrong re-entrant jet, and the bubble doesn’t collapse at this time.From Figure 3.11, we can see that the density computed by New/simplified hybrid WENOand Classical WENO methods are similar, however, New/simplified hybrid WENO methodhas higher efficiency than Classical WENO scheme for saving almost 27.13% computationtime. In addition, we find New/simplified hybrid WENO method with the new identificationskill can save 12.00% CPU time than the old one in Old hybrid WENO method by calcula-tion, meanwhile, there are only 19.05% and 19.21% points where the WENO procedures arecomputed in the two hybrid WENO methods at the final time step, respectively, and thelocations of WENO reconstruction at the final time step are shown in the bottom of Fig-ure 3.12, which illustrate that the new identification skill in New/simplified hybrid WENOmethod can identify the regions of the extreme points as the old one in Old hybrid WENO,but New/simplified hybrid WENO method with the new one has higher efficiency, and thenew one has simpler implementation procedure.24 Y -4 -3 -2 -1 0 1 2-3-2-10123 X Y -4 -3 -2 -1 X Y -4 -3 -2 -1 0 1 2-3-2-10 X Y -4 -3 -2 -1 0 1 2-3-2-10 X Y -4 -3 -2 X Y -4 -3 -2 - X Y -4 -3 -2 -1 0 1 2-3-2-10123 X Y -4 -3 -2 -1 0 1 2-3-2-10123 Figure 3.11: Example 3.8. The results computed by Classical WENO method (left) andNew/simplified hybrid WENO method (right). 30 equally spaced density contours from 1.0to 1200. From top to bottom are t = 0 . t = 0 . t = 0 .
357 and t = 0 . × In this paper, we combine the new simplified hybrid WENO method with the modifiedghost fluid method [28] to simulate the compressible two-medium flow problems, whichadapts between the linear upwind approximation and WENO reconstruction automaticallyin terms of the regions of the extreme points for the big reconstruction polynomial, andwe have an improvement about the identification technique for the regions of the extremepoints of the big reconstruction polynomial. This new switch principle is not only simplefor we doesn’t need to adjust the parameters, but also it is effective as we just need toknow the range of the extreme point for the big reconstruction polynomial, rather thansolving the exact location of the extreme point as the old one in hybrid WENO schemes[50, 48]. Comparing with the classical WENO scheme [17], the new simplified hybrid WENOscheme is more efficient with less numerical errors in the smooth region and less computationcosts, meanwhile, the new simplified hybrid WENO scheme with MGFM is robust and non-26scillatory to simulate these two-medium flow problems. In general, these numerical resultsall show the good performances of the new simplified hybrid WENO scheme with the modifiedghost fluid method.
The next finite difference hybrid WENO method for Hamilton-Jacobi equations is mainlydeveloped from the fifth order WENO scheme introduced by Jiang and Peng [16] to solvethe Hamilton-Jacobi equations (2.4), (2.5) and (2.6) in Section 2.2.We first consider one dimensional Hamilton-Jacobi equation (cid:26) φ t + H ( x, t, φ, φ x ) = 0 ,φ ( x,
0) = φ ( x ) . (5.1)The computing domain is divided by uniform grid points { x i } , and the semi-discrete formof (5.1) is dφ i ( t ) dt = − ˆ H ( x i , t, φ i , φ + x i , φ − x i ) , (5.2)where φ i ( t ) is represented as φ ( x i , t ), and φ ± x i are linear or WENO approximations for ∂φ ( x i ) ∂x .ˆ H is a numerical flux to approximate H , and we use the Lax-Friedrichs (LF) flux here as:ˆ H ( x, t, φ, u + , u − ) = H (cid:18) x, t, φ, u + + u − (cid:19) − α ( u + , u − ) u + − u − , where α is max u | H ( u ) | , where H is represented as the partial derivative of H with respectto φ x .Next, we only introduce the reconstruction procedures for φ − x i , and the reconstruction for φ + x i is mirror symmetric with respect to x i of that for φ − x i . Firstly, we can easily obtain thefourth degree polynomial p ( x ) to approximate φ x in terms of the requirements:1∆ x (cid:90) x i + l x i − l p ( x ) dx = 1∆ x (cid:90) x i + k x i − k φ x dx = φ x i + k − φ x i − k ∆ x , k = − , ..., . To increase the efficiency, if all extreme points of p ( x ) are located outside of the bigspatial stencil, φ − x i is taken as p ( x i ) directly, otherwise we would use the next classical27ENO procedures [16, 17], and the method to identify the regions of the extreme points for p ( x ) has been detailedly introduced in Section 2.2.Now, we would give a brief review of the WENO reconstruction for φ − x i . Similarly, weobtain three quadratic polynomials p l ( x ) to approximate φ x , satisfying1∆ x (cid:90) x i − k + l x i − k + l p l ( x ) dx = 1∆ x (cid:90) x i − k + l x i − k + l φ x dx = φ x i − k + l − φ x i − k + l ∆ x , k = − , , , l = 1 , , . For saving space, the explicit values of p l ( x i ), the linear weights γ l , the smoothness indicators β l , and the nonlinear weights ω l are not present here, and these expressions can be seen in[16, 17]. Finally, the WENO reconstruction of φ − x i is approximated by φ − x i = (cid:88) l =1 ω l p l ( x i ) . For the two dimensional Hamilton-Jacobi equation (cid:26) φ t + H ( x, y, t, φ, φ x , φ y ) = 0 ,φ ( x, y,
0) = φ ( x, y ) . (5.3)The computing domain is divided by uniform grid points { ( x i , y j ) } , and the semi-discreteform of (5.3) is dφ i,j ( t ) dt = − ˆ H ( x i , y j , t, φ i,j , φ + x,i,j , φ − x,i,j , φ + y,i,j , φ − y,i,j ) , (5.4)where φ i,j ( t ) is represented as φ ( x i , y j , t ). φ ± x,i,j and φ ± y,i,j are linear or WENO approximationsfor ∂φ ( x i ,y j ) ∂x and ∂φ ( x i ,y j ) ∂y , respectively. ˆ H is a numerical flux to approximate H , and we usethe Lax-Friedrichs (LF) flux here as:ˆ H ( x i , y j , t, φ i,j , u + , u − , v + , v − ) = H ( x i , y j , t, φ i,j , u + + u − , v + + v − − α u + − u − − β v + − v − , where α is max u | H ( u, v ) | and β is max v | H ( u, v ) | . H and H are represented as thepartial derivative of H with respect to φ x and φ y , respectively. Finally, φ ± x,i,j and φ ± y,i,j arereconstructed by dimension by dimension as one dimension case.For the time discretization, the semi-discrete schemes (5.2) and (5.4) are discretized bythe third order Runge-Kutta method [39] in (2.3) of Section 2.2.28 eferences [1] R. Abgrall, How to prevent oscillations in multicomponent flow calculations: a quasiconservative approach, J. Comput. Phys., 125 (1996), 150-160.[2] R. Abgrall and S. Karni, Computations of compressible multifluids, J. Comput. Phys.,169 (2001), 594-623.[3] M. Castro, B. Costa and W.S. Don, High order weighted essentially non-oscillatoryWENO-Z schemes for hyperbolic conservation laws, J. Comput. Phys., 230 (2011),1766-1792.[4] T.-J. Chen and C.H. Cooke, On the Riemann problem for liquid or gas-liquid media,Int. J. Numer. Meth. Fluids, 18 (1994), 529-541.[5] Y. Chen and S. Jiang, A non-oscillatory kinetic scheme for multi-component flows withequation of state for a stiffened gas, J. Comput. Math., 29 (2011), 661-683.[6] R. P. Fedkiw, T. Aslam, B. Merriman and S. Osher, A non-oscillatory Eulerian approachto interfaces in multimaterial flows (the ghost fluid method), J. Comput. Phys., 152(1999), 457-492.[7] J. Glimm, J.W. Grove, X.L. Li, K.-M. Shyue, Y. Zheng and Q. Zhang, Three-dimensional front tracking, SIAM J. Sci. Comput., 19 (1998), 703-727.[8] J.-F. Haas and B. Sturtevant, Interaction of weak shock waves with cylindrical andspherical gas inhomogeneities, J. Fluid Mech., 181 (1987), 41-76.[9] A. Harten, Preliminary results on the extension of ENO schemes to two-dimensionalproblems, in Proceedings, International Conference on Nonlinear Hyperbolic Prob-lems, Saint-Etienne, 1986, Lecture Notes in Mathematics, edited by C. Carasso et al.et al.