Harten-Lax-van Leer-Contact Research Articles

Modelling of wall-bounded cavitating flow and spray mixing in multi-component environments using the PC-SAFT equation of state

This work introduces a numerical multiphase model for multi-component mixtures, utilizing tabulated data for physical and transport properties across a spectrum of conditions from near-vacuum pressures to supercritical states. The property data are derived using Perturbed Chain Statistical Associating Fluid Theory (PC-SAFT), vapor-liquid equilibrium (VLE) calculations, entropy scaling methodologies, and Group Contribution (GC) methods. These techniques accurately reflect the thermodynamic behaviors of real fluids, avoiding the empirical estimation of Equation of State (EoS) input parameters. Implemented in OpenFOAM, the fluid dynamics solver is designed to address the three-dimensional Navier-Stokes equations for multi-component mixtures. The methodology integrates operator splitting to manage hyperbolic and parabolic steps distinctively. Hyperbolic terms are solved using the HLLC (Harten-Lax-van Leer-Contact) solver with temporal integration performed via a third-order Strong-Stability-Preserving Runge–Kutta (SSP-RK3) method. Viscous stress tensor contributions in the momentum equation are handled through an implicit velocity correction equation, while parabolic terms in the energy equation are explicitly solved. The simulation efficiency is further enhanced by adaptive Local Time Stepping and the Immersed Boundary (IB) method, which addresses interactions between the fluid and solid boundaries. Turbulence is resolved using the Wall Adaptive Large Eddy (WALE) model. Applied to high-pressure diesel fuel spray injections into non-reacting (nitrogen) gas environments, the model has been validated against Engine Combustion Network (ECN) data for the Spray-C configuration, featuring a fully cavitating multi-hole orifice. Results demonstrate that the model achieves accurate predictions across a broad range of tested conditions without the need for tuning or calibration parameters.

High-order accurate positivity-preserving and well-balanced discontinuous Galerkin schemes for ten-moment Gaussian closure equations with source terms

This paper proposes novel high-order accurate discontinuous Galerkin (DG) schemes for the one- and two-dimensional ten-moment Gaussian closure equations with source terms defined by a known potential function. Our DG schemes exhibit the desirable capability of being well-balanced (WB) for a known hydrostatic equilibrium state while simultaneously preserving positive density and positive-definite anisotropic pressure tensor. The well-balancedness is built on carefully modifying the solution states in the Harten–Lax–van Leer–contact (HLLC) flux, and appropriate reformulation and discretization of the source terms. Our novel modification technique overcomes the difficulties posed by the anisotropic effects, maintains the high-order accuracy, and ensures that the modified solution state remains within the physically admissible state set. We provide the rigorous positivity-preserving analyses of our WB DG schemes, based on several key properties of the admissible state set, the HLLC flux and the HLLC solver, as well as the geometric quasilinearization (GQL) approach in Wu and Shu (2023) [52], which was originally applied to analyze the admissible state set and the physical-constraints-preserving schemes for the relativistic magnetohydrodynamic equations in Wu and Tang (2017) [54], to address the difficulties arising from the nonlinear constraints on the pressure tensor. Moreover, the proposed WB DG schemes satisfy the weak positivity for the cell averages, implying the use of a simple scaling limiter to enforce the physical admissibility of the DG solution polynomials at certain points of interest. Extensive numerical experiments are conducted to validate the preservation of equilibrium states, accuracy in capturing small perturbations to such states, robustness in solving problems involving low density or low pressure, and high resolution for both smooth and discontinuous solutions.

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.

A systematic analysis of three-dimensional Riemann problems for verification of compressible-flow solvers

Numerical simulation is a well-established way of analyzing compressible flows. Due to high computational demands of solvers for such flow problems, their verification is typically limited to two-dimensional (2D) cases. However, 2D simulations suppress fundamental three-dimensional (3D) aspects of the flow evolution in (compressible) turbulent flows or shock-bubble and shock-drop interactions. With increase of computational power, 3D simulations become more feasible for routine analyses. The verification of 3D simulation frameworks is often limited to transformations of lower dimensional test cases in the 3D space. There is a lack of strictly 3D reference test cases for gas dynamics. In this work, we present a set of genuine 3D Riemann problems in order to validate and verify numerical solvers for compressible flows. The problems are designed such that each octant’s constant initial data connects two neighboring states by an elementary wave only. The problem design is inspired by well-established 2D Riemann problems most prominently posted by Lax and Liu (1998). In contrast to the twenty published 2D cases, more than 300 distinct combinations can be found in 3D. We provide example solutions for the particularly interesting ones of these case combinations and show how the cases help to expose shortcomings of numerical solvers. We provide reference data from computations with an open-source compressible multiresolution flow solver. For the reference solutions, we employ the Harten-Lax-van Leer contact (HLLC) Riemann solver and a weighted essentially non-oscillatory (WENO) reconstruction stencil of fifth order. The reference solutions use an effective resolution of one billion cells. We additionally make the full compute pipeline of this work publicly available, so interested researchers may reproduce and extend the current work.

A shock-stable rotated-hybrid Riemann solver on rectangular and triangular grids

The rotated Riemann solver is robust against the carbuncle phenomenon, especially for multidimensional computation. Moreover, hybrid techniques are usually used to enhance the stability of an accurate scheme by combining an accurate scheme with a diffusive scheme. This paper proposes a rotated-hybrid Riemann solver named the rotated-HLLC+ scheme. The scheme is developed by hybridizing the Harten–Lax–van Leer contact (HLLC) scheme with the advection upstream splitting method based on a flux vector splitting (AUSMV+) scheme by following the rotated Riemann solver approach. The unit vector n1 is calculated from the velocity-difference vector, and the unit vector n2 is the orthogonal vector. The linearized analysis suggests that the HLLC scheme should be used in the direction of n1 and the AUSMV+ scheme in the direction n2. In this way, the hybrid scheme becomes shock-stable with less numerical dissipation. Moreover, the pressure-based method is used to detect the shock wave. Several numerical experiments suggest that the pressure cutoff parameter εp=0.01 may be generally suitable and provide a stable solution with little additional numerical dissipation. The last two numerical examples show that the computational performance of the rotated-HLLC+ scheme is comparable to the HLLC scheme for the weak shock reflection over convex double wedges. However, the scheme is approximately 9% faster than the HLLC scheme for the double Mach reflection of a strong shock wave. The proposed scheme gives fast, stable, and accurate solutions on rectangular and triangular grids.

Numerical simulation of shock wave propagation over a dense particle layer using the Baer–Nunziato model

The present study examines the possibility of numerical simulation of a strong shock wave propagating over the surface of a dense layer of particles poured onto an impermeable wall using the Baer–Nunziato two-phase flow model. The setting of the problem follows the full-scale experiment. The mathematical model is based on a two-dimensional system of Baer–Nunziato equations and takes into account intergranular stresses arising in the solid phase of particles. The computational algorithm is based on the Harten–Lax–van Leer–Contact method with a pressure relaxation procedure. The developed algorithm proved to be workable for two-phase problems with explicit interfacial boundaries and strong shock waves. These issues are typical of problems arising from the interaction of a shock wave with a bed or a layer of particles. A comparison with the simulations and full-scale experiments of other authors is carried out. A reasonable agreement with the experiment is obtained for the angles of the transmitted compaction wave and granular contact, including their dependency on the intensity of the propagating shock wave. The granular contact angle increases with the incident shock wave Mach number, while the transmitted compaction wave angle decreases. An explanation is given of the phenomenon of the decrease in thickness of the compacted region in the layer with the increase in intensity of the propagating shock wave. The main reason is that the maximal value of the particle volume fraction in the plug of compacted particles in the layer rises with the increase in shock wave intensity.

Numerical Simulation of Sediment Transport in Unsteady Open Channel Flow

This paper presented a two-dimensional, well-balanced hydrodynamic and sediment transport model based on the solutions of variable-density shallow-water equations (VDSWEs) and the Exner equation for bed change for simulating sediment transport in unsteady flows. Those equations are solved in a coupled way by the first-order Godunov-type finite volume method. The Harten–Lax–van Leer–Contact (HLLC) Riemann solver is extended to find the local Riemann fluxes to maintain the exact balance between the momentum term and the bed slope term. A well-balanced scheme is superior to an unbalanced scheme to minimize numerical dispersion as demonstrated by the synthetic standing contact-discontinuity test case. Following this, the model is employed to simulate two laboratory experiments and a field case, the 1996 Lake Ha! Ha! flood event in Canada. The results of water surface elevations and bed surface profiles agree well with the measurements. The accuracy and robustness of the numerical schemes make the model a good candidate for practical engineering applications.

A cell-centered spatiotemporal coupled method for the compressible Euler equations

A cell-centered spatiotemporal coupled method is developed to solve the compressible Euler equations. The spatial discretization is performed using an improved weighted essentially non-oscillation scheme, where the Harten–Lax–van Leer–contact approximate Riemann solver is used for computing the numerical fluxes. A two-stage fourth-order scheme is adopted to carry out time advancement for unsteady problems. The proposed method is featured by spatiotemporal coupling time-stepping that can be generalized without using the case-dependent generalized Riemann problem solver. A number of one- and two-dimensional test cases are presented to demonstrate the performance of the proposed method for solving the compressible Euler equations on structured grids. The numerical results indicate that the novel method can achieve relatively large Courant–Friedrichs–Lewy (CFL) number compared to other studies that implement the two-stage fourth-order scheme, and that it is more capable of capturing small-scale flow structures than the Runge–Kutta (RK) method.

A unified consistent source term computational algorithm for the [formula omitted]-based compressible multi-fluid flow model

Using the criterion of one-fluid preservation, we developed a consistent algorithm for the source term to solve the γ-based compressible multi-fluid flow model with three approximate Riemann solvers, namely the Lax–Friedrichs (LxF), Kurganov, and Harten–Lax–van Leer contact (HLLC) solvers. The consistent algorithm comprises a standard Godunov solver with a high-order reconstruction and a consistent source term integration part. We prove that the present algorithm is consistent with Abgrall's criterion of the moving-material-interface property in the finite volume method framework. The cell boundary velocity for the source term discretization is found to be the same as HLLC solvers from five-equation model, but for the first time for LxF and Kurganov Riemann solvers. We develop 12 compressible multi-fluid solvers by combining four reconstruction schemes and three approximate Riemann solvers together with their consistent-source-term integration algorithms. All 12 solvers can maintain the one-fluid preservation and moving-material-interface properties numerically in several one- and two-dimensional example cases. The simulation of an underwater explosion demonstrates that a boundary-variation-diminishing reconstruction predicts an interface with width controlled within approximately three cells which is independent of the Riemann solver. The simulation of an explosion near a free surface demonstrates that the proposed model can simulate a severe compressible multi-fluid flow involving an interface across which there are large differences in density, pressure and parameters in the equation of state. In conclusion, the proposed consistent algorithm provides a unified framework for one kind of non-conservative hyperbolic system into a conservative hyperbolic system and a source term with velocity divergence, where the former can be computed by classical Godunov-type algorithms and the latter can be solved by the proposed consistent algorithm.

An interface sharpening technique for the simulation of underwater explosions

In this study, a simple and effective interface sharpening technique is proposed to simulate underwater explosions. The five-equation model is used as the numerical model for underwater explosions. The multi-dimensional problem can be split into multiple one-dimensional problems using the Strang splitting. The MUSCL–Hancock Method (MHM) is employed to solve the governing equations in one-dimensional space, and the HLLC (Harten–Lax–van Leer contact) scheme is adopted to solve the Riemann problem approximately. The numerical results show that the five-equation model cannot accurately simulate underwater explosions due to the diffusion problem. To overcome this problem, a new interface sharpening model based on the five-equation model is obtained by adding the anti-diffusion source term to the volume-fraction transport equation. A fractional step method is employed to solve equations of the interface sharpening model in two separate steps. The model is validated by conducting numerical cases. On this basis, the one-dimensional and two-dimensional underwater explosion simulations are conducted. The numerical results indicate that the interface sharpening model based on the five-equation model can accurately simulate underwater explosions.

A new two-phase shallow water hydro-sediment-morphodynamic model based on the HLLC solver and the hybrid LTS/GMaTS approach

In this paper, we present a new computationally-efficient and high-resolution depth-averaged two-phase flow model for hydro-sediment-morphydynamic processes, featuring an advance over existing models in terms of accuracy and efficiency of numerical solution. Under the framework of finite volume method (FVM) on unstructured grids, the Harten-Lax-van Leer-Contact (HLLC) approximate Riemann solver is proposed to compute inter-cell fluxes by applying the classical upwind HLLC approach to the water-sediment mixture and the sediment phases separately, in contrast to previous two-phase flow models using centered schemes. Moreover, to improve computational efficiency, the local-time-stepping (LTS) approach is implemented, the first attempt in the field of two-phase flow modelling. After a convergence rate study, the model is tested against a series of flow-sediment-bed evolutions induced by (1) two refilling processes of dredged trenches, (2) two instantaneous dam-break flooding flows, and (3) one levee breaching process by overtopping flows. It features encouraging performance when compared to a two-phase flow model based on a centered scheme and global time bound, characterized by more accurate results and much less computational cost. The present modelling framework shows promise in practical shallow water hydro-sediment-morphodynamic modelling applications.

Development of a new OpenFOAM solver for plasma cutting modeling

A new OpenFOAM solver is presented, for the simulation of plasma cutting torches. The mathematical model that is introduced is based on the compressible Navier–Stokes equations coupled via source terms to the electric current conservation equation. Due to the conservative and hyperbolic nature of the model, a Godunov-type scheme is used for the first time in the context of plasma cutting simulation. The numerical method consists of a second-order Total Variation Diminishing (TVD) integration with flux Harten–Lax–van Leer-Contact (HLLC) Riemann solver for the flow conservation equations, coupled with a Laplace solver for the current conservation equation. An efficient formulation for the equation of state, accurately taking into account the plasma properties, is also presented. The solver is validated through a set of canonical test cases (shock tubes and 2D Riemann problems) and it is used to simulate a three-dimensional plasma cutting torch. Good agreement is found with the literature, with an improvement in the ability to deal with the shocks occurring during plasma cutting.

A novel combination of mathematical modelling approaches for simulating nearly hypersonic, viscous, reacting magnetogasdynamic flows

This work focuses on the mathematical combination of magnetogasdynamics (MGD) and the Dunn and Kang chemical kinetic model with a multi-component Harten–Lax–van Leer Contact (HLLC) local Riemann solver using high-order interpolation schemes. Due to the nature of nearly hypersonic, unsteady, viscous, reacting MGD flows, the gas dissociates which is taken into account with 11 species and 26 chemical reactions along with the temperature dependent transport properties through the Wilke model. In the present work, the electromagnetic effects are considered by solving a Poisson equation simultaneously with the chemical reaction equations. The non-linear convective terms are treated numerically with third-, fifth- and seventh-order Weighted Essentially Non-Oscillatory (WENO) interpolation schemes along with second-, fourth- and sixth-order central finite difference schemes for the viscous terms. The combination of these numerical and physical models are analysed for stability and validated for benchmark problems such as (a) Hartmann flow as MGD validation test case and (b) the shock-tube problem for testing the chemical reaction model in which case the reaction rate controlling temperature suggested by Park has also been taken into account. Furthermore, a complex physical problem of an external aerodynamic flow over the cube for the Mach number varying from 1.25 to 3 have also been investigated, because a limited number of studies is available for understanding the flow behaviour of sharp cornered blunt bodies such as cubes. The knowledge of shock structures of these objects has become imperative to predict the atmospheric re-entry trajectories of de-orbiting CubeSats, space debris or fragments of launch vehicles. The main physical finding is that the shock structures are observed to be much more resolved compared to previous works, especially in the wake region. The knowledge gap is bridged by detailing the shock structures and showing the dependence of the shock standoff distance on total enthalpy of the flow.

Fully compressible multiphase model for computation of compressible fluid flows with large density ratio and the presence of shock waves

We propose a fully compressible multiphase model for simulating compressible interfacial flows. The mathematical model is formulated on the basis of the typical conservation laws for mixtures, including mass, momentum, and energy balance. The mass balance equation in each phase is considered to discretely conserve the mass of each phase; therefore, it is always mass conservative for each phase. The numerical solution procedure for this model is developed using a high-resolution shock-capturing finite-volume method based on the Harten–Lax–van Leer-contact (HLLC) Riemann solver and total variation diminishing (TVD) time-integration schemes for simulating unsteady multiphase flows and capturing strong shocks. The finite-volume Riemann solver is combined with a monotonic upstream-centered scheme for conservation laws (MUSCL) and compressive limiters to maintain a sharp interface. The numerical framework yields a conceptually simple but robust and versatile numerical procedure. Numerical results of benchmark Riemann problems demonstrate the ability of the proposed method to obtain sharp interfaces, free of spurious oscillations, high accuracy, conservation properties, and thermodynamic consistency over various Mach numbers.

Oscillation-free nodal discontinuous spectral element method for the simulation of compressible multicomponent flows

The nodal discontinuous spectral element method (DSEM) is improved to avoid pressure and velocity oscillations from material interfaces in the solution of multi-component compressible flows. A novel double flux method compatible with the DSEM scheme is developed to account for variations in the specific heats ratio at material interfaces. The double flux method consists of a modified Harten-Lax-van Leer-Contact (HLLC) Riemann solver and an averaging technique to freeze the temperature-dependent specific heats ratio during the solution update to remove oscillations. The double flux model is combined with a conservative method of calculating the convective flux function to avoid unnecessary energy conservation error at shock fronts. Numerical diffusion terms are also introduced for capturing a shock wave, contact discontinuity, and material interface, without inducing oscillations at the material interface. A series of one- and two-dimensional tests are conducted to evaluate the efficacy of the numerical scheme in capturing sharp discontinuities, as well as vortices, without introducing significant conservation error. According to the benchmark results, the new approach maintains the pressure, velocity, and temperature equilibriums across material interfaces. Artificial diffusion is added to the equation in a controlled manner that enables a higher Courant-Friedrichs-Lewy (CFL) number for high polynomial orders. The present scheme's performance is validated by comparing the spatiotemporal evolution of characteristic points in the shock-bubble interaction and Richtmyer-Meshokov's instability problems.

A graphics processing unit-based robust numerical model for solute transport driven by torrential flow condition

Solute transport simulations are important in water pollution events. This paper introduces a finite volume Godunov-type model for solving a 4×4 matrix form of the hyperbolic conservation laws consisting of 2D shallow water equations and transport equations. The model adopts the Harten-Lax-van Leer-contact (HLLC)-approximate Riemann solution to calculate the cell interface fluxes. It can deal well with the changes in the dry and wet interfaces in an actual complex terrain, and it has a strong shock-wave capturing ability. Using monotonic upstream-centred scheme for conservation laws (MUSCL) linear reconstruction with finite slope and the Runge-Kutta time integration method can achieve second-order accuracy. At the same time, the introduction of graphics processing unit (GPU)-accelerated computing technology greatly increases the computing speed. The model is validated against multiple benchmarks, and the results are in good agreement with analytical solutions and other published numerical predictions. The third test case uses the GPU and central processing unit (CPU) calculation models which take 3.865 s and 13.865 s, respectively, indicating that the GPU calculation model can increase the calculation speed by 3.6 times. In the fourth test case, comparing the numerical model calculated by GPU with the traditional numerical model calculated by CPU, the calculation efficiencies of the numerical model calculated by GPU under different resolution grids are 9.8–44.6 times higher than those by CPU. Therefore, it has better potential than previous models for large-scale simulation of solute transport in water pollution incidents. It can provide a reliable theoretical basis and strong data support in the rapid assessment and early warning of water pollution accidents.

A carbuncle cure for the Harten-Lax-van Leer contact (HLLC) scheme using a novel velocity-based sensor

A hybrid numerical flux scheme is proposed by adapting the carbuncle-free modified Harten-Lax-van Leer contact (HLLCM) scheme to smoothly revert to the Harten-Lax-van Leer contact (HLLC) scheme in regions of shear. This hybrid scheme, referred to as the HLLCT scheme, employs a novel, velocity-based shear sensor. In contrast to the non-local pressure-based shock sensors often used in carbuncle cures, the proposed shear sensor can be computed in a localized manner meaning that the HLLCT scheme can be easily introduced into existing codes without having to implement additional data structures. Through numerical experiments, it is shown that the HLLCT scheme is able to resolve shear layers accurately without succumbing to the shock instability.

Porous Shallow Water Modeling for Urban Floods in the Zhoushan City, China

Typhoon-induced intense rainfall and urban flooding have endangered the city of Zhoushan every year, urging efficient and accurate flooding prediction. Here, two models (the classical shallow water model that approximates complex buildings by locally refined meshes, and the porous shallow water model that adopts the concept of porosity) are developed and compared for the city of Zhoushan. Specifically, in the porous shallow water model, the building effects on flow storage and conveyance are modeled by the volumetric and edge porosities for each grid, and those on flow resistance are considered by adding extra drag in the flow momentum. Both models are developed under the framework of finite volume method using unstructured triangular grids, along with the Harten-Lax-van Leer-Contact (HLLC) approximate Riemann solver for flux computation and a flexible dry-wet treatment that guarantee model accuracy in dealing with complex flow regimes and topography. The pluvial flooding is simulated during the Super Typhoon Lekima in a 46 km2 mountain-bounded urban area, where efficient and accurate flooding prediction is challenged by local complex building geometry and mountainous topography. It is shown that the computed water depth and flow velocity of the two models agree with each other quite well. For a 2.8-day prediction, the computational cost is 120 min for the porous model using 12 cores of the Intel(R) Xeon(R) Platinum 8173M CPU @ 2.00 GHz processor, whereas it is as high as 17,154 min for the classical shallow water model. It indicates a speed-up of 143 times and sufficient pre-warning time by using the porous shallow water model, without appreciable loss in the quantitative accuracy.

A positivity-preserving high-order weighted compact nonlinear scheme for compressible gas-liquid flows

We present a robust, highly accurate, and efficient positivity- and boundedness-preserving diffuse interface method for the simulations of compressible gas-liquid two-phase flows with the five-equation model by Allaire et al. using high-order finite difference weighted compact nonlinear scheme (WCNS) in the explicit form. The equation of states of gas and liquid are given by the ideal gas and stiffened gas laws respectively. Under a mild assumption on the relative magnitude between the ratios of specific heats of the gas and liquid, we can construct limiting procedures for the fifth order incremental-stencil WCNS (WCNS-IS) with the first order Harten–Lax–van Leer contact (HLLC) flux such that positive partial densities and squared speed of sound can be ensured in the solutions, together with bounded volume fractions and mass fractions. The limiting procedures are discretely conservative for all conservative equations in the five-equation model and can also be easily applied to any other conservative finite difference or finite volume scheme. Numerical tests with liquid water and air are reported to demonstrate the robustness and high accuracy of the WCNS-IS with the positivity- and boundedness-preserving limiters even under extreme conditions.

HLLC-type methods for compressible two-phase flow in ducts with discontinuous area changes

In this work, the Harten-Lax-van Leer Contact (HLLC) approximate Riemann solver is extended to two-phase flow through ducts with discontinuous cross-sections. Two main strategies are explored regarding the treatment of the non-conservative term arising in the governing equations. In the first, labelled HLLC+S, the non-conservative term is discretized separately. In the second, labelled HLLCS, the non-conservative term is incorporated in the Riemann solver. The methods are assessed by numerical tests for single and two-phase flow of CO2, the latter employing a homogeneous equilibrium model where the thermodynamic properties are calculated using the Peng–Robinson equation of state. The methods have different strengths, but in general, HLLCS is found to work best. In particular, it is demonstrated to be equally accurate and more robust than existing methods for non-resonant flow. It is also well-balanced for subsonic flow in the sense that it conserves steady-state flow.