Model.xml

Fluid mechanics models used in the computation are specified in the model.xml file. It is mandatory located in the folder of the current case. A typical form of this file is:

<?xml version = "1.0" encoding = "UTF-8" standalone = "yes"?>
<model>
  <flowModel name="PressureVelocityEq" numberPhases="2" alphaNull="false"/>
  <EOS name="IG_air.xml"/>
  <EOS name="SG_water.xml"/>
</model>

Flow model

<flowModel name="PressureVelocityEq" numberPhases="2" numberTransports="1" alphaNull="false"/>

The <flowModel> markup is mandatory to specify the mathematical model to solve during the computation. This markup may contain the following attributes:

  • name: Name of the mathematical flow model. This attribute can take the values: Euler [Eul57], PressureVelocityEq (previously named Kapila) [SPB09][KMB+01], VelocityEq (previously named multiP) [SCazeP+23], VelocityEqTotEnergy, NonEqDEM [AS03], NonEq [FS15], TemperaturePressureVelocityEq (previously named ThermalEq) [LMSN14], EulerHomogeneous, ExtendedNavierStokesKorteweg1D [DFL26].

  • numberPhases: Integer number corresponding to the number of phases present in the simulation. The total amount of equations is related to this number. This attribute is not necessary when the values of name are Euler or EulerHomogeneous.

  • numberTransports: This attribute is optional and is set to 0 as default. It can be used to add specific variables advected in the flow (color functions).

  • alphaNull: For PressureVelocityEq, VelocityEq, VelocityEqTotEnergy and NonEqDEM models, the volume fraction can either be null (true) or not (false) and this choice is determined by this parameter. Default value is false.

Remark:

If EulerHomogeneous model is chosen, two other attributes are to be used: liquid and vapor to specify the number corresponding to the liquid phase and the vapor phase. They are phase 0 for the first and 1 for the second.

<flowModel name="EulerHomogeneous" liquid="0" vapor="1"/>
<EOS name="SG_waterLiq.xml"/>
<EOS name="IG_waterVap.xml"/>

If ExtendedNavierStokesKorteweg1D model is chosen, other attributes are needed, in accordance with the work of Didierlaurent et al. [DFL26]. From this generic model, one can choose different equations of state to obtain particular models lying in different fields: Defocusing Nonlinear Schrödinger equation [DFG19], dynamics of a bistable tape spring [BFL20], or the dynamics of phase transition by considering the van der Waals EOS [DFL26].

<flowModel name="ExtendedNavierStokesKorteweg1D" lambda="84000" beta="2.15e-14" tau="1." eps="1."/>
<EOS name="VDW_dodecane.xml"/>

If NonEqDEM model is selected, two additional attributes may be specified. RiemannSolver to determine which Riemann solver will be used during simulation, either HLLC or Exact, and contactAreas to specify which contact area algorithm will be used, either Symmetric or NonSymmetric.

<flowModel name="NonEqDEM" numberPhases="2" alphaNull="false" RiemannSolver="HLLC" contactAreas="Symmetric"/>
<flowModel name="NonEqDEM" numberPhases="2" alphaNull="false" RiemannSolver="Exact" contactAreas="NonSymmetric"/>

If NonEq model is selected, an additional attribute must be set: solver to specify the solver of the Riemann problem. The HLLAC (HLLC-type) solver presented in [LDGH16], the HLLC4waves solver presented in [FS15], the HLL- or Rusanov-type solvers presented in [SA99] (with different closure laws) are implemented in ECOGEN to solve the total-non-equilibrium model [SA99]. They are referenced in Section Non-equilibrium model.

<flowModel name="NonEq" numberPhases="2" alphaNull="false" solver="HLLAC"/>
<flowModel name="NonEq" numberPhases="2" alphaNull="false" solver="HLLC4waves"/>
<flowModel name="NonEq" numberPhases="2" alphaNull="false" solver="HLL"/>

If PressureVelocityEq model is selected, an additional attribute might be set: threshold. If the volume fraction of a phase is superior to 1 - threshold, the pressure-relaxation procedure won’t be applied, otherwise, it will be applied. If this attribute is not specified, the default value is 1.e-5.

<flowModel name="PressureVelocityEq" numberPhases="2" numberTransports="0" alphaNull="false" threshold="1.e-6"/>

Equations of state (EOS)

<EOS name="IG_air.xml"/>

The model.xml input file must contain as many <EOS> markups as number of phases specified in the Flow model markup. Each phase is described thanks to relations and parameters (see Section Materials for more details). The values of these parameters are specified in a separate file: the name attribute contains the name of this file which must be placed in the ECOGEN/libEOS/ folder. Some fluid files are already present in the ECOGEN package.

Global accuracy order of the numerical scheme

By default, the numerical scheme of the fluid-flow solver is first order in time and space. ECOGEN can use a second-order scheme in time and space using the MUSCL method with TVD limiters (see [Tor13] for example for an overview of the method). Depending on whether the mesh is Cartesian (with/without AMR) or unstructured, the MUSCL method used and its parameters are different.

Cartesian mesh (with/without AMR)

For a cartesian mesh (with/without AMR), second-order scheme is based on a TVD slope limiter; see [SPLMD20] for details. In this case, the optional markup <secondOrder> can be inserted in the model.xml input file as in the following example:

<secondOrder>
  <globalLimiter>minmod</globalLimiter>
  <interfaceLimiter>superbee</interfaceLimiter>                            <!-- optional node -->
  <globalVolumeFractionLimiter>minmod</globalVolumeFractionLimiter>        <!-- optional node -->
  <interfaceVolumeFractionLimiter>thinc</interfaceVolumeFractionLimiter>   <!-- optional node -->
</secondOrder>

The <secondOrder> markup must contain the node <globalLimiter>. The other nodes are optional. The slope limiters available in ECOGEN are the following: minmod [Roe86], vanleer [VL74], vanalbada [VAVLR97], mc [VL77], superbee [Roe86], koren [VK93] and thinc [SX14] (only for volume fraction at the interface location). Markup significations are:

  • <globalLimiter>: Applied everywhere and on all variables unless it is overwritten by the following optional limiters.

  • <interfaceLimiter>: Applied on all variables but only at the interface location. By default is equal to the global limiter.

  • <globalVolumeFractionLimiter>: Applied everywhere but only on the volume-fraction and transport equations (THINC is only applied on the volume fraction) unless it is overwritten by the interface volume-fraction limiter. By default is equal to the global limiter.

  • <interfaceVolumeFractionLimiter>: Applied only at the interface location and on the volume-fraction and transport equations (THINC is only applied on the volume fraction). By default is equal to the interface limiter.

Unstructured mesh

For an unstructured mesh, second-order scheme is based on a TVD gradient limitation; see [CSN17] for details. In this case, the chosen limiter is a global limiter applied everywhere and on all variables. Currently the gradient limiters are the following: minmod [Roe86] and superbee [Roe86]. The optional markup <secondOrder> can be inserted in the model.xml input file as in the following example:

<secondOrder>
  <globalLimiter>minmod</globalLimiter>
</secondOrder>

Important

Since the MUSCL method relies on gradients computation, the gradient method chosen should be compatible with unstructured meshes (see the markup <gradient>).

Gradient method

In ECOGEN, to compute gradients, it is possible to use:

  • finite-difference-like gradient on cartesian mesh (without AMR) (<method>finite-difference</method>)

  • Green-Gauss gradient on cartesian and unstructured mesh (<method>green-gauss</method>).

By default (without the XML markup <gradient>), the gradients will be computed using the Green-Gauss scheme. To define explicitly the gradient method, one can use:

<gradient>
  <method>finite-difference</method> <!-- For finite difference like gradient scheme -->
</gradient>

Note

The Green-Gauss method for the gradient computation must be selected when:

  • AMR is used,

  • periodic boundary conditions are used,

  • unstructured meshes, in case second-order scheme is set and/or additional physics are used (see Section Additional physics).

Advected additional variables

<transport name="color"/>

The model.xml input file must contain as many <transport> markups as number of transports specified in the Flow model markup. Each transported variable is described by its name. The default number of advected variable is 0.

Relaxation procedures

An additional markup <relaxation> may be used to impose some specific equilibrium between the phases depending on the flow model used. The attribute type specifies the type of equilibrium:

  • U: A velocity equilibrium is imposed at every location of the flow. The attribute speed can be added to specify the speed at which the relaxation operates (only infinite [SA99] for now). Default is infinite.

<relaxation type="U" speed="infinite"/>

  • P: A pressure equilibrium is imposed at every location of the flow. The attribute speed can be added to specify the speed at which the relaxation operates (infinite [SPB09] or finite [BPS23, SCazeP+23]). Default is infinite. In this case, an additional attribute might be set: threshold to determine when the pressure-relaxation procedure is bypassed. If a volume fraction exceeds 1 - threshold, the procedure is not applied; otherwise, it is executed. If a finite relaxation is chosen, the rate and solver have to be specified. The rate is a real number while the solver can either be Euler [SCazeP+23] or LSODA [Hin83, Pet83].

<relaxation type="P" speed="infinite" threshold="1.e-6"/>

<relaxation type="P" speed="finite" rate="1.e-5" solver="Euler"/>

  • PT: Both pressure and thermal equilibrium are imposed at every location of the flow. It does not require additional attributes.

<relaxation type="PT"/>

  • PTMu: A thermodynamical equilibrium is imposed at every location of the flow. It must be associated with the node <dataPTMu> whose attributes are liquid and vapor to specify the name of the EOS of the liquid and the vapor phase. Hereafter the complete node when PTMu is used:

<relaxation type="PTMu">
  <dataPTMu liquid="SG_waterLiq.xml" vapor="IG_waterVap.xml"/>
</relaxation>

Source terms

The additional <sourceTerms> markup can be used to numerically integrate some source terms in the equations. The attribute type selects the source term:

  • heating: Related to a thermal energy heating/cooling. This attribute requires the <dataHeating> node with the attribute volumeHeatPower; a real number corresponding to the power by volume unit added to the flow (unit: W/m3 (SI)).

<sourceTerms type="heating" order="EULER">
  <dataHeating volumeHeatPower="1.e6"/>
</sourceTerms>

  • gravity: If the gravity is considered. The node <dataGravity> must be present with the attributes x, y and z giving the coordinates for the gravity acceleration vector in real numbers (unit: m/s2 (SI)).

<sourceTerms type="gravity" order="EULER">
  <gravity x="0." y="-9.81" z="0."/>
</sourceTerms>

  • MRF: For a simulation in a moving reference frame. Allow to compute solution in a rotating frame [CazePD+24]. The node <omega> requires the attributes x, y and z giving the coordinates for the rotating vector in real numbers (unit: rad/s (SI)). The node <timeToOmega> is optional and allow to specify a progressing acceleration (linear) to the final rotating velocity (requires the attribute tf for acceleration time).

<sourceTerms type="MRF" order="RK4">
  <omega x="0." y="0." z="1."/>
  <timeToOmega tf="1.e-3"/>  <!-- Optional: If activated, the angular velocity increase linearly to omega in during tf -->
</sourceTerms>

  • Acoustic waves: Acoustic waves can be modelled with sinusoidal or Gaussian profils. This follows the work of [MC17] for single-phase flows and has been extended for multiphase system of equations. The node typeAcousticWave requires you to mention either a sinusoidal or a Gaussian type of pulses. For both, nodes sourceLocation and pulseDirection are required and have the attributes x, y and z giving the coordinates of the source and the direction of the pulses in real numbers (unit: m (SI) and none), respectively. For pulseDirection, the sum of the attributes must be equal to 1. Additionally, the pressureAmplitude (unit: Pa (SI)) and delayedTime (unit: s (SI)) of the pulses are required for both. The sinusoidal pulses require the numberOfPulse (can be a real number) and frequency (unit: Hz (SI)), while the Gaussian pulses require the pulseWidth (unit: m (SI)). XML blocks can be written for sinusoidal or Gaussian pulses as:

<sourceTerms type="acousticWave" order="EULER">
  <typeAcousticWave type="planeSinusoidalPulse">
    <dataPlaneSinusoidalPulse pressureAmplitude="10000." delayedTime="5.e-7" numberOfPulse ="4." frequency="3.e5">
      <sourceLocation x="10.e-3" y="10.e-3" z="0.5"/>
      <pulseDirection x="1." y="0." z="0."/>
    </dataPlaneSinusoidalPulse>
  </typeAcousticWave>
</sourceTerms>

<sourceTerms type="acousticWave" order="EULER">
  <typeAcousticWave type="planeGaussianPulse">
    <dataPlaneGaussianPulse pressureAmplitude="10." delayedTime="20.e-6" pulseWidth="5.e-6">
      <sourceLocation x="10.03e-3" y="0." z="0."/> <!-- Must be placed at: cell_center + 1/4 * cell_size_in_wanted_direction -->
      <pulseDirection x="1." y="0." z="0."/>
    </dataPlaneGaussianPulse>
  </typeAcousticWave>
</sourceTerms>

Source-term integration can be done up to 4th order. Available integration schemes are:

  • EULER (1st order)

  • Runge–Kutta 2 RK2 (2nd order)

  • Runge–Kutta 4 RK4 (4th order)

Scheme selection is done through order attribute.

Symmetry terms

Both cylindrical (2D) and spherical (1D) symmetries are implemented. The additional <symmetryTerm> markup can be used. It requires the attribute type that can take the value cylindrical or spherical. It also requires an additional node to specify the radial axis:

  • Cylindrical:

<symmetryTerm type="cylindrical">
  <dataSymCyl radialAxis="X"/>
</symmetryTerm>

  • Spherical:

<symmetryTerm type="spherical">
  <dataSymSpher radialAxis="X"/>
</symmetryTerm>

Additional physics

Depending on the model chosen in section Flow model, additional physical effects can be added. This is the case for surface tension, viscosity, turbulence and conductive heat transfers. These additional physical effects are obtained thanks to the additional markup additionalPhysics with the attribute type that can take different value according to the chosen effect.

Surface tension

This physical effect is obtained by using the type surface tension and its treatment is detailed in [SPD+17]. Then it requires the node dataSurfaceTension with following attributes:

  • transport: This is the name of advected variable used as color function for surface-tension terms. This advected variable has been precised in the section Global accuracy order of the numerical scheme. The name should be the same.

  • sigma: A real number for the surface-tension coefficient (unit: N/m (SI)).

<additionalPhysics type="surfaceTension" >
  <dataSurfaceTension transport="color" sigma="72.e-3"/>
</additionalPhysics>

Viscosity

This physical effect is obtained by using the type viscosity and it needs the attribute mu of the <physicalParameters> to be set in the EOS files (Materials).

<additionalPhysics type="viscosity"/>

Turbulence: static Smagorinsky

This physical effect is obtained by using the type turbulence. One can choose the turbulence model by changing the <model> node. The current implemented LES model is the static Smagorinsky [Sma63]. It needs the attributes C_s, C_I and Pr_T of the <smagorinsky> node to be set in the model file (Model.xml). The value of C_s depends on the problem, but here are some common values: C_s=0.16 [PMPC00], C_s=0.0085 [KYW+16]. Additionally, it needs the attributes mu and lambda of the <physicalParameters> to be set in the EOS files (Materials) to include the molecular viscosity and conductivity effects. Finally, one can set the attribute ducros to true/false in order to activate/deactivate the Ducros et al. flow sensor [DFN+99]. This feature is currently available for the Euler flow model.

<additionalPhysics type="turbulence">
        <turbulence model = "smagorinsky">
                <smagorinsky C_s="0.12" C_I="0.09" Pr_T="0.71" ducros="false"/>
        </turbulence>
</additionalPhysics>

Warning

When enabling the turbulence model, one must not enable viscosity and conductivity, as they are already taken into account. Setting the attributes mu and lambda of the <physicalParameters> is enough.

Strang splitting (only for Sources and Additional physics)

For a system of balance laws of this form:

\begin{align} &\cfrac{\partial \mathbf{Q}}{\partial t} + \textrm{div}( \mathbf{F}(\mathbf{Q}) ) + \mathbf{G}(\mathbf{Q}) \textrm{div}(\mathbf{u}) = \mathbf{S}(\mathbf{Q}), \end{align}

it is common to split the hyperbolic part and the source part:

\begin{array}{*{20}{l}} \mathscr{H}(t) : \quad \cfrac{\partial \mathbf{Q}}{\partial t} + \textrm{div}( \mathbf{F}(\mathbf{Q}) ) + \mathbf{G}(\mathbf{Q}) \textrm{div}(\mathbf{u}) = 0,\\ \mathscr{S}(t) : \quad \cfrac{\partial \mathbf{Q}}{\partial t} = \mathbf{S}(\mathbf{Q}), \end{array}

where \(\mathscr{H}(t)\) is the hyperbolic operator and \(\mathscr{S}(t)\) the source operator. Each operator is discretized by an appropriate numerical method of any order of accuracy. Such splitting procedure leads to a splitting error independant of the numerical discretization of the operators. To obtain the solution at the next step, one have to combine the numerical solutions of each operator in some ordered sequence. Depending on the sequence, the splitting error (and thus the accuracy) change. For a first order accurate in time, the Godunov splitting can be used:

\begin{align} &\mathbf{Q}^{n+1} = \mathscr{S}(\Delta t) \, \mathscr{H}(\Delta t) \, \mathbf{Q}^{n} \end{align}

For a second order accurate in time, the Strang splitting can be considered [Str68] [LeVeque02]:

\begin{align} &\mathbf{Q}^{n+1} = \mathscr{S} \left( \cfrac{\Delta t}{2} \right) \, \mathscr{H}(\Delta t) \, \mathscr{S} \left( \cfrac{\Delta t}{2} \right) \, \mathbf{Q}^{n} \end{align}

If the nondiscretized operators \(\mathscr{H}\) and \(\mathscr{S}\) commute \((\mathscr{H} \mathscr{S} - \mathscr{S} \mathscr{F} = 0)\), then the Godunov splitting is also second order accurate in time. By default, the Godunov splitting is used in ECOGEN. To use the Strang splitting, put the following line in model.xml:

<splitting type="strang"/>

Others

In dev…

Low-Mach preconditioning

In case of low-speed flows, a low-Mach preconditioning can be applied to the selected flow model. The method is detailed in [LMNS13]. Currently this option is compatible with the flow models: Euler, UEq and PUEq.

<lowMach state="true"/>

Note that the minimum reference Mach number can be chosen if required (default is 0.01). When the local Mach number of the flow is lower than this threshold, it is replaced by the threshold value in order to keep a reasonable computation time.

<lowMach state="true" machRefMin="5.e-2"/>

For an example of use, see the test case presented in the Section Low-Mach subsonic flow.

1D geometry with smooth cross section variation

In case of a 2D geometry with smooth cross section variation, it is possible to use a 1D geometry and reproduce 2D effects with this attribute as detailed in the first appendix of [LMNS13]. This option is only compatible with unstructured meshes defined with flow direction along the X-axis. Boundary condition contributions in other directions than the X-axis should use nullFlux boundary condition (see InitialConditions.xml). Note that this option is available for the flow models: Euler, UEq, PUEq and EulerHomogeneous. For an example of use, see the test case presented in the Section Low-Mach subsonic flow.

<geometry smoothCrossSection1d="true"/>