Multi-material DG hydrodynamics

Contents

Inciter supports multiple hydrodynamics schemes. This page describes governing equations and the DG method for multi-material flows.

Governing equations for multi-material hydrodynamics

The stiff velocity equilibrium multi-material hydrodynamic equations are considered in this work. In this system, although all the materials are advected with the same velocity $u_j$ , they possess different pressures and internal energies. Finite pressure relaxation source terms are added to the material volume fraction and material total energy equations based on the material bulk-modulii. The Eulerian form of this multi-material system of equations is,

\[ \begin{split} \frac{\partial \boldsymbol{U}}{\partial t} + \frac{\partial \boldsymbol{F}_j}{\partial x_j} + \boldsymbol{D} = \boldsymbol{S}, \end{split} \]

where,

\[ \begin{split} & \boldsymbol{U} = \begin{bmatrix} \alpha_k \\ \alpha_k \rho_k \\ \overline{\rho} u_i \\ \alpha_k \rho_k E_k \end{bmatrix}, \quad \boldsymbol{F}_j = \begin{bmatrix} 0 \\ \alpha_k \rho_k u_j \\ \overline{\rho} u_i u_j + \overline{p} \delta_{ij} \\ \alpha_k \rho_k H_k u_j \end{bmatrix}, \\ & \boldsymbol{D} = \begin{bmatrix} u_j \frac{\partial \alpha_k}{\partial x_j} \\ 0 \\ 0 \\ Y_k u_j \frac{\partial \overline{p}}{\partial x_j} - u_j \frac{\partial (\alpha_k p_k)}{\partial x_j} \end{bmatrix}, \quad \boldsymbol{S} = \begin{bmatrix} S_{\alpha,k} \\ 0 \\ 0 \\ -\overline{p} S_{\alpha,k} \end{bmatrix}. \end{split} \]

where $k=1,2,...,m$ and $m$ is the number of materials. $\alpha_k$ is the volume-fraction of material- $k$ . Bulk properties such as density $\overline{\rho}$ , pressure $\overline{p}$ , and internal energy $\overline{\rho e}$ are defined as,

\[ \begin{split} \overline{\phi} = \sum_k \alpha_k \phi_k, \end{split} \]

where $\phi_k$ is the material density $\rho_k$ , material pressure $p_k$ , or material internal energy $\rho_k e_k$ as required. The specific total energy of material- $k$ is, $E_k = e_k + u_j u_j/2$ , and its specific total enthalpy is, $H_k = E_k + p_k/\rho_k$ , where $e_k$ is the specific internal energy of the material. $Y_k = \alpha_k \rho_k/\overline{\rho}$ is the mass fraction of material- $k$ . The source term $S_{\alpha,k}$ consists of the finite pressure relaxation. Note that although the material total energy equations are written in a non-conservative form, the mixture total energy equation ( $\rho E = \sum_k \alpha_k \rho_k E_k$ ) is a conservative equation (non-conservative terms cancel when summed over all materials). Thus, total energy of the system is conserved. The system is closed by specifying equations of state for each material, usually in the form of $p_k = p_k(\rho_k, e_k)$ .

Finite-rate pressure relaxation

A mixed-cell pressure closure $S_{\alpha,k}$ has to be specified for the above system of equations to be complete. This term accounts for volume fraction redistribution due to differential compaction. This term attempts to model the different amount of compression in materials caused by unequal material compressibilities, by redistributing pressure-induced volume changes based on material bulk modulii. This redistribution of volume fractions results in a finite amount of relaxation between material pressures. The functional form is similar to Tipton's closure. The volume redistribution is given by,

\[ \begin{split} S_{\alpha,k} = \frac{1}{\tau} (p_k-p^*)\frac{\alpha_k}{\mathcal{K}_k}, \end{split} \]

where $p^*$ is the equilibrium pressure that the multi-material cell is expected to reach after sufficient time:

\[ \begin{split} p^* = \frac{\sum_k \left( p_k \frac{\alpha_k}{\mathcal{K}_k} \right)} {\sum_k \frac{\alpha_k}{\mathcal{K}_k}}, \end{split} \]

$\mathcal{K}_k = \rho_k a_{k}^2$ is the material's bulk modulus, and $\tau$ is the pressure-equilibration time-scale,

\[ \begin{split} \tau = c_{\tau} \max_k \left( \frac{h}{a_k} \right). \end{split} \]

The constant $c_{\tau}$ adjusts the pressure equilibration rate of the materials in mixed cells. It controls how fast the pressure relaxation takes place in a mixed cell, as compared to the sound-crossing speed in that cell $(h/a_k)$ . The value of $c_{\tau}$ can be specified in the input file.

Equations of state

Currently, the stiffened-gas equation of state (SG-EoS) is used to close the PDE system. This EoS can be used to describe materials ranging from ideal gases to stiff liquid-like fluids. The internal energy, temperature and speed of sound for material- $k$ using the SG-EoS are:

\[ \begin{split} \rho_k e_k &= \frac{p_k+P_{\text{c}_k}}{\gamma_k-1} + P_{\text{c}_k}, \\ T_k &= \left( \frac{\gamma_k}{\gamma_k-1} \right) \frac{p_k+P_{\text{c}_k}}{\rho_k C_{\text{p}_k}}, \\ a_k &= \sqrt{\gamma_k \frac{p_k+P_{\text{c}_k}}{\rho_k}}, \end{split} \]

where $P_{\text{c}_k}$ , $\gamma_k$ , $C_{\text{p}_k}$ , and $T_k$ are the stiffness parameter, heat capacity ratio, specific heat at constant pressure, and temperature for material-$k$ respectively.

More complex equations of state are a topic of future work.

Temporal discretization

The temporal and spatial derivatives are discretized separately resulting in a method of lines. A 3rd order total variation diminishing (TVD) Runge-Kutta scheme is used to discretize the time derivatives.

Spatial discretization

A modal Reconstructed Discontinuous Galerkin (rDG) method is used to discretize the multi-material system in space. The specific rDG method implemented in Quinoa uses Dubiner's orthogonal basis functions (see the DG page for more details). Reconstruction to higher order is via a least-squares procedure based on node-neighbors. Fluxes through cell-faces are calculated using multi-material Riemann solvers. Non-linear instabilities are suppresed using a vertex-based limiter. Further details can be found in the Inciter papers.

Algebraic interface reconstruction

Tangent of Hyperbola for INterface Capturing (THINC) is used to algebraically reconstruct material interfaces in the domain. THINC assumes a smooth distribution of volume fractions across the interface, using the $\tanh$ function. This results in sharp interface capturing, without expensive geometric procedures. The THINC scheme can be tuned to capture material interfaces in 2-4 cells. Due to the smooth nature of the $\tanh$ function, stability concerns are mitigated.

For further details about the numerical method, see Inciter papers.

Example of multi-material setup

For an example test set up see the example gas impact problem.