Riemann Solver Research Articles

Numerical simulation of phase transition with the hyperbolic Godunov-Peshkov-Romenski model

In this paper, a thermodynamically consistent numerical solution of the interfacial Riemann problem for the first-order hyperbolic continuum model of Godunov, Peshkov and Romenski (GPR model) is presented. In the presence of phase transition, interfacial physics are governed by molecular interaction on a microscopic scale, beyond the scope of the macroscopic continuum model in the bulk phases. The developed approximate two-phase Riemann solvers tackle this multi-scale problem, by incorporating a local thermodynamic model to predict the interfacial entropy production. Using phenomenological relations of non-equilibrium thermodynamics, interfacial mass and heat fluxes are derived from the entropy production and provide closure at the phase boundary. We employ the proposed Riemann solvers in an efficient sharp interface level-set Ghost-Fluid framework to provide coupling conditions at phase interfaces under phase transition. As a single-phase benchmark, a Rayleigh-Bénard convection is studied to compare the hyperbolic thermal relaxation formulation of the GPR model against the hyperbolic-parabolic Euler-Fourier system. The novel interfacial Riemann solvers are validated against molecular dynamics simulations of evaporating shock tubes with the Lennard-Jones shifted and truncated potential. On a macroscopic scale, evaporating shock tubes are computed for the material n-Dodecane and compared against Euler-Fourier results. Finally, the efficiency and robustness of the scheme is demonstrated with shock-droplet interaction simulations that involve both phase transfer and surface tension, while featuring severe interface deformations.

A staggered Lagrangian magnetohydrodynamics method based on subcell Riemann solver

This paper uses a general formalism to derive staggered Lagrangian method for 2D compressible magnetohydrodynamics (MHD) flows. A subcell method is introduced to discretize the MHD system and some Riemann problems over subcells are solved at the cell center and grid node respectively. In these solvers, only the fast-waves in all jumping relations are considered and thus the solution structure is simple. The discrete conservations of mass, momentum and energy are preserved naturally in the proposed numerical method. In order to meet the thermodynamic Gibbs relation in isentropic flows, an adaptive Riemann solver is implemented at the cell center, in which a criterion is proposed to reduce overheating errors in the rarefying problems and maintains the excellent shock-capturing ability simultaneously. It is worth to be noticed that the divergence-free condition is naturally satisfied in the Lagrangian method. Various numerical tests are presented to demonstrate the accuracy and robustness of the algorithm.

A corrected Raviart–Thomas Spectral Difference scheme stable for arbitrary order of accuracy on triangular and tetrahedral meshes

A theoretical stability estimate is presented for the Raviart–Thomas Spectral Difference scheme for triangular and tetrahedral meshes. It reveals the difficulty of guaranteeing the existence of stable flux points for an arbitrary order of accuracy. To obtain a numerical method stable for any order of accuracy and any space dimension, the present article introduces a new correction term of the flux approximation to cancel the numerical error caused by the Riemann solver exactly. This modification enables the introduction of a new family of so-called Corrected Raviart–Thomas Spectral Difference stable for any polynomial order of approximation. Stability estimates are provided for the L2 norm and a weighted broken Sobolev norm. Numerical experiments in 2D and 3D are conducted for solving linear and non-linear hyperbolic equations with a polynomial order up to p=8 for triangular elements and up to p=5 for tetrahedral elements. They confirm the accuracy, stability, and robustness of the newly introduced numerical scheme. The new method is shown to have a computational cost similar to the Raviart–Thomas Spectral Difference with – for a given mesh – a higher accuracy.

A family of TENOA-THINC-MOOD schemes based on diffuse-interface method for compressible multiphase flows

In this work, we develop a family of hybrid algorithms to simulate compressible multi-phase flows by solving the quasi-conservative five-equation model. The general numerical framework utilizes targeted essentially non-oscillatory schemes with adaptive dissipation (TENOA) to boost the resolution of high-frequency waves, whereas the Tangent of Hyperbola for INterface Capturing (THINC) scheme is activated near the shock and contact discontinuities as well as material interfaces to maintain their sharpness, based on a novel indicator. Also, the newly derived THINC reconstruction scheme is both simpler and better at preserving the flow symmetry. The robustness of our numerical approach is further fortified through the adaptation of a modified a posteriori multidimensional optimal order detection (MOOD) technique. Moreover, to mitigate carbuncle phenomenon encountered in two-dimensional simulations, a hybrid Riemann solver is also proposed based on the low Mach modification around shock for the Harten-Lax-van Leer contact (HLLC) approximate Riemann solver. Numerical results of the challenging one- and two-dimensional benchmark cases demonstrate the superiority of the proposed methods regarding oscillation-free, interface-sharpening, low-dissipation and robustness properties.

Tangential artificial viscosity to alleviate the carbuncle phenomenon, with applications to single-component and multi-material flows

This paper describes a novel approach to alleviate the carbuncle phenomenon which consists in adding to any carbuncle prone Riemann solver an extra viscosity term in tangential momentum flux and its contribution to the energy conservation equation. This term contains one numerical parameter only, a scalar viscosity, which is reduced using a face-based shear detector to preserve shear waves.The idea stems from the investigation of some of the existing Riemann solvers, also presented in the paper. Indeed, when splitting the numerical flux into the face normal and tangential components, we observe that all the carbuncle free Riemann solvers present in the tangential part a numerical viscosity which scales with the sound speed when the normal flow velocity becomes zero. Opposite, in the carbuncle prone solvers this viscosity scales with the normal flow velocity. In particular the carbuncle free HLLCM scheme proposed by Shen et al. can be written by adding to the carbuncle prone HLLC scheme a tangential artificial viscosity term. Then the same can be done for any other Riemann solver, which renders the approach easy to implement in CFD codes for compressible flows.Numerical experiments shows the efficiency of the approach in computing carbuncle free single-component and multi-material flows.

المحاكاة العددية للتدفقات المضغوطة المكو نة من عنصرين مع المعادلة العامة للحالة

Numerical simulations of compressible two-component flows with a general equation of state (EOS), written in the form of the Mie-Gru¨neisen EOS, have been conducted using the diffuse interface method on a structured mesh. The unsteady and inviscid one-dimensional sixequation model of Kapila is employed to describe compressible two-component flows. The model is hyperbolic and non-conservative. The solution of the hyperbolic equations, including the non-conservative equation for the volume fraction evolution, is obtained using an extended Harten-Lax-Leer (HLL) approximate Riemann solver. A general formulation for various EOSs is proposed to enable simulations of a wide range of applications. In which the two constituents are either governed by the same EOS or different types of EOS. In this model, both fluids have the same velocity, but each fluid has its own pressure. Therefore, the pressure relaxation process is performed instantaneously to drive both constituents' pressures towards equilibrium. Several computational test problems were conducted in one and two dimensions to show the ability of the numerical method to simulate such problems. The computed results are presented without introducing spurious pressure oscillations in the method.

A method to represent the strain hardening effect in the hyper-elastic model within a fully Eulerian framework

We propose a method to treat the strain hardening effect in the hyperelastic model. Coupled with the level set interface capturing method and a real ghost fluid method, it is formulated in a fully Eulerian framework. The HLLD approximate Riemann solver is used to obtain the interface state. We use the elastic inverse deformation gradient tensor to track the material deformation state. Once the solid starts to yield, along the normal direction in stress space, an iterative algorithm is used to relax the elastic inverse deformation gradient tensor to target value where the equivalent stress meets the requirement of consistency condition in plasticity. This approach is intuitive and can avoid the use of the physical delay time, which is difficult to calibrate from experiments. We use the Wilkins flying plate impact, cylindrical shell collapse, Taylor impact and the conical-nosed projectile penetration test cases to validate the method. The results show good agreement with the theoretical analysis and experiments. This approach can also apply without further extension when thermal softening effect or damage models are considered.

A well-balanced finite volume solver for the 2D shallow water magnetohydrodynamic equations with topography

In this paper, a second-order finite volume Non-Homogeneous Riemann Solver is used to obtain an approximate solution for the two-dimensional shallow water magnetohydrodynamic (SWMHD) equations considering non-flat bottom topography. We investigate the stability of a perturbed steady state, as well as the stability of energy in these equations after a perturbation of a steady state using a dispersive analysis. To address the elliptic constraint ∇⋅hB=0, the GLM (Generalized Lagrange Multiplier) method designed specifically for finite volume schemes, is used. The proposed solver is implemented on unstructured meshes and verifies the exact conservation property. Several numerical results are presented to validate the high accuracy of our schemes, the well-balanced, and the ability to resolve smooth and discontinuous solutions. The developed finite volume Non-Homogeneous Riemann Solver and the GLM method offer a reliable approach for solving the SWMHD equations, preserving numerical and physical equilibrium, and ensuring stability in the presence of perturbations.

Numerical simulating the blood flow model via nonhomogeneous Riemann solver scheme

This paper examines a one-dimensional (1D) model that appears in arterial blood flow. The mathematical model for blood flow via arteries is similar to that of unstable incompressible flows in thin-walled collapsible tubes. We present the Riemann invariants of the suggested model, which is one of the fundamental components of this work. For numerical modeling of blood flow model, we present a nonhomogeneous Riemann solver (NHRS) technique. Next, we demonstrate the simulation of how pressure, velocity, and cross section area waveforms propagate through arteries. Specifically, we present numerical test cases with various initial data sets. In addition, we compare the NHRS scheme to the classic Rusanov, Lax–Friedrichs, and Roe schemes. Theoretical models for thin-walled collapsible tubes are applicable to a wide range of physiological events and may be used to build clinical devices for actual biomedical science. The NHRS method’s accuracy and efficiency are demonstrated by the numerical tests.

Locally divergence-free well-balanced path-conservative central-upwind schemes for rotating shallow water MHD

We develop a new second-order flux globalization based path-conservative central-upwind (PCCU) scheme for rotating shallow water magnetohydrodynamic equations. The new scheme is designed not only to maintain the divergence-free constraint of the magnetic field at the discrete level but also to satisfy the well-balanced (WB) property by exactly preserving some physically relevant steady states of the underlying system. The locally divergence-free constraint of the magnetic field is enforced by following the method recently introduced in Chertock et al. (2024) [19]: we consider a Godunov-Powell modified version of the studied system, introduce additional equations by spatially differentiating the magnetic field equations, and modify the reconstruction procedures for magnetic field variables. The WB property is ensured by implementing a flux globalization approach within the PCCU scheme, leading to a method capable of preserving both still- and moving-water equilibria exactly. In addition to provably achieving both the WB and divergence-free properties, the new method is implemented on an unstaggered grid and does not require any (approximate) Riemann problem solvers. The performance of the proposed method is demonstrated in several numerical experiments that confirms robustness, a high resolution of obtained results, and a lack of spurious oscillations.

Embedment of WENO-Z reconstruction in Lagrangian WLS scheme implemented on GPU for strongly-compressible multi-phase flows

It is a notoriously challenge problem for Lagrangian particle methods to solve compressible flows with strong discontinuities and large volume variations. To overcome the problems of numerical oscillations, an accurate and stable Lagrangian Weighted-Least-Square (LWLS) method embedded with WENO-Z scheme (LWLS-WENO-Z) is proposed. To further enhance the numerical accuracy and stability, the particle shifting technique, the multi-level local-refinement technique, and the interface repulsive force are appropriately employed in the numerical scheme. Moreover, the smooth and stable coupled particle method between the LWLS and Roe's Riemann solver (LWLS-RR) is considered for the comparisons, further demonstrating the merit of the proposed LWLS-WENO-Z scheme. To reduce the simulation time, parallel computing with a CUDA-based program for the proposed scheme is implemented on a single GPU. In the numerical experiments, several 1D/2D benchmarks are solved to test the accuracy of the proposed method. Particularly, the challenge 2D underwater explosion and the cavitation bubble collapse and jetting are numerically investigated. Compared to other reference solutions, the present coupled particle method demonstrates robustness and superior accuracy in simulating these strongly-compressible multi-phase flows.

A genuinely three-dimensional Riemann solver based on self-similar variables in curvilinear coordinates

A genuinely three-dimensional Riemann solver named MuSIC-3DC (Multidimensional, Self-similar, strongly-Interacting, Consistent, Three-dimensional, Curvilinear coordinates) in curvilinear coordinates is proposed to simulate three-dimensional complex engineering problems. Following Balsara's idea, a three-dimensional Riemann problem at each grid vertex is studied. In order to accurately capture contact discontinuities, the linear variations are assumed in the strongly-interacting zone. Also, the self-similar variables are adopted to simplify the derivation according to the self-similar evolution of the Riemann problem, and the strongly-interacting state could be determined by combining the contributions of all the two-dimensional Riemann problems at the surfaces bounded the strongly-interacting state. Besides, the one-dimensional Riemann problem at the center of each cell interface is solved by the one-dimensional HLLE scheme. After these, the numerical flux at each cell interface is obtained by assembling the three-dimensional fluxes at the vertices and the one-dimensional flux at the center of the cell interface. Several numerical test cases are simulated to validate the proposed solver. The spherical blast wave and three-dimensional Riemann problems indicate that the MuSIC-3DC scheme is capable of capturing three-dimensional self-similar flow structures accurately and suppressing the mesh imprinting phenomenon effectively. The transonic case demonstrates the MuSIC-3DC scheme's high resolution in capturing shock waves. In the hypersonic cases, the MuSIC-3DC scheme is more accurate in predicting stagnation and reattachment heating loads. Moreover, the MuSIC-3DC scheme is capable of obtaining more physical surface heating distribution by considering the information traveling both normal and transverse to the cell interfaces.

A Study of Fifth-Order WENO Reconstruction for Genuinely Two-Dimensional Convection-Pressure Flux Split Riemann Solver

A Study of Fifth-Order WENO Reconstruction for Genuinely Two-Dimensional Convection-Pressure Flux Split Riemann Solver

Cholla-MHD: An Exascale-capable Magnetohydrodynamic Extension to the Cholla Astrophysical Simulation Code

We present an extension of the massively parallel, GPU native, astrophysical hydrodynamics code Cholla to magnetohydrodynamics (MHD). Cholla solves the ideal MHD equations in their Eulerian form on a static Cartesian mesh utilizing the Van Leer + constrained transport integrator, the HLLD Riemann solver, and reconstruction methods at second and third order. Cholla’s MHD module can perform ≈260 million cell updates per GPU-second on an NVIDIA A100 while using the HLLD Riemann solver and second order reconstruction. The inherently parallel nature of GPUs combined with increased memory in new hardware allows Cholla’s MHD module to perform simulations with resolutions ∼5003 cells on a single high-end GPU (e.g., an NVIDIA A100 with 80 GB of memory). We employ GPU direct Message Passing Interface to attain excellent weak scaling on the exascale supercomputer Frontier, while using 74,088 GPUs and simulating a total grid size of over 7.2 trillion cells. A suite of test problems highlights the accuracy of Cholla’s MHD module and demonstrates that zero magnetic divergence in solutions is maintained to round off error. We also present new testing and CI tools using GoogleTest, GitHub Actions, and Jenkins that have made development more robust and accurate and ensure reliability in the future.

Hypersonic similarity for steady compressible full Euler flows over two-dimensional Lipschitz wedges

We establish the optimal convergence rate to the hypersonic similarity law, which is also called the Mach number independence principle, for steady compressible full Euler flows over two-dimensional slender Lipschitz wedges. Mathematically, it can be formulated as the comparison of the entropy solutions in BV∩L1 between the two initial-boundary value problems for the compressible full Euler equations with parameter τ>0 and the hypersonic small-disturbance equations (the scaled compressible full Euler equations with parameter τ=0) with curved characteristic boundaries. We establish the L1–convergence estimate of these two solutions with the optimal convergence rate, which justifies the Van Dyke's similarity theory rigorously for the compressible full Euler flows. This is the first mathematical result on the comparison of two solutions of the compressible Euler equations with characteristic boundary conditions. To achieve this, we first employ the special structures of the two systems and establish the global existence and the L1–stability of the entropy solutions via the wave-front tracking scheme under the smallness assumption on the total variation of both the initial data and the tangential slope function of the wedge boundary. Based on the L1–stability properties of the approximate solutions to the scaled equations with parameter τ>0, a uniform Lipschtiz continuous map with respect to the initial data and the wedge boundary is obtained, which is the first time for the characteristic boundary conditions. Next, we compare the solutions given by the Riemann solvers of the two systems by taking the boundary perturbations into account case by case. Then, for a given fixed hypersonic similarity parameter, as the Mach number tends to infinity, by employing the Lipschitz continuous properties of the map, we establish the desired L1–convergence estimate with the optimal convergence rate. Finally, we show the optimality of the convergence rate by investigating a special solution.

A face-centred finite volume method for laminar and turbulent incompressible flows

This work develops, for the first time, a face-centred finite volume (FCFV) solver for the simulation of laminar and turbulent viscous incompressible flows. The formulation relies on the Reynolds-averaged Navier–Stokes (RANS) equations coupled with the negative Spalart–Allmaras (SA) model and three novel convective stabilisations, inspired by Riemann solvers, are derived and compared numerically. The resulting method achieves first-order convergence of the velocity, the velocity-gradient tensor and the pressure. FCFV accurately predicts engineering quantities of interest, such as drag and lift, on unstructured meshes and, by avoiding gradient reconstruction, the method is less sensitive to mesh quality than other FV methods, even in the presence of highly distorted and stretched cells. A monolithic and a staggered solution strategies for the RANS-SA system are derived and compared numerically. Numerical benchmarks, involving laminar and turbulent, steady and transient cases are used to assess the performance, accuracy and robustness of the proposed FCFV method.

Macroscopic modeling of urban flood inundation through areal-averaged Shallow-Water-Equations

An areal-averaged form of classical Shallow-Water-Equations is developed in conjunction with Finite-Volume-Method for capturing sub-grid bed variation. The averaging mechanism treats sub-grid obstacles through depth-dependent-area-averaged porosity at the macroscopic level. This porosity assumes a binary distribution (0,1) for a resolution fine enough to treat bed-variation separately, resulting in convergence of the developed framework to classical form. An attempt has been made to incorporate the unresolved fine-scale flow-information (e.g., micro-scale and cross-scale interaction components) in terms of the macroscopic variables through a non-linear closure model. An augmented approximated Riemann solver incorporates varying source–sink terms within interfacial fluxes along with discontinuous porosity and bed variation. The model is applied to three test-cases ranging from wave-interaction with trapezoidal porous block to dam-break flows through obstacle(s) with varying grid configurations. The coarse-scale formulation, along with closure, produces a reasonably accurate solution with minimal computational overhead.

Improved local time-stepping schemes for storm surge modeling on unstructured grids

This paper presents improved explicit local time-stepping (LTS) schemes of both first and second order accuracy for storm surge modeling. The two-dimensional shallow water equations are numerically solved on unstructured triangular meshes using finite volume method with Roe’s approximate Riemann solver. The LTS algorithms are designed based on explicit Euler and strong stability preserving Runge–Kutta time integration methods. A single-layer interface prediction–correction scheme is adopted to combine coarse and fine time discretization, further enhancing the stability of the LTS schemes, particularly at higher LTS levels and during long time simulations. An ideal numerical test validates the efficiency of the improved LTS models, revealing their capability to improve computational speed while preserving conservation properties and reducing accuracy loss as LTS levels increase. We further apply the LTS models to cross-scale simulations of storm surges in the Northwest Pacific. Results show that compared to the global time-stepping (GTS) models, the LTS models significantly boost computational speed by up to 37%, all while delivering equally reliable computational outcomes. With expanding high-resolution coastal data and the need for high-resolution modeling, the improved LTS models show great potential for cross-scale storm surge modeling.

Numerical modelling of large elasto-plastic multi-material deformations on Eulerian grids

Abstract The paper deals with an Eulerian model for heterogeneous multiphase (multi-material) elasto-plastic media based on the diffuse interface method in the one-dimensional uniaxial strain approximation. In this approximation, an equilibrium hypoelastic Wilkins bi-material model is derived. This is carried out on the basis of an analogy between the elasto-plastic Wilkins and hydrodynamic Euler models. For the obtained model, a Godunov-type numerical method is developed based on the approximate Riemann solver HLLC. The results of the proposed Eulerian diffuse interface model are compared with reference solutions on moving grids with explicit tracking interfaces. It is shown that the present diffuse interface numerical model describes with good accuracy strong shock-wave processes in heterogeneous multiphase elasto-plastic media.

A Genuinely Two-Dimensional Approximate Riemann Solver with Stress Continuity for Hypo-Elastic Solids

A Genuinely Two-Dimensional Approximate Riemann Solver with Stress Continuity for Hypo-Elastic Solids