An efficient implicit scheme for the multimaterial Euler equations in Lagrangian coordinates

TL;DR

The paper develops a fully implicit finite-volume scheme for 1D multimaterial Euler equations in Lagrangian mass coordinates, enabling large time steps in stiff layered flows by reducing the problem to a single implicit discrete wave equation for the pressure with an SPD structure. It couples a nonconservative smooth predictor with a conservative non-oscillatory corrector, augmented by filtering to suppress pressure and density oscillations across material interfaces, and supports high-order time integration via SDIRK schemes. A mass-constrained, smoothly varying mesh in the mass coordinate is introduced to handle extreme density ratios without introducing spurious diffusion. The method is validated on a broad suite of tests, including single- and multi-material Riemann problems, shock–bubble interactions, and highly stratified systems, showing robust accuracy and efficiency and compatibility with homogenised Kapila-type models in appropriate limits. These contributions enable reliable, interface-preserving simulations of stratified metamaterials and layered multimaterial flows at scales where explicit-Lagrangian methods become prohibitive.

Abstract

Stratified fluids composed of a sequence of alternate layers show interesting macroscopic properties, which may be quite different from those of the individual constituent fluids. On a macroscopic scale, such systems can be considered a sort of fluid metamaterial. In many cases each fluid layer can be described by Euler equations following the stiffened gas equation of state. The computation of detailed numerical solutions of such stratified material poses several challenges, first and foremost the issue of artificial smearing of material parameters across interface boundaries. Lagrangian schemes completely eliminate this issue, but at the cost of rather stringent time step restrictions. In this work we introduce an implicit numerical method for the multimaterial Euler equations in Lagrangian coordinates. The implicit discretization is aimed at bypassing the prohibitive time step restrictions present in flows with stratified media, where one of the materials is particularly dense, or rigid (or both). This is the case for flows of water-air mixtures, air-granular media, or similar high density ratio systems. We will present the novel discretisation approach, which makes extensive use of the remarkable structure of the governing equations in Lagrangian coordinates to find the solution by means of a single implicit discrete wave equation for the pressure field, yielding a symmetric positive definite structure and thus a particularly efficient algorithm. Additionally, we will introduce simple filtering strategies for counteracting the emergence of pressure or density oscillations typically encountered in multimaterial flows, and will present results concerning the robustness, accuracy, and performance of the proposed method, including applications to stratified media with high density and stiffness ratios.

Figures (17)

• Figure 1: Motivation for the filtering techniques introduced in this work. On the left: central fluxes for the density update introduce spurious oscillations in correspondence of contact discontinuities. The oscillations are eliminated by numerical diffusion-like filtering (Section \ref{['sec:densityfiltering']}). On the right: pressure oscillations appear at material interfaces, across which the parameters of the equation of state are discontinuous (the interface positions are highlighted with dashed lines). The oscillations disappear with filtering (Section \ref{['sec:pressurefiltering']}). Both simulations were run with SDIRK2 time-stepping and $N=1000$ cells. For the left panel $k_\mathrm{CFL}=5.0$ and for the right panel $k_\mathrm{CFL}=1000$.
• Figure 2: Illustration of the graded mesh generation procedure. In the left panel, we show the cell mass content $\Delta m$ as a function of the auxiliary (reference) coordinate $z$, along with the mesh spacing that would be induced by fixed uniform value of $\Delta x$, and the mass density function $\mathop{\mathrm{d\!}}\nolimits m/\mathop{\mathrm{d\!}}\nolimits z$ from which $\Delta m$ is integrated (graphically rescaled to match the plot of $\Delta m$). In the central panel, we show how the smooth transient in reference coordinates translates to Eulerian coordinates (that is, as a function of $x$). In the right panel we show how several copies of the half-layer-couples are connected in a sequence of four layer couples. The regions corresponding to the central panel are highlighted. Note the presence of a boundary effect on the mesh spacing due to the fact that the first and last half-layer of the domain has a constant mesh spacing, unlike the internal ones.
• Figure 3: Behaviour of the proposed implicit Lagrangian numerical method at different uniform mesh sizes, from $N=250$ to $N=16000$ cells. The number of points found at contact (0 to 1) and at shocks (2-3) is largely independent of mesh resolution, like it would be for an explicit Lagrangian scheme. Major oscillations are absent, despite the relatively large CFL number ($k_\mathrm{CFL}=10.0$) and the presence of mild shockwaves. The simulations employ a second order diagonally implicit Runge--Kutta time integrator (SDIRK2).
• Figure 4: Behaviour of the proposed implicit Lagrangian numerical method at different CFL values, ranging from $k_\mathrm{CFL} = 1.0$ to $k_\mathrm{CFL} = 50.0$. The first timestep of the simulations at $k_\mathrm{CFL}=25.0$ and $k_\mathrm{CFL}=50.0$ uses $k_\mathrm{CFL}^0=k_\mathrm{CFL}/10$. Some oscillations and additional artificial diffusion appear for the higher CFL numbers but the overall structure of the solution is well preserved. The simulations employ a third order diagonally implicit Runge--Kutta time integrator (SDIRK3) and a constant-mass mesh of $N=1000$ cells. The run at $k_\mathrm{CFL}=50.0$ is carried out in 16 timesteps.
• Figure 5: Numerical results for the Lax shock tube problem and for the modified Sod problem with transsonic rarefaction (Toro 1). The simulations employ a second order diagonally implicit Runge--Kutta time integrator (SDIRK2) and a constant-mass mesh of $N=1000$ cells.