Open Journal of Fluid Dynamics
Vol.4 No.1(2014), Article ID:44079,25 pages DOI:10.4236/ojfd.2014.41004

Sharpening Diffuse Interfaces with Compressible Flow Solvers

Nicolas Favrie1, Sergey Gavrilyuk1, Boniface Nkonga2, Richard Saurel1

1Aix-Marseille University, Marseille, France

2J.A. Dieudonné Mathematics Lab., University Nice-Sophia Antipolis, Nice, France

Email: nicolas.favrie@univ-amu.fr, sergey.gavrilyuk@univ-amu.fr, boniface.nkonga@unice.fr, richard.saurel@univ-amu.fr

Copyright © 2014 by authors and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

Received 17 January 2014; revised 17 February 2014; accepted 24 February 2014

Abstract

Diffuse interfaces appear with any Eulerian discontinuity capturing compressible flow solver. When dealing with multifluid and multimaterial computations, interfaces smearing results in serious difficulties to fulfil contact conditions, as spurious oscillations appear. To circumvent these difficulties, several approaches have been proposed. One of them relies on multiphase flow modelling of the numerically diffused zone and is based on extended hyperbolic systems with stiff mechanical relaxation (Saurel and Abgrall, 1999 [4] , Saurel et al., 2009 [6] ). This approach is very robust, accurate and flexible in the sense that many physical effects can be included: surface tension, phase transition, elastic-plastic materials, detonations, granular effects etc. It is also able to deal with dynamic appearance of interfaces. However it suffers of an important drawback when long time evolution is under interest as the interface becomes more and more diffused. The present paper addresses this issue and provides an efficient way to sharpen interfaces. A sharpening flow model is used to correct the solution after each time step. The sharpening process is based on a hyperbolic equation that produces a steady shock in finite time at the interface location. This equation is embedded in a “sharpening multiphase model” redistributing volume fractions, masses, momentum and energy in a consistent way. The method is conservative with respect to the masses, mixture momentum and mixture energy. It results in diffused interfaces sharpened in one or two mesh points. The method is validated on test problems having exact solutions.

Keywords:Material Interfaces; Multifluid; Multiphase; Multimaterial; Shocks; Non-Conservative; Hyperbolic Equations

1. Introduction

Numerical smearing of interfaces (or contact discontinuities) appears with any Eulerian discontinuity capturing compressible flow solver. In the context of the Euler equations, it results in at least four points in the capturing zone and the number of points increases during time evolution with most methods. When dealing with multifluid and multimaterial computations, diffuse interfaces result in serious difficulties to match interface conditions, as spurious oscillations appear. To circumvent these difficulties, several approaches have been proposed.

Lagrangian methods and Front Tracking schemes are aimed to remove artificial smearing but result in other difficulties.

Interface reconstruction methods (Hirt and Nichols, 1981 [1] , Youngs, 1989 [2] ) are aimed to restore sharp interfaces with the help of a volume fraction reconstruction algorithm but do not consider mass, momentum and energy redistribution. Moreover, these methods are more suitable for non-compressible fluids.

Level-set methods (Fedkiw et al., 1999 [3] ) are now very popular as they result in sharp interface representation and are quite easy to implement. They present however robustness and conservation issues.

Another option relies on multiphase flow modelling of the eight="37.3749995231628 " /> it reads,

. (4.9)

This equation expresses that the mixture density is going to sharpen during the resolution of System (4.3). However, as the mass transport has already been solved during the diffuse interface model resolution (3.1 or 3.7) no extra mass flux must appear. Equation (4.9) violates the local mass conservation. However, when integrated between the left and right side of the diffuse interface, mass is preserved. Indeed,

becomes,

.

The same remark holds for phase’s mass, momentum and energy:

(4.10)

and consequently for the phase’s total energies,

.

The method is thus conservative in the weak sense.

4.6. Summary

As the sharpening function gradient is not necessarily positive, the general 1D model reads:

(4.11)

Jump conditions are the same as previously with obvious modifications related to the sign.

5. Riemann Problem and Numerical Scheme

Flows computations proceed in two steps. First the flow model (System 3.7) is solved with Saurel et al. (2009) [6] method as summarized in Section 5.1. During this step, the interface becomes diffused. Thus, during the second step, the sharpening System (4.11) is used to remove most of the numerical diffusion at material interfaces. This system is hyperbolic and conservative. Conventional ingredients of hyperbolic conservation laws can thus be used.

5.1. Numerical Method for the Flow Model

Numerical resolution of the 6-equation augmented model (3.7) in the limit of stiff pressure relaxation has been addressed in Saurel et al. (2009) [6] . This method is particularly robust and quite simple. Computational examples and validations are given in Saurel et al. (2009) [6] and other references cited in introduction, where more sophisticated physical effects are considered.

The system to consider during numerical resolution thus involves 7 equations: System (3.7) and Equation (3.8).

Many approximate Riemann solvers can be considered to determine the cell boundary states, but the HLLC solver of Toro et al. (1994) [24] seems the simplest and the most efficient.

Saurel et al. (2009) [6] method can be summarized as follows:

• At each cell boundary solve the Riemann problem of System (3.7) with HLLC solver.

• Evolve all flow variables with a Godunov type scheme (or its higher order variants).

• Determine the relaxed pressure and especially the volume fraction.

• Compute the mixture pressure with Equation (3.4).

• Reset the internal energies with the computed pressure with the help of their respective EOS.

• Go to the first item for the next time step.

The main drawback of this method is related to the excessive numerical diffusion of interfaces, especially for long time evolutions. The sharpening System (4.11) is addressed in this aim. Its numerical resolution is detailed hereafter.

5.2. Approximate Riemann Solver

Consider a cell boundary separating a left state (L) and a right state (R). At each cell boundary an initial discontinuity is present and three waves are emerging:

• Two contact discontinuities with velocities:

,

• A shock wave propagating with speed:

.

For the sake of simplicity we still assume that. Moreover, the function that was previously a constant in the entire domain becomes now a local constant:

. (5.1)

It is important to consider as a local constant. Indeed, as the volume fraction may vary during the pressure relaxation step, volume fraction jumps may change while the sharpening function will not. This is the reason System (4.11) involves a volume fraction equation and another equation for the sharpening function.

For a given, corresponding to, three cases have to be considered:

*, *, *

Let us examine the first situation, the others being similar. The corresponding configuration is depicted in the Figure 4.

The shock speed is determined immediately from the initial data,

.

The sharpening function solution is unchanged compared to the scalar case of Section 2:

(5.2)

Consequently, the volume fractions are determined as:

System (4.4) can be cast under compact form as:

(5.3)

with.

The approximate Riemann problem solution for variables thus read:

, (5.4)

It is worth to mention that in this approximate solution, the solution for is uncoupled of those of and.

Therefore, when, whatever the situation is (or or), the approximate Riemann problem solution is always given by (5.2) and (5.4). The determination of the exact Riemann problem solution is possible, but not necessary for numerical purpose. For the Godunov method, only the flux along the axis is needed. Consequently the solution simplifies considerably and is summarised hereafter (below, superscript * represents the solution along the trajectory):

ComputeIf,

, , ,.

If,

, , , (5.5)

ComputeIf, ,.

If, ,.

The fluxes of System (4.11) are thus easily deduced. They are used in the Godunov (1959) method that reads:

(5.6)

with,

The fluxes are computed from the Riemann problem solution (5.4) and contain the appropriate sign function:

,.

The method is stable under conventional CFL restriction based on the fastest wave speed:

Theoretically, the interface is ideally sharpened when the steady state solution of (4.11) is obtained. For practical computations, a CFL number of 0.9 is used and a single time step is done.

5.3. Relaxation Step

In Section 4.5 we have shown that the sharpening method was not creating spurious pressure and velocity oscillations when the sharpening process was starting from initial data corresponding to mechanical equilibrium state. However, when a shock or an expansion wave interacts with an interface, the velocities and pressures are different on both sides of the interface at the discrete level. Thus, the mass redistribution achieved with the sharpening method is associated to momentum and energy redistribution that are going to produce non-equilibrium states on both sides of the interface. In order to restore mechanical equilibrium conditions, which consequence will be to fulfill the interface conditions of equal normal velocities and equal pressures, the following correction is done:

• The center of mass velocity is computed at each mesh point as,

and both phase velocity are reset to the centre of mass velocity,

.

• A correction is also needed for the internal energies and volume fractions. Indeed, as the mixture becomes out of pressure equilibrium, a pressure relaxation step is done, exactly in the same way as in Saurel et al. (2009) [6] . The volume fractions at pressure equilibrium are determined. Then the mixture pressure is computed with (3.4) and the internal energies are reset using the phase equation of state.

It means that System (4.11) has been complemented by velocity and pressure relaxation terms, in the same way as System (3.7) and that the limit solution when velocities and pressures are relaxed is used to reset the variables. Such correction is obviously conservative with respect to the masses, mixture momentum and mixture total energy. See for example Saurel and Abgrall (1999) [4] for the details in the context of a slightly different two-phase flow model.

6. One-Dimensional Illustrations

We now consider method validation on various test problems of increasing difficulty. A single pseudo time step with CFL = 0.9 is used. With such a choice the convergence is reached at each time step and the interface is computed in one point.

6.1. Advection of an Interface in a Uniform Pressure and Velocity Flow

A volume fraction discontinuity associated to a mixture density discontinuity is moving in a uniform pressure and velocity flow at 100 m/s. Initially the discontinuity is located at in a 1 m length tube. This discontinuity separates two nearly pure fluids:

• Liquid water on the left defined by, governed by the stiffened gas EOS parameters with parameters,

• Air on the right defined by and the EOS parameters and.

In the left chamber, the water volume fraction is set to and in the right chamber its value is

, with. The uniform pressure is set equal to.

A mesh involving 100 cells is used. The results are shown in the Figure 5 at time where the conventional method, without interface sharpening (left column) is compared to the version with sharpening correction (right column).

The sharpening algorithm clearly improves the results as the interface is handled in two points instead of 4 - 5 points. The pressure and velocity fields are maintained invariant.

6.2. Single Phase Shock Tube

The aim of this example is to show that the diffuse interface model with the sharpening algorithm is able to improve solutions of the Euler equations. Let’s consider a shock tube filled with air in both chambers. At initial time a pressure discontinuity is present at in a 1 m length tube. In the left chamber the pressure is set to 2 atm and to 1 atm in the right one. The density is initially uniform as well as the velocity that is zero in the entire domain.

In the Figure 6 the results obtained with a higher order extension of the Godunov scheme with MUSCL reconstruction and Superbee limiter are compared to those of the sharpening method with the two-phase model using two times the same ideal gas EOS with. Results are shown at time on a mesh

Figure 5. Advection of a volume fraction discontinuity in a uniform pressure and velocity flow. The diffuse interface model without interface sharpening (left column) is solved with the Saurel et al. (2009) [6] method extended to higher order with MUSCL reconstruction and Superbee limiter. The sharpening algorithm is employed for the results of the right column with CFL 0.9 with a single pseudo time step after each step of the diffuse interface solver..            

involving 100 cells.

A small oscillation due to the Superbee limiter is present on the velocity graph at shock front when the sharpening algorithm is used. These oscillations are absent with other limiters. Except regarding this discrepancy, the sharpening method improves the results.

6.3. Liquid Water-Air Shock Tube

A 1 m long shock tube containing two chambers separated by an interface at the location is considered. Each chamber contains a nearly pure fluid. The initial density of water is and the stiffened gas EOS parameters are and. The initial density of air is and EOS parameters are and. The left chamber contains a very small volume fraction of air and the initial pressure is set equal to 1 GPa. The right chamber contains the same fluids but the volume fractions are reversed. The initial pressure is set equal to 0.1 MPa. In both chambers the velocity is equal to zero. On Figure 7 the results are shown at time. A 1000 cells mesh is used.

A single time step with CFL = 0.9 is done with the sharpening solver after each hyperbolic time step of the diffuse interface model. The interface separating liquid and gas is clearly sharpened, as shown on the mixture density and volume fractions graphs.

6.4. Liquid Water-Air Shock Tube in Extreme Conditions

The same shock tube problem is solved, but initially, the left chamber pressure is set to 1 TPa (10 Mbars). The initial discontinuity is located at.This test illustrates the robustness and convergence of the new algorithm. The results are shown in the Figure 8 at time with a mesh using 1000 cells.

Again, the sharpening algorithm clearly improves the solution. An oscillation appears on the mixture density graph with the sharpening method, with similar behavior as Lagrangian and antidiffusion schemes (Kokh and Lagoutière, 2010) [14] . This effect is due to convergence errors that appear during the first time steps, when shock formation occurs. Our method is the only one able to deal with such pressure and density ratio with only one point in the interface.

6.5. Epoxy-Spinel Shock Tube

This test addresses method ability to deal with interfaces separating mixtures of materials. It highlights the need for an evolution equation for the sharpening variable, in addition to the volume fraction variable, both

Figure 6. Single phase shock tube. Comparison of a conventional Euler solver (MUSCL reconstruction with Superbee limiter) in symbols (left column) with the sharpening algorithm (with MUSCL reconstruction and Superbee limiter) in the single phase limit (square symbols) shown in the right column for the same variables. The volume fraction is shown for the sharpening algorithm.                                                                         

Figure 7. Water-air shock tube. The diffuse interface model without interface sharpening (left column) is solved with the Saurel et al. (2009) [6] method extended to higher order with MUSCL reconstruction with Superbee limiter. The solution with the sharpening algorithm and MUSCL reconstruction with Superbee limiter is shown on the right column.         

present in System (4.11).

A 1 m long shock tube containing two chambers separated by an interface at the location is considered. Each chamber contains a mixture of epoxy and spinel. The initial density of epoxy is and the stiffened gas EOS parameters are and. The initial density of spinel is and EOS parameters are and. The left chamber contains 70% of epoxy in volume and 30% spinel. The initial pressure is set equal to 2 GPa. The right chamber contains the same fluids but the volume fractions are reversed. The initial pressure is set equal to 0.1 MPa. In both chambers the velocity is equal to zero. On Figure 9 the results are shown at time. A 400 cells mesh is used.

In this example, the volume fraction and the sharpening function are not anymore linearly dependent. Their dependency is varying with time and space, showing the importance of Relation (5.1). The interface smearing is lowered with the sharp interface method.

7. More than Two Phases

The multiphase sharpening system can be rewritten as follows:

,  

,

,  

,

.

These systems are solved independently for an arbitrary number of fluids. The Godunov method presented in Section 5 results in the fulfilment of the constraints and.

Figure 8. Liquid water-air shock tube in extreme conditions. The diffuse interface model without interface sharpening (left column) is solved with the Saurel et al. (2009) [6] method and extended to higher order with MUSCL reconstruction and Superbee limiter. The solution with the sharpening algorithm and the same MUSCL reconstruction with Superbee limiter is shown on the right column.                                                                     

Figure 9. Epoxy-spinel shock tube problem. The diffuse interface model (square symbols) without interface sharpening solved with Saurel et al. (2009) [6] method is compared to the solution with the sharpening algorithm (cross symbols) and the exact solution. The same MUSCL reconstruction method with Van Leer limiter is used in both computations.       

8. Conclusions

An interface sharpening method has been presented. It is based on Harten (1977 [16] , 1978 [17] ) artificial compression method, and here it is extended to multiphase formulations of diffuse material interfaces (Kapila et al., 2001 [5] , Saurel et al., 2009 [6] ). The key point of the new method is that interface sharpening is achieved in a consistent way with masses, momentum and energies redistribution.

It preserves spurious pressure and velocity oscillations and handles interfaces in about two mesh points. It is also able to handle interfaces separating mixtures of materials.

Multidimensional extension is under examination.

Acknowledgements

This work was partially supported by A * MIDEX and ANR, the grant ANR-11-LABX-0092 and ANR-11- IDEX-0001-02.

References

  1. Hirt, C.W. and Nichols, B.D. (1981) Volume of Fluid (VOF) METHOD for the Dynamics of Free Boundaries. Journal of Computational Physics, 39, 201-255. http://dx.doi.org/10.1016/0021-9991(81)90145-5
  2. Youngs, D.L. (1989) Modelling Turbulent Mixing by Rayleigh-Taylor Instability. Physica D: Nonlinear Phenomena, 37, 270-287. http://dx.doi.org/10.1016/0167-2789(89)90135-8
  3. Fedkiw, R.P., Aslam, T., Merriman, B. and Osher, S. (1999) A Non-Oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method). Journal of Computational Physics, 152, 457-492. http://dx.doi.org/10.1006/jcph.1999.6236
  4. Saurel, R. and Abgrall, R. (1999) A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows. Journal of Computational Physics, 150, 425-467. http://dx.doi.org/10.1006/jcph.1999.6187
  5. Kapila, A.K., Menikoff, Bdzil, J.B., Son, S.F. and Stewart, D.S. (2001) Two-Phase Modeling of Deflagration to Detonation Transition in Granular Materials: Reduced Equations. Physics of Fluids, 13, 3002-3024. http://dx.doi.org/10.1063/1.1398042
  6. Saurel, R., Petitpas, F. and Berry, R.A. (2009) Simple and Efficient Relaxation for Interfaces Separating Compressible Fluids, Cavitating Flows and Shocks in Multiphase Mixtures. Journal of Computational Physics, 228, 1678-1712. http://dx.doi.org/10.1016/j.jcp.2008.11.002
  7. Perigaud, G. and Saurel, R. (2005) A Compressible Flow Model with Capillary Effects. Journal of Computational Physics, 209, 139-178. http://dx.doi.org/10.1016/j.jcp.2005.03.018
  8. Braconnier, B. and Nkonga, B. (2009) An All-Speed Relaxation Scheme for Interface Flows with Surface Tension. Journal of Computational Physics, 228, 5722-5739. http://dx.doi.org/10.1016/j.jcp.2009.04.046
  9. Saurel, R., Petitpas, F. and Abgrall, R. (2008) Modeling Phase Transition in Metastable Liquids. Application to Flashing and Cavitating Flows. Journal of Fluid Mechanics, 607, 313-350. http://dx.doi.org/10.1017/S0022112008002061
  10. Favrie, N., Gavrilyuk, S. and Saurel, R. (2009) Solid-Fluid Diffuse Interface Model in Cases of Extreme Deformations. Journal of Computational Physics, 228, 6037-6077. http://dx.doi.org/10.1016/j.jcp.2009.05.015
  11. Favrie, N. and Gavrilyuk, S. (2012) Diffuse Interface Model for Compressible Fluid-Compressible Elastic-Plastic Solid Interaction. Journal of Computational Physics, 231, 2695-2723. http://dx.doi.org/10.1016/j.jcp.2011.11.027
  12. Petitpas, F., Saurel, R., Franquet, E. and Chinnayya, A. (2009) Modelling Detonation Waves in Condensed Materials: Multiphase CJ Conditions and Multidimensional Computations. Shock Waves, 19, 377-401. http://dx.doi.org/10.1007/s00193-009-0217-7
  13. Saurel, R., Favrie, N., Petitpas, F., Lallemand, M.H. and Gavrilyuk, S. (2010) Modelling Dynamic and Irreversible Powder Compaction. Journal of Fluid Mechanics, 664, 348-396. http://dx.doi.org/10.1017/S0022112010003794
  14. Kokh, S. and Lagoutière, F. (2010) An Anti-Diffusive Numerical Scheme for the Simulation of Interfaces between Compressible Fluids by Means of a Five-Equation Model. Journal of Computational Physics, 229, 2773-2809. http://dx.doi.org/10.1016/j.jcp.2009.12.003
  15. So, K.K., Hu, X.Y. and Adams, N.A. (2012) Anti-Diffusion Interface Sharpening Technique for Two-Phase Compressible Flow Simulations. Journal of Computational Physics, 231, 4304-4323. http://dx.doi.org/10.1016/j.jcp.2012.02.013
  16. Harten, A. (1977) The Artificial Compression Method for Computation of Shocks and Contact Discontinuities. I. Single Conservation Laws. Communications on Pure and Applied Mathematics, 30, 611-638. http://dx.doi.org/10.1002/cpa.3160300506
  17. Harten, A. (1978) The Artificial Compression Method for Computation of Shocks and Contact Discontinuities: III. Self-Adjusting Hybrid Schemes. Mathematics of Computation, 32, 363-389.
  18. Shukla, R.K., Pantano, C. and Freund, J.B. (2010) An Interface Capturing Method for the Simulation of Multi-Phase Compressible Flows. Journal of Computational Physics, 229, 7411-7439. http://dx.doi.org/10.1016/j.jcp.2010.06.025
  19. Allaire, G., Clerc, S. and Kokh, S. (2002) A Five-Equation Model for the Simulation of Interfaces between Compressible Fluids. Journal of Computational Physics, 181, 577-616. http://dx.doi.org/10.1006/jcph.2002.7143
  20. Massoni, J., Saurel, R., Nkonga, B. and Abgrall, R. (2002) Proposition de Méthodes et Modèles Eulériens pour les Problèmes à Interfaces Entre Fluides Compressibles en Présence de Transfert de Chaleur: Some Models and Eulerian Methods for Interface Problems between Compressible Fluids with Heat Transfer. International Journal of Heat and Mass Transfer, 45, 1287-1307. http://dx.doi.org/10.1016/S0017-9310(01)00238-1
  21. Wood, A.B. (1930) A Textbook of Sound. G. Bell and Sons LTD, London.
  22. Saurel, R., Le Metayer, O., Massoni, J. and Gavrilyuk, S. (2007) Shock Jump Relations for Multiphase Mixtures with Stiff Mechanical Relaxation. Shock Waves, 16, 209-232. http://dx.doi.org/10.1007/s00193-006-0065-7
  23. Saurel, R., Gavrilyuk, S. and Renaud, F. (2003) A Multiphase Model with Internal Degree of Freedom: Application to Shock-Bubble Interaction. Journal of Fluid Mechanics, 495, 283-321. http://dx.doi.org/10.1017/S002211200300630X
  24. Toro, E.F., Spruce, M. and Speares, W. (1994) Restoration of the Contact Surface in the HLL Riemann Solver. Shock Waves, 4, 25-34. http://dx.doi.org/10.1007/BF01414629