Module materials 

Materials: fluid models and interface responses

The crate merlin-materials holds what a material is (its constitutive law) and what one side of an interface answers about it. It depends on nothing; the systems that use it (merlin-systems) and the numerics that consume those (merlin-finite-volume) are elsewhere (crates). Nothing in it knows a grid.

Fluid models

A FluidModel closes the Euler equations with $p=p(\rho,e)$, the pressure as a function of the density and the specific internal energy (Menikoff and Plohr 1989, Section II). From $p$, its inverse $e(\rho,p)$ and the two partial derivatives $p_\rho$ (at fixed $e$) and $p_e$ (at fixed $\rho$) everything else follows:

c^2=p_\rho+\frac{p}{\rho^2}\,p_e,\qquad \Gamma=\frac{p_e}{\rho},\qquad
\left(\frac{\partial(\rho e)}{\partial\rho}\right)_p=e-\frac{\rho\,p_\rho}{p_e}.

The last is the energy that the density jump of the entropy wave carries beside its kinetic energy; it enters the eigenvectors of Euler<D, M> and is exactly zero for the ideal and stiffened gases, which say so. Three models exist:

Model	Law	Used for
Ideal	$p=(\gamma-1)\rho e$	gases
Stiffened	$p=(\gamma-1)\rho e-\gamma p_\infty$ (Ivings, Causon and Toro 1998)	liquids: water $\gamma=4.4$, $p_\infty=6\times10^8$ Pa
MieGruneisen	$p=p_{\text{ref}}(\rho)+\Gamma(\rho)\rho\,(e-e_{\text{ref}}(\rho))$ (Menikoff and Plohr 1989, Section VI)	metals under shock

The Mie-Gruneisen reference curve is the shock Hugoniot of the material at rest built from $U_s=c_0+s\,u_p$ ($p_H=\rho_0c_0^2\eta/(1-s\eta)^2$, $e_H=p_H\eta/(2\rho_0)$, $\eta=1-\rho_0/\rho$), or a constant pressure with zero energy, and $\Gamma$ is constant or $\rho\Gamma$ is. With a constant reference pressure $-\gamma p_\infty$ and $\Gamma=\gamma-1$ it is the stiffened gas, and with $p_\infty=0$ the ideal gas. The choices inside the model are values matched at run time; the model itself is a type parameter, so a call costs nothing and the solver is compiled once for each fluid model and system shape (see limitations).

Checked by the derived quantities against finite differences of the model’s own pressure (the speed of sound as the derivative along an isentrope, $p_\rho$ and $p_e$, the entropy-wave energy), the inverse against the pressure, Mie-Gruneisen reproducing the stiffened and ideal gases to $10^{-14}$, copper at rest having $c_0$ and $\Gamma_0$, and states built from $U_s=c_0+s\,u_p$ satisfying the model. The eigenvectors of Euler<D, M> satisfy $LR=I$ and diagonalise a finite-difference flux Jacobian for each model in 1, 2 and 3D.

Tests that enforce this

Enforced by crates/merlin-materials/tests/fluid/consistency.rs: the_ideal_gas_is_consistent, the_stiffened_gas_is_consistent, the_mie_gruneisen_model_is_consistent, the_derivative_check_can_fail; crates/merlin-materials/tests/fluid/mie_gruneisen.rs: a_constant_reference_curve_reproduces_the_stiffened_and_ideal_gases, copper_at_rest_has_the_bulk_sound_speed_and_the_reference_gruneisen_coefficient, the_model_satisfies_the_hugoniot_of_its_shock_speed_relation; crates/merlin-systems/tests/equations/eigenvectors.rs: the_gas_eigenvectors_invert_each_other_and_diagonalise_the_flux_jacobian, the_stiffened_gas_eigenvectors_invert_each_other_and_diagonalise_the_flux_jacobian, the_mie_gruneisen_eigenvectors_invert_each_other_and_diagonalise_the_flux_jacobian.

The response of one side of an interface

In a Riemann problem normal to an interface, the waves on each side depend only on that side’s own state and on the load the interface carries (Fedkiw, Aslam, Merriman and Osher 1999; Liu, Khoo and Yeo 2003). So each side can be described alone, by one function: how fast does its material move when the interface pushes on it with a normal traction $t$? That function is the side’s response, $v(t)$. The interface is where the two responses agree.

One side, in its own frame. A response is written for one side with the normal pointing out of that side, at the other. $v(t)$ is the velocity of the side’s material along that normal and $t$ the normal traction (negative in compression; for a fluid $t=-p$). The left and the right of an interface are therefore the same code, and two sides that move together have $v_\text{in}+v_\text{out}=0$ (in one solve the inside is the side whose ghost cell is being filled and the outside the other).

An InterfaceResponse provides:

It gives	Meaning	Used for
$v(t)$ and $dv/dt$	the velocity under traction $t$, and its slope; absent for a side with no wave curve	the exact solve
an acoustic impedance	the slope at the side’s own state (for a fluid $\rho c$)	the linear solves, and the first guess of the exact one
the impedance of a wave	the slope $1/W$ of the chord of the wave curve from the state to the traction $t$: for a shock the Rankine-Hugoniot mass flux, $W^2=(p-p_k)/(1/\rho_k-1/\rho)$; for a rarefaction the chord of the isentrope	the HLLC-type solve
the range of $t$	what the material can carry (a gas cannot be pulled, $p\ge0$; a stiffened gas to $p=-p_\infty$)	the bracket, and separation

The state behind a fluid’s waves (density, velocity and pressure once $t$ and $v$ are known, which the ghost states are built from) is a method of the fluid, FluidResponse::state_at_interface, and not of the trait, so that the trait has no associated type and can be used as &dyn InterfaceResponse: only the solution is shared between materials.

A fluid’s response is its wave curve. For a fluid in the state $(\rho_k,u_n,p_k)$ with sound speed $c_k$ (FluidResponse; Toro 2009, Chapter 4):

v(p)=u_n-f(p),\qquad
f(p)=\begin{cases}
\sqrt{(p-p_k)\left(\dfrac{1}{\rho_k}-\dfrac{1}{\rho}\right)} & p>p_k \quad\text{(a shock)}\\[2.2ex]
\displaystyle\int_{\rho_k}^{\rho}\frac{c\,d\rho}{\rho} & p\le p_k \quad\text{(a rarefaction)}
\end{cases}

where $\rho$ is the density behind the wave: on the Hugoniot for a shock, and on the isentrope for a rarefaction,

e(\rho,p)-e_k=\tfrac12\,(p+p_k)\left(\frac{1}{\rho_k}-\frac{1}{\rho}\right)
\qquad\text{and}\qquad
de=\frac{p}{\rho^{2}}\,d\rho .

$f$ increases, $f(p_k)=0$, and the two branches meet with the same slope $1/(\rho_kc_k)$. The ideal and stiffened gases have $f$ in closed form (the stiffened gas in the shifted pressure $p+p_\infty$); any other model gets it from fluid::numerical, and both assume the Hugoniot is convex and the fluid stable ($c^2>0$; Menikoff and Plohr 1989, Section V). A pressure the model does not reach gives NaN.

How the numerical curve is computed

• A shock ends at the density on the Hugoniot relation above, found by a bracketed root find (Illinois regula falsi); then $f=\sqrt{(p-p_k)(1/\rho_k-1/\rho)}$, with the slope from differentiating the relation implicitly. A shock weaker than $10^{-9}$ of $\rho_kc_k^2$ is a sound wave: the root find cannot see its density jump, and the acoustic limit is exact to the same fraction.
• A rarefaction is integrated along the isentrope in $\ln\rho$ by an adaptive Dormand-Prince pair (relative tolerance $10^{-13}$ per step), and the density where the pressure reaches $p$ is found by Newton’s method with $dp/d\ln\rho=\rho c^2$.

What was measured

• The shock branch satisfies the three Rankine-Hugoniot conditions, the energy with the internal energy of the model itself, so the check does not rest on the formulas of the curve. The rarefaction branch keeps the Riemann invariant $u+2c/(\gamma-1)$ and the isentrope. The two branches meet with equal value and slope, and the slope is the derivative.
• The numerical curve reproduces the closed forms of the ideal and stiffened gases to $3\times10^{-12}$ at worst (the function, its slope and the density, for pressure ratios $10^{-3}$ to $10^3$), and so does a Mie-Gruneisen model configured to be a stiffened gas.
• A shock into copper at rest follows $U_s=c_0+s\,u_p$ to $10^{-9}$. The rarefaction of copper released from 48 GPa satisfies the velocity integral and the first law along the isentrope, by a plain quadrature in the pressure that does not use the integration that produced it.

Tests that enforce this

Enforced by crates/merlin-materials/tests/interface/wave_curves.rs: the_shock_branch_satisfies_rankine_hugoniot, the_rarefaction_branch_keeps_the_riemann_invariant_and_the_isentrope, the_two_branches_meet_at_the_state_with_equal_value_and_slope, the_slope_of_the_wave_curve_is_its_derivative, the_rarefaction_to_the_lowest_pressure_has_a_finite_velocity; crates/merlin-materials/tests/interface/numerical_curves.rs: the_numerical_curve_reproduces_the_closed_forms_of_the_ideal_and_stiffened_gases, a_mie_gruneisen_model_that_is_a_stiffened_gas_has_its_curves_on_both_branches, a_shock_into_copper_at_rest_follows_the_linear_shock_speed_relation, the_rarefaction_of_copper_satisfies_the_velocity_and_energy_integrals; crates/merlin-materials/tests/fluid/shock.rs: the_mass_flux_of_every_model_satisfies_the_hugoniot_energy_relation, the_generic_root_find_agrees_with_the_closed_forms, a_vanishing_shock_has_the_acoustic_mass_flux.

The interface law and the matching solve

The law. An InterfaceLaw says which quantities must agree across the interface. Phase 9 has one, Slip:

v_\text{in}+v_\text{out}=0,\qquad t_\text{in}=t_\text{out},\qquad\text{the tangential velocity may jump.}

The two sides move together, they carry the same load, and they may slide past each other. For two inviscid fluids this is not a choice but the contact condition of the Euler equations (Rankine-Hugoniot across a contact: $[u_n]=0$, $[p]=0$, and a vortex sheet for the tangential velocity), so the law follows from the pair of systems and is not a key of the problem file. It stops being right with viscosity, with surface tension, and between solids, where later phases add arms.

The solve. With one traction $t$ for both sides, the interface state is the $t$ at which the two velocities cancel. That is the crossing in the figure above, and solve_interface_problem finds it:

g(t)\;=\;v_\text{in}(t)+v_\text{out}(t)\;=\;0 .

Each $v$ falls as the load grows more compressive, so $g$ increases with $t$ and has at most one root. For the same gas on both sides this is Toro’s pressure equation, $f_L(p)+f_R(p)+u_R-u_L=0$, and the root is his star pressure.

A root exists if and only if $g$ is not negative at the least compressive traction both sides can carry. Otherwise the sides pull apart faster than either can follow, and they separate (a vacuum between gases, a cavity in a liquid).

In one solve the inside is the side whose ghost cell is being filled and the outside the other; the solution is written in the frame of the inside, and seen_from_inside and seen_from_outside hand each side its values in its own frame (the same traction, the velocity with its sign changed), so no other code decides a sign.

Three solvers, named by the caller (InterfaceSolver):

Solver	What it does	Cost
Acoustic	each side replaced by its tangent at its state, $v_k(t)=v_{k0}+(t-t_{k0})/Z_k$, and $g=0$ solved in closed form	one linear solve
Hllc (the default)	the acoustic solve gives a trial traction; each side reports the impedance $W$ of its wave to it (the chord of its wave curve from its state to the trial) and the linear solve is made again: the HLLC star pressure with the pressure-based wave speeds (Toro 2009, Section 10.5.2) for any fluid model; a secant step from the acoustic trial	two linear solves
Exact	Newton’s method on the wave curves, safeguarded by bisection inside a bracket, stopping when the traction changes by less than $10^{-13}$ of itself	1 to 13 iterations

All three check the solution against the range of traction of both sides. The exact solve is the reference and the one the rest of this page measures.

Outcome	When
a root	the normal case: the traction and velocity of the interface
InterfaceError::Separation	$g<0$ at the least compressive traction; carries the speed of recession (for the linear solvers, that of the linear model)
NoConvergence	the exact iteration’s velocity sum jumps over zero
Inadmissible	a response is not a number, or a linear solution is outside the range
NoWaveCurve	the exact solve of a side that offers no wave curve

Nothing panics. For a solid at rest the pressure is zero, so the scale of the search is taken from the estimate.

Where the approximate solvers are inaccurate (measured against the exact one; air against air at equal density, and air against water, an impedance ratio of about 4000):

• Weak jumps. The acoustic solve is a tangent at each state, so the error of its star pressure is second order in the jump (observed 1.97 to 2.00) and that of its velocity third order between equal densities (the quadratic terms of the two sides cancel) and second order against water. Hllc is one secant step from that trial, which squares the error it starts from: the pressure is fourth order between gases (3.95) and third against water (3.0), and the velocity third order in both. At a 1 % pressure jump in air its pressure error is 2e-6 of the acoustic one; its velocity between equal gases is about twice the acoustic error, at the same order, the one place where Acoustic is better.
• Strong jumps. For a pressure ratio of 1000 the error of the star pressure is 0.93 (Acoustic) and 0.31 (Hllc) between gases, 0.030 and 0.0028 against water: Hllc is better at every ratio tried, from 0.001 to 1000, and neither is good for a huge jump between gases, where the exact solve takes 7 iterations.
• Strong expansions. The acoustic line of a rarefaction reaches zero pressure at a recession of $2c/\gamma$ between equal gases (1.43 c for air), and the true isentrope only at $4c/(\gamma-1)$ (10 c): for a recession between them (the 123 problem recedes at 3.4 c) the exact solve finds the solution and Acoustic reports a separation that is not there. The chord of the isentrope runs to the true value at zero pressure, so Hllc has the true limit: it solves at 9.9 c and reports the separation at 10.1 c, as the exact solve does. At 3.4 c its pressure is valid but crude, about 12 times the exact one (0.66 against 0.055 in the units of the test), the price of one correction.

What was measured

• The same gas on both sides. The solve gives the star pressure and velocity of Toro’s five tests to $10^{-12}$ of the exact solver of this code and to the digits of his tables, and the densities behind the two waves to $10^{-11}$. It takes 6 to 8 iterations, and 13 for the near-vacuum “123” problem, where the slope of the wave curve is steep and the iteration falls back to bisection.
• Different fluids. That solver cannot do these; the solve agrees to $10^{-10}$ with a plain bisection on Toro’s functions (air against a heavier gas, air against a helium-like gas, gas against water). A stiffened gas on both sides agrees to $10^{-11}$ with the exact solution in the shifted pressure.
• One gas, random states. 6000 random pairs of one gas (three gases, every pattern of waves, jumps up to a million in pressure) agree with the exact solver to $10^{-12}$, except for three pairs a hair inside the vacuum limit (a star pressure below $10^{-18}$ of the pressures), where the wave curves are vertical and the solve says NoConvergence or agrees more loosely.
• Symmetry and separation. Swapping the sides gives the same traction and the opposite velocity, for every solver. Two gases receding faster than the rarefaction to a vacuum separate, and two liquids cavitate at $p=-p_\infty$, with the gap the closed form gives.
• Solids. Two copper plates meeting, with or without a stationary target, match $p^*=\rho_0(c_0+s\,u_p)u_p$ found from the shock-speed relation alone.

Tests that enforce this

Enforced by riemann/interface_matching.rs: the_star_values_of_toros_tests_match_the_exact_solver_and_the_published_values, the_densities_behind_the_waves_match_the_exact_solver, two_different_fluids_match_an_independent_bisection, a_stiffened_gas_matches_the_exact_solution_in_the_shifted_pressure; crates/merlin-materials/tests/interface/matching.rs: swapping_the_sides_gives_the_same_traction_and_the_opposite_velocity, two_linear_responses_are_solved_in_closed_form, gases_receding_faster_than_the_rarefaction_to_a_vacuum_separate, liquids_receding_fast_cavitate_at_minus_p_inf, a_velocity_sum_that_jumps_is_a_failure_to_converge_and_a_nan_is_inadmissible; crates/merlin-materials/tests/interface/numerical_curves.rs: copper_plates_meeting_match_the_shock_speed_relation; crates/merlin-materials/tests/interface/solvers.rs: the_wave_impedance_is_the_chord_of_the_wave_curve_for_a_shock_and_for_a_rarefaction, every_solver_is_exact_for_linear_sides, a_contact_is_left_alone_by_every_solver, swapping_inside_and_outside_gives_the_same_traction_and_the_opposite_velocity, the_approximate_solvers_converge_at_the_orders_the_chord_step_predicts, hllc_corrects_the_pressure_by_orders_of_magnitude_and_the_velocity_between_gases_by_nothing, hllc_beats_the_acoustic_solve_where_the_jump_is_strong_and_neither_is_good_for_a_huge_one, every_solver_reports_separation_for_receding_gases, hllc_has_the_true_separation_limit_in_an_expansion_and_the_acoustic_solve_does_not, the_exact_solver_refuses_a_side_without_a_wave_curve; riemann/interface_matching.rs: random_pairs_of_one_gas_match_the_exact_solver_in_every_wave_pattern.

## What this does not cover

• The responses are not connected to a solver yet. Nothing in the merlin program calls the matching solve: it is library code that the coupled solver and the ghost states of phase 9.2 and 9.3 will use, and its validation here is against exact and independent solutions of the Riemann problem.
• The exact Riemann flux is still the ideal-gas solver. Rebuilding it on the matching solve, for any model, is deferred (item 9.1b of the roadmap); the matching solve already reproduces its star values.
• Only the normal component. The vector response (tangential velocity and traction, shear) comes with the solids of phase 10.
• No lower pressure for a solid. A Mie-Gruneisen model has no tension limit, so two plates that recede from each other end in Inadmissible, not Separation; a solid’s strength is phase 10.

References

• Menikoff, R. and Plohr, B. J. (1989), The Riemann problem for fluid flow of real materials, Rev. Mod. Phys. 61: the equation of state (Section II), the Hugoniot and isentrope (Section III), convexity (Section V), the Mie-Gruneisen law (Section VI).
• Toro, E. F. (2009), Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed., Springer, Chapter 4: the wave functions $f_K$ and the pressure equation.
• Ivings, M. J., Causon, D. M. and Toro, E. F. (1998), On Riemann solvers for compressible liquids, Int. J. Numer. Meth. Fluids 28: the stiffened gas.
• 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), J. Comput. Phys. 152; Liu, T. G., Khoo, B. C. and Yeo, K. S. (2003), Ghost fluid method for strong shock impacting on material interface, J. Comput. Phys. 190.
• Dormand, J. R. and Prince, P. J. (1980), A family of embedded Runge-Kutta formulae, J. Comput. Appl. Math. 6.
• Marsh, S. P. (ed.) (1980), LASL Shock Hugoniot Data: the shock-speed relation of copper.

References were cited from memory and are to be checked against the sources before they are relied on.