## 1 Introduction

Stiff algebraic source terms, accounting for mechanical relaxation and phase transition in two-phase flow models of the Baer–Nunziato type

baernunziato ; pelanti ; spb , are one of the key difficulties in computing solutions to these systems of hyperbolic partial differential equations (PDE). Their accurate solution is relevant for the study of droplet dynamics with Baer–Nunziato models. These weakly compressible phenomena can be accurately described by the reduced models that assume instantaneous pressure and velocity equilibrium like the one forwarded by Kapila et al. kapila . Solving more general sets of equations like baernunziato ; pelanti ; spb in the stiff relaxation limit gives results that are similar to those obtained from the instantaneous equilibrium model, while allowing more modelling flexibility, since less physical assumptions have to be made.A simple computational strategy for dealing with stiff sources is the *splitting*
approach strang ; fractionalstep .
The procedure consists of two steps: at each timestep, first one solves the homogeneous part of the PDE

(1) |

for example with a path-conservative castro2006 ; pares2006 MUSCL–Hancock muscl method, obtaining a preliminary solution

and then one can use this state vector as initial condition for the Cauchy problem

(2) |

of which the solution will then yield the updated quantities at the new time level

. This way, the problem is reduced to the integration of a system of ordinary differential equations (ODE), and general-purpose ODE solvers or more specialised tools can be employed for this task.

It is often the case that the time scales associated with relaxations sources are much shorter than those given by the stability condition of the PDE scheme, thus one must be able to deal with source terms that are potentially stiff. In order to integrate stiff ODEs with conventional explicit solvers, one has to impose very severe restrictions on the maximum timestep size, and for this reason implicit methods are commonly preferred torobook . Unfortunately implicit solvers, are, on a per-timestep basis, much more expensive than explicit integrators, and they still might require variable sub-timestepping in order to avoid under-resolving complex transients in the solution.

In this work, we will develop a technique for constructing a solver for stiff finite-rate mechanical relaxation sources, specifically those encountered in models of the Baer–Nunziato type.

The proposed method overcomes the issues typical of explicit solvers with three concurrent strategies: first, the update formula is based on exponential integration expint1 ; expint2 , in order to mimic at the discrete level the behaviour of the differential equation; second, information at the new time level is taken into account by iteratively updating a linearisation of the ODE system, this is achieved without resorting to a fully implicit method like those introduced in butcher1964 , and for which one would need to solve a system of nonlinear algebraic equations at each timestep ; third, the method incorporates a simple and effective adaptive timestepping criterion, which is crucial for capturing abrupt changes in the state variables and dealing with the different time scales that characterise the equations under investigation.

## 2 Model equations

We are interested in the solution of two-phase flow models of the Baer–Nunziato family, which can be written in the general form (1), with a vector of conserved variables defined as

(3) |

a conservative flux and a non-conservative term written as

(4) |

and a source term vector written as

(5) |

Here we indicate with and the volume fractions of the first phase and of the second phase respectively, with and the phase densities, and indicate the velocity vectors, and are the partial energy densities. The pressure fields are denoted with and , and the interface pressure and velocity are named and . Finally, the parameters and control the time scales for friction and pressure relaxation kinetics respectively.

In the following, we will study the system of ordinary differential equations arising from the source term (5) only, that is, the one constructed as given in equation (2) and specifically its one-dimensional simplification in terms of the primitive variables , with an initial condition . Since no source is present in the mass conservation equations, they have a trivial solution, that is, and remain constant in time; for compactness, these quantities will be included in our analysis as constant parameters, rather than as variables of the ODE system.

The one-dimensional ODE system is written as

(6) | |||

(7) | |||

(8) | |||

(9) | |||

(10) |

The choices for interface pressure and velocity are and . Finally, one can verify that, using the stiffened-gas equation of state for both phases, we have , , , and .

## 3 Description of the numerical method

The methodology is described in the following with reference to a generic nonlinear first order Cauchy problem

(11) |

for which the ODE can be linearised about a given state and time as

(12) |

Here we defined the Jacobian matrix of the source and analogously the source vector evaluated at the linearisation state is . We then introduce the vector

(13) |

which will be used as an indicator for the adaptive timestepping algorithm and may be constructed for example listing all of the components of the matrix together with all the components of the vector and the state , or only with a selection of these variables, or any other relevant combination of the listed variables, that is, any group indicative of changes in the nature or the magnitude of the linearised source operator.

It is then necessary to compute an accurate analytical solution of the non-homogeneous linear Cauchy problem

(14) |

We will denote the analytical solution of the IVP (14) as . As for , the semicolon separates the variable on which and continuously depend ( or ) from the parameters used in the construction of the operators. The state vector at a generic time level is written as , the variable timestep size is .

### 3.1 Timestepping

Marching from a start time to an end time is carried out as follows. First, an initial timestep size is chosen, then, at each time iteration, the state at the new time level is computed by means of the iterative procedure described below. The iterative procedure will terminate by computing a value for , together with a new timestep size

based on an estimator which is embedded in the iterative solution algorithm. There is also the possibility that, due to the timestep size

being too large, the value of be flagged as not acceptable. In this case, the procedure will return a new shorter timestep size for the current timestep and a new attempt at the solution for will be carried out. Specifically, in practice we choose the new timestep size to be half of the one used in the previous attempt.### 3.2 Iterative computation of the timestep solution

At each iteration (denoted by the superscript ) we define an average state vector to be formally associated with an intermediate time level . For the first iteration we need a guess value for , with the simplest choice being . Then the coefficients are computed as

(15) |

In a joint way, one can build the affine source operator

(16) |

Then one can solve analytically

(17) |

by computing

(18) |

It is then checked that the state vector be physically admissible: in our case this means verifying that internal energy of each phase be positive and that the volume fraction be bounded between 0 and 1. Also one can check for absence of floating-point exceptions. Additionally, one must evaluate

(19) |

This vector of coefficients will not be employed for the construction of an affine source operator , but only for checking the validity of the solution obtained from the approximate problem (17) by comparing the coefficients vector to , as well as comparing the coefficients used in the middle-point affine operator for the initial coefficients . At the end of the iterative procedure, one will set , so that this will be the new reference vector of coefficients for the next timestep. The convergence criterion for stopping the iterations is implemented by computing

(20) |

and checking if , with and given tolerances, or if the iteration count has reached a fixed maximum value . Note that in principle any norm may be used to compute the error metric given in equation (20), as this is just a measure of the degree to which was corrected in the current iteration. Moreover, we found convenient to limit the maximum number of iterations allowed, and specifically here we set , but stricter bounds can be used. For safety, we decide to flag the state vector as not admissible, as if a floating-point exception had been triggered, whenever the iterative procedure terminates by reaching the maximum iteration count.

After the convergence has been obtained, in order to test if the IVP (11) is well approximated by its linearised version (17), we compute

(21) | |||

(22) |

and we verify if The user should specify a tolerance as well as the floor value , which is used in order to prevent that excessive precision requirements be imposed in those situations when all the coefficients are so small than even large relative variations expressed by equations (21) and (22) do not affect the solution in a significant manner. If we confirm the state vector at the new time level to be and a new timestep size is computed as

(23) |

otherwise the solution of the IVP (17) is attempted again with a reduced timestep size, specifically one that is obtained by halving the timestep used in the current attempt. The same happens if at any time the admissibility test on fails.

### 3.3 Analytical solution of the linearised problem

The general solution to an initial value problem like (17) can be written as

(24) |

Note that, in addition to evaluating the matrix exponential , one must also compute the inverse Jacobian matrix . Computation of matrix exponentials can be carried out rather robustly in double precision arithmetic with the aid of the algorithms of Higham matexp2005 and Al-Mohy and Higham matexp2009 ; matexp2011 , while inversion of the Jacobian matrix can be an arbitrarily ill-conditioned problem, to be carefully treated or avoided if possible.

For this reason we propose the following strategy for choosing a more suitable linearisation and computing analytical solutions of the linearised problem for the ODE system (6)–(10). First, it is easy to see that the velocity sub-system (equations for and ) can be fully decoupled from the other equations, as the partial densities and remain constant in the relaxation step.

Then the solution of the velocity sub-system can be immediately obtained as

(25) | |||

(26) |

with .

In a second step, the pressure sub-system (8)–(9) is linearised as

(27) | |||

(28) |

where and are constant coefficients directly obtained from equations (8)–(9). This way, at the cost of suppressing the dependence on in the Jacobian of the pressure sub-system, the homogeneous part of equations (27)–(28) has the same simple structure found in the velocity sub-system, with the addition of a non-homogeneous term, which is known, as and already have been computed. The solution can again be evaluated using standard scalar exponential functions, which are fast and robust, compared to matrix exponentials and especially so, because one no longer needs to perform the inversion of the Jacobian matrix of the full system.

Finally, the solution to equation (10) can be integrated analytically from the expressions of and . Full coupling of the system is restored in the successive iterations by recomputing the constant coefficients and using an updated midpoint value for . See Figure 1 for a graphical description of the proposed simplified solution structure.

## 4 Test problems

We provide validation of the proposed method first by computing solutions to the ODE system (6)–(10) and comparing the results with a reference solution obtained from a sixth order, fully implicit, Runge–Kutta–Gauss–Legendre method butcher1964 (labeled RKGL3) employing adaptive timestepping (test problems A1 and A2, Figures 2 and 3). Furthermore, test problem A1 is employed also for carrying out a convergence study of the scheme (Figure 4), showing that second order convergence is easily achieved. The initial data for the ODE tests are, for test A1,

(29) |

while for test A2,

(30) |

The parametric data are, for test A1,

(31) | ||||

and for test A2,

(32) | ||||

Then, we show an application of the method in the solution of the mixture-energy-consistent formulation of the six-equation reduced Baer–Nunziato model forwarded in pelanti . For these simulations the interface pressure is computed as

(33) |

The first two shock-tube problems (from cavitation_shocktube ; pelanti ), show that the method is able to deal with very stiff () sources, and in particular in Figure 5 (RP1, a liquid-vapour dodecane shock tube featuring a strong right-moving shockwave) we show mesh convergence of the solution by comparing two runs, both employing the HLLEM Riemann solver proposed in hllem , on two different meshes consisting of 2 000 uniform control volumes and 20 000 control volumes respectively, with a computational domain delimited by . In Figure 6 (RP2, two diverging rarefaction waves in liquid water) we then show that, with very stiff relaxation (), the solution matches the one computed by solving directly the five-equation instantaneous equilibrium model kapila , again using a mesh consisting of 2 000 uniform cells for the six-equation model and a mesh of 20 000 uniform cells for the reference solution, and in particular, rarefaction waves propagate at the same speed for both models. All tests are run using a second order path-conservative MUSCL-Hancock scheme with . The first Riemann Problem (RP1) is set up with uniform liquid and vapour densities and , uniform velocity , a jump in pressure given by , , almost pure liquid on the left side of the initial discontinuity (), and almost pure vapor on the right side (). The discontinuity is initially found at , and the end time is . The parameters of the stiffened gas EOS are , , , .

The second Riemann Problem (RP2) is initialised with constant liquid and vapour densities , , constant pressure , constant liquid volume fraction , and a jump in velocity (initially located at ) such that and . The final time is and for this test the parameters of the equation of state , , , .

Finally, in Figure 7 we show the behaviour of the solution of a third Riemann problem (RP3) with several different values of the pressure relaxation parameter (ranging from to ), highlighting the vast range of solution structures that can be obtained not only with stiff relaxation (the pressure profiles and coincide) or in total absence of it, but also with finite values of the relaxation time scale. For RP3, the initial data on the left are

(34) | ||||

while on the right one has

(35) | ||||

The initial jump is located at , the domain is and the final time is . The parameters of the stiffened gas EOS are , , , .

## 5 Conclusions

We presented a technique for integrating ordinary differential equations associated with stiff relaxation sources and promising results have been shown for a set of test problems. The method can efficiently resolve very abrupt variations in the solution and adapt to multiple timescales. A key feature of the algorithm is that it can avoid delicate linear algebra operations entirely, thus improving the robustness of the scheme. Future applications will include liquid-gas and liquid-solid phase transition, strain relaxation for nonlinear elasticity gpr and the computation of material failure in elasto-plastic and brittle solids.

###### Acknowledgements.

The authors of this work were supported by the German Research Foundation (DFG) through the project GRK 2160/1 “Droplet Interaction Technologies”. This is a pre-print of the following work: G. Lamanna, S. Tonini, G.E. Cossali and B. Weigand (Eds.), “Dropet Interaction and Spray Processes”, 2020, Springer, Heidelberg, Berlin. Reproduced with permission of Springer Nature Switzerland AG. DOI: 10.1007/978-3-030-33338-6.## Bibliography

- [1] A. H. Al-Mohy and N. J. Higham. A new scaling and squaring algorithm for the matrix exponential. SIAM Journal on Matrix Analysis and Applications, 31(3):970–989, 2009.
- [2] A. H. Al-Mohy and N. J. Higham. Computing the action of the matrix exponential, with an application to exponential integrators. SIAM Journal on Scientific Computing, 33(2):488–511, 2011.
- [3] M. R. Baer and J. W. Nunziato. A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials. International Journal of Multiphase Flow, 12(6):861–889, 1986.
- [4] J. C. Butcher. Implicit Runge–Kutta processes. Mathematics of Computation, 18(85):50–64, 1964.
- [5] M. Castro, J. M. Gallardo, and C. Parés. High order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. Applications to shallow-water systems. Mathematics of Computation, 75(255):1103–1134, 2006.
- [6] J. Certaine. The solution of ordinary differential equations with large time constants. Mathematical Methods for Digital Computers, pages 128–132, 1960.
- [7] M. Dumbser and D. S. Balsara. A new efficient formulation of the HLLEM Riemann solver for general conservative and non-conservative hyperbolic systems. Journal of Computational Physics, 304(C):275–319, 2016.
- [8] M. Dumbser, I. Peshkov, E. Romenski, and O. Zanotti. High order ader schemes for a unified first order hyperbolic formulation of continuum mechanics: Viscous heat-conducting fluids and elastic solids. Journal of Computational Physics, 314:824–862, 2016.
- [9] N. J. Higham. The scaling and squaring method for the matrix exponential revisited. SIAM Journal on Matrix Analysis and Applications, 26(4):1179–1193, 2005.
- [10] A. K. Kapila, R. Menikoff, J. B. Bdzil, S. F. Son, and D. S. Stewart. Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations. Physics of Fluids, 13(10):3002–3024, 2001.
- [11] C. Parès. Numerical methods for nonconservative hyperbolic systems: a theoretical framework. SIAM Journal on Numerical Analysis, 44(1):300–321, 2006.
- [12] M. Pelanti and K.-M. Shyue. A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves. Journal of Computational Physics, 259:331–357, 2014.
- [13] D. A. Pope. An exponential method of numerical integration of ordinary differential equations. Communications of the ACM, 6(8):491–493, 1963.
- [14] R. Saurel, F. Petitpas, and R. A. Berry. Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures. Journal of Computational Physics, 228(5):1678–1712, 2009.
- [15] Richard Saurel, Fabien Petitpas, and Remi Abgrall. Modelling phase transition in metastable liquids: application to cavitating and flashing flows. Journal of Fluid Mechanics, 607:313––350, 2008.
- [16] G. Strang. On the Construction and Comparison of Difference Schemes. SIAM Journal on Numerical Analysis, 5:506–517, 1968.
- [17] E. F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction, Third edition. Springer-Verlag, Berlin, 2009.
- [18] B. van Leer. Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method. Journal of Computational Physics, 32(1):101–136, 1979.
- [19] N. N. Yanenko. The method of fractional steps. Springer, 1971.

Comments

There are no comments yet.