Harten Lax Research Articles

Uniformly high-order bound-preserving OEDG schemes for two-phase flows

This paper proposes and rigorously analyzes novel high-order bound-preserving oscillation-eliminating discontinuous Galerkin (BP-OEDG) schemes with the Harten–Lax–van Leer (HLL) numerical flux for the Kapila five-equation two-phase flow model. The evolution equation for the volume fraction is formulated as a conservative advection equation with an additional non-conservative source term. Utilizing a technical splitting approach, we design a uniformly high-order discretization for this source term, incorporating an “upwind” discretization of the non-conservative product at cell interfaces based on the HLL wave speeds. This method inherently ensures the symmetry of the volume fractions, which is crucial for establishing the bound-preserving property. The positivity of partial densities and internal energy is rigorously proven using the geometric quasilinearization (GQL) approach, which transforms the nonlinear pressure positivity constraint into equivalent linear constraints. To suppress potential spurious oscillations, we incorporate a scale-invariant and linearity-invariant oscillation elimination (OE) procedure that damps the DG modal coefficients after each Runge–Kutta stage, as proposed in [M. Peng, Z. Sun and K. Wu, OEDG: Oscillation-eliminating discontinuous Galerkin method for hyperbolic conservation laws, Math. Comput. (2024), https://doi.org/10.1090/mcom/3998]. This OE procedure, acting as a post-processing filter based on the jumps of the DG solution at cell interfaces, is easy to implement, maintains the Abgrall equilibrium condition around an isolated material interface, and preserves the high-order accuracy of the DG schemes. The effectiveness and robustness of the proposed high-order BP-OEDG schemes are demonstrated through several benchmark numerical experiments.

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 finite volume model for maintaining stationarity and reducing spurious oscillations in simulations of sewer system filling and emptying

This paper introduces a model for simulating the unsteady dynamics of sewer systems filling and emptying, offering greater accuracy and stability. This article presents two novel contributions: first, the HLLS scheme (Harten–Lax–van Leer + Source term) is adapted to ensure the preservation of stationary conditions not only in free surface flows, as originally conceived, but also in pressurized flows and mixed flow scenarios. Second, a new method is proposed for the treatment of open channel flow cells near pressurization or adjacent to pressurized cells to minimize spurious oscillations when utilizing the two-component pressure approach (TPA) model. To verify the new model's effectiveness, it was tested for various conditions against the outcomes of the Open Source Field Operation and Manipulation (OpenFOAM) computational fluid dynamics (CFD) model. Furthermore, to demonstrate the model's potential for simulating real systems, the model was applied to three sewer systems that closely resemble real-world conditions, each of which had been intentionally modified for confidentiality purposes. The results show that the improved model successfully maintains stationary conditions within a sloped pipe across various flow conditions, while also preventing spurious oscillations at mixed flow interfaces even when using a pressure wave speed of 1000 m s − 1 .

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.

A low diffusion flux-split scheme for all Mach number flows

A low diffusion version of the Harten–Lax–Leer convective pressure split (HLL-CPS) scheme for resolving the shear layers and the flow features at low Mach numbers is presented here. The low diffusion HLL-CPS scheme is obtained by reconstructing the velocities at the cell interface with the face normal Mach number and a pressure function. Asymptotic analysis of the modified scheme shows a correct scaling of the pressure at low Mach numbers and a significant reduction in numerical dissipation. The robustness of the HLL-CPS scheme for strong shock is improved by reducing the contribution of the contact wave in the vicinity of the shock. The improvement in robustness for strong shock is demonstrated analytically through linear perturbation and matrix stability analyses. A set of numerical test cases are solved to demonstrate the efficacy of the proposed scheme over a wide range of Mach numbers.

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.

An efficient HLL-based scheme for capturing contact-discontinuity in scalar transport by shallow water flow

An efficient numerical scheme is developed for coupled hydrodynamic and scalar transport systems to guarantee the conservation and the positivity-preserving properties for water depth and scalar concentration. A second-order well-balanced positivity-preserving central-upwind scheme based on the finite volume method is adopted to discretize both the Saint-Venant system and the advective fluxes in the scalar transport system. In particular, an anti-diffusion modification is augmented in an ad hoc manner to the Harten–Lax–van Leer approximate Riemann solver for the scalar transport system, with the aim of significantly reducing the numerical diffusion near contact discontinuities. The proposed scheme is validated through seven numerical experiments, wherein the advection and diffusion processes of scalar transport are considered either separately or in combination. Accuracy is measured based on the difference between the numerical results and analytical solutions, and the performance of the proposed scheme is assessed by comparing it with that of the existing scheme. Convergence analysis confirms that the proposed scheme is accurate to the second order. The stationary steady state with wet-dry fronts is simulated, which verifies that the proposed scheme can preserve the exact C-property. The results for the pure diffusion case reveal that the original scheme has limitations in precisely predicting scalar diffusion when the diffusion coefficient is small. By contrast, the proposed scheme accurately approximates scalar diffusion, even at low diffusivity. The results under flow conditions with various Froude and Péclet numbers confirm that the proposed scheme is more accurate than the existing scheme in solving scalar transport problems over a wide range of shallow water flows.

Accurate Numerical Modeling for 1D Open-Channel Flow with Varying Topography

In this paper, an accurate numerical model is presented for one-dimensional open-channel flows with varying topographies; the model is specifically applied to rectangular channels with variable widths. A pressure term is introduced in the shallow water momentum equation to construct a new conservation term, and the resulting non-conservation term is included in the source term to characterize the actual topography changes in the open-channel flow. Based on a Harten Lax and van Leer (HLL) Riemann solver for a homogeneous system, an upwind scheme is introduced into the model in which the flux is determined via randomly selecting a local Riemann solution state. This two-step random choice method enables the scheme to reach second-order accuracy in space. A Runge–Kutta scheme is introduced into the discretization of the source term of the system to achieve second-order accuracy. The proposed model is validated via a selection of steady and transient hydraulic problems with reference solutions. When compared with published experimental results, the predictions of the proposed model show a high degree of accuracy.

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 shock-capturing meshless method for solving the one-dimensional Saint-Venant equations on a highly variable topography

Abstract The Saint-Venant equations are numerically solved to simulate free surface flows in one dimension. A Riemann solver is needed to compute the numerical flux for capturing shocks and flow discontinuities occurring in flow situations such as hydraulic jump, dam-break wave propagation, or bore wave propagation. A Riemann solver that captures shocks and flow discontinuities is not yet reported to be implemented within the framework of a meshless method for solving the Saint-Venant equations. Therefore, a wide range of free surface flow problems cannot be simulated by the available meshless methods. In this study, a shock-capturing meshless method is proposed for simulating one-dimensional (1D) flows on a highly variable topography. The Harten–Lax–van Leer Riemann solver is used for computing the convective flux in the proposed meshless method. Spatial derivatives in the Saint-Venant equations and the reconstruction of conservative variables for flux terms are computed using a weighted least square approximation. The proposed method is tested for various numerically challenging problems and laboratory experiments on different flow regimes. The proposed highly accurate shock-capturing meshless method has the potential to be extended to solve the two-dimensional (2D) shallow water equations without any mesh requirements.

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.

A New Numerical Implementation for Solar Coronal Modeling by an HLL Generalized Riemann Problem Solver

In this paper, we employ a Harten–Lax–van Leer (HLL) generalized Riemann problem (HLL-GRP) solver within the framework of a finite volume method to model 3D solar coronal structures for the first time. Based on the rotational invariance of magnetohydrodynamics (MHD) equations, the HLL-GRP solver is successfully implemented into 3D MHD simulations. To constrain the divergence of the magnetic field, the locally divergence-free weighted-least-squares-based essentially nonoscillatory reconstruction and the properly discretized Godunov–Powell source term are applied. To keep density and pressure positive, a positivity-preserving limiter is added to the reconstructed polynomials of density and pressure. We first test a 3D blast wave problem to preliminarily validate the effectiveness of the proposed scheme on Cartesian structured grid. Then, we further run our code on a six-component grid to numerically study the steady-state coronal structures of Carrington rotation 2218 during the solar minimum phase. A comparison with the two-stage Runge–Kutta scheme is performed for both the 3D blast wave problem and solar coronal problem. Numerical results of large-scale solar coronal structures are basically consistent with the observational characteristics, indicating the robustness of the proposed model.

A finite-volume scheme for modeling compressible magnetohydrodynamic flows at low Mach numbers in stellar interiors

Fully compressible magnetohydrodynamic (MHD) simulations are a fundamental tool for investigating the role of dynamo amplification in the generation of magnetic fields in deep convective layers of stars. The flows that arise in such environments are characterized by low (sonic) Mach numbers (ℳson ≲ 10−2). In these regimes, conventional MHD codes typically show excessive dissipation and tend to be inefficient as the Courant–Friedrichs–Lewy (CFL) constraint on the time step becomes too strict. In this work we present a new method for efficiently simulating MHD flows at low Mach numbers in a space-dependent gravitational potential while still retaining all effects of compressibility. The proposed scheme is implemented in the finite-volume SEVEN-LEAGUE HYDRO (SLH) code, and it makes use of a low-Mach version of the five-wave Harten–Lax–van Leer discontinuities (HLLD) solver to reduce numerical dissipation, an implicit–explicit time discretization technique based on Strang splitting to overcome the overly strict CFL constraint, and a well-balancing method that dramatically reduces the magnitude of spatial discretization errors in strongly stratified setups. The solenoidal constraint on the magnetic field is enforced by using a constrained transport method on a staggered grid. We carry out five verification tests, including the simulation of a small-scale dynamo in a star-like environment at ℳson ~ 10−3. We demonstrate that the proposed scheme can be used to accurately simulate compressible MHD flows in regimes of low Mach numbers and strongly stratified setups even with moderately coarse grids.

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.

Flood Prediction with Two-Dimensional Shallow Water Equations: A Case Study of Tongo-Bassa Watershed in Cameroon

As a result of urbanization, combined with the anthropogenic effects of climate change, natural events such as floods are showing increasingly adverse impacts on human existence. This study proposes a new model, based on shallow water equations, that is able to predict these floods and minimize their impacts. The first-order finite volume method (FVM), the Harten Lax and van Leer (HLL) scheme, and the monotone upwind scheme for conservation laws (MUSCL) are applied in the model. In addition, a virtual boundary cell approach is adopted to achieve a monotonic solution for both interior and boundary cells and flux computations at the boundary cells. The model integrates the infiltration parameters recorded in the area, as well as the Manning coefficient specific to each land-cover type of the catchment region. The results provided were mapped to highlight the potential flood zones and the distribution of water heights throughout the catchment region at any given time, as well as that at the outlet. It has been observed that when standard infiltration and the Manning parameters were selected, the floodable surface increased, as expected, with the increasing rainfall intensity and duration of the simulation. With sufficient infiltration, only a portion of the water tends to stagnate and flow off on the surface toward the outlet. A sensitivity analysis of certain parameters, such as rainfall data and the final infiltration coefficient in the lower watershed of the littoral region, was conducted; the results show that the model simulates well the general character of water flow in the watershed. Finally, the model’s validation using field-collected parameters during the flood of 25 July 2017 and 18 to 22 July 2016 in the Grand Ouaga basin in Burkina reveals Nash–Sutcliffe values of 0.7 and 0.73, respectively.

A physical-constraint-preserving finite volume WENO method for special relativistic hydrodynamics on unstructured meshes

This paper presents a highly robust third-order accurate finite volume weighted essentially non-oscillatory (WENO) method for special relativistic hydrodynamics on unstructured triangular meshes. We rigorously prove that the proposed method is physical-constraint-preserving (PCP), namely, always preserves the positivity of the pressure and the rest-mass density as well as the subluminal constraint on the fluid velocity. The method is built on a highly efficient compact WENO reconstruction on unstructured meshes, a simple PCP limiter, the provably PCP property of the Harten–Lax–van Leer flux, and third-order strong-stability-preserving time discretization. Due to the relativistic effects, the primitive variables (namely, the rest-mass density, velocity, and pressure) are highly nonlinear implicit functions in terms of the conservative variables, making the design and analysis of our method nontrivial. To address the difficulties arising from the strong nonlinearity, we adopt a novel quasilinear technique for the theoretical proof of the PCP property. Three provable convergence-guaranteed iterative algorithms are also introduced for the robust recovery of primitive quantities from admissible conservative variables. We also propose a slight modification to an existing WENO reconstruction to ensure the scaling invariance of the nonlinear weights and thus to accommodate the homogeneity of the evolution operator, leading to the advantages of the modified WENO reconstruction in resolving multi-scale wave structures. Extensive numerical examples are presented to demonstrate the robustness, expected accuracy, and high resolution of the proposed method.

Case studies on simulations of flow-induced vibrations of a cooled circular cylinder: Incompressible flow solver for moving mesh problem

The paper presents an in-house code developed for simulating incompressible flow involving moving boundaries. The solver is based on Arbitrarily Lagrangian-Eulerian formulation with Harten Lax and van Leer with Contact for artificial compressibility (ALE-HLLC-AC) method. In the finite volume formulation, the moving boundary problem is taken care of by the arbitrarily Lagrangian-Eulerian formulation, where radial basis function-based interpolation is employed to move the meshes. The Harten Lax and van Leer with -Contact scheme developed to evaluate the convective fluxes in artificial compressibility formulation. High order accuracy is achieved over an unstructured data structure by developing solution-dependent weighted least squares-based gradient calculations. Flow-induced vibrations caused due to the cyclic lift forces over a cooled circular cylinder are simulated and compared with available literature results. For flow-induced vibration of circular cylinder: fluid-structure interaction problem, the elastically mounted cylinder wall temperature is maintained lower than the flowing fluid. The cylinder vibrates transverse to the direction of fluid flow, and its oscillating movement is simulated using a numerical mass-spring-damper system. Both heat transfer scenarios, i.e., without gravity and gravity acting against the flow direction, are simulated. The unsteady, laminar, incompressible flow analysis is carried out at Richardson number (−1, 0) for the cooled circular cylinder with various reduced velocities Ur (3–20). Prandtl number and Reynolds number are Pr = 7.1 and Re = 150, with the fixed mass ratio of 2 and zero damping coefficient. The validation of the results is presented in the literature. It has been noticed that the thermal boundary layer affects the fluid characteristics. The vibration of the cylinder reached the highest amplitude of cylinder dimension, more dominantly in the range of reduced velocities (4–15), showing a galloping kind of flow nature.

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 High–Order WENO Scheme Based on Different Numerical Fluxes for the Savage–Hutter Equations

The study of rapid free surface granular avalanche flows has attracted much attention in recent years, which is widely modeled using the Savage–Hutter equations. The model is closely related to shallow water equations. We employ a high-order shock-capturing numerical model based on the weighted essential non-oscillatory (WENO) reconstruction method for solving Savage–Hutter equations. Three numerical fluxes, i.e., Lax–Friedrichs (LF), Harten–Lax–van Leer (HLL), and HLL contact (HLLC) numerical fluxes, are considered with the WENO finite volume method and TVD Runge–Kutta time discretization for the Savage–Hutter equations. Numerical examples in 1D and 2D space are presented to compare the resolution of shock waves and free surface capture. The numerical results show that the method proposed provides excellent performance with high accuracy and robustness.