Shock-capturing Research Articles

A Class of New High-order Finite-Volume TENO Schemes for Hyperbolic Conservation Laws with Unstructured Meshes

The recently proposed high-order TENO scheme [Fu et al., Journal of Computational Physics, 305, pp.333-359] has shown great potential in predicting complex fluids owing to the novel weighting strategy, which ensures the high-order accuracy, the low numerical dissipation, and the sharp shock-capturing capability. However, the applications are still restricted to simple geometries with Cartesian or curvilinear meshes. In this work, a new class of high-order shock-capturing TENO schemes for unstructured meshes are proposed. Similar to the standard TENO schemes and some variants of WENO schemes, the candidate stencils include one large stencil and several small third-order stencils. Following a strong scale-separation procedure, a tailored novel ENO-like stencil selection strategy is proposed such that the high-order accuracy is restored in smooth regions by selecting the candidate reconstruction on the large stencil while the ENO property is enforced near discontinuities by adopting the candidate reconstruction from smooth small stencils. The nonsmooth stencils containing genuine discontinuities are explicitly excluded from the final reconstruction, leading to excellent numerical stability. Different from the WENO concept, such unique sharp stencil selection retains the low numerical dissipation without sacrificing the shock-capturing capability. The newly proposed framework enables arbitrarily high-order TENO reconstructions on unstructured meshes. For conceptual verification, the TENO schemes with third- to sixth-order accuracy are constructed. Without parameter tuning case by case, the performance of the proposed TENO schemes is demonstrated by examining a set of benchmark cases with broadband flow length scales.

Large eddy simulation of two separated hypersonic shock/turbulent boundary layer interactions

We present a detailed statistical description of turbulence in the high-fidelity Large Eddy Simulation (LES) of hypersonic shock-and-turbulent boundary layer interactions with strong separation. Our high-fidelity LES method uses a 4th-order, bandwidth-optimized Weighted Essentially Non-Oscillatory (WENO) scheme combined with a robust shock-capturing method and dynamic mixed model formulations for the closure of subgrid scale terms. Our mixed models, which contain a scale-similar term in addition to the purely dissipative one, are shown to accurately and significantly improve separation prediction compared to purely dissipative closures, such as the dynamic eddy viscosity formulation.

A Wave-Targeted Essentially Non-Oscillatory 3D Shock-Capturing Scheme for Breaking Wave Simulation

A new three-dimensional high-order shock-capturing model for the numerical simulation of breaking waves is proposed. The proposed model is based on an integral contravariant form of the Navier–Stokes equations in a time-dependent generalized curvilinear coordinate system. Such an integral contravariant form of the equations of motion is numerically integrated by a new conservative numerical scheme that is based on three elements of originality: the time evolution of the state of the system is carried out using a predictor–corrector method in which exclusively the conserved variables are used; the point values of the conserved variables on the cell face of the computational grid are obtained using an original high-order reconstruction procedure called a wave-targeted essentially non-oscillatory scheme; the time evolution of the discontinuity on the cell faces is calculated using an exact Riemann solver. The proposed model is validated by numerically reproducing several experimental tests of breaking waves on computational grids that are significantly coarser than those used in the literature to validate the existing 3D shock-capturing models. The results obtained with the proposed model are also compared with those obtained with a previously published model, which is based on second-order total variation diminishing reconstructions and an approximate Riemann solver usually adopted in the existing 3D shock-capturing models. Through the above comparison, the main drawbacks of the existing 3D shock-capturing models and the ability of the proposed model to simulate breaking waves and wave-induced currents are shown. The proposed 3D model is able to correctly simulate the wave height increase in the shoaling zone and to effectively predict the location of the wave breaking point, the maximum wave height, and the wave height decay in the surf zone. The validated model is applied to the simulation of the interaction between breaking waves and an emerged breakwater. The numerical results show that the proposed model is able to simulate both the large-scale circulation patterns downstream of the barrier and the onset of quasi-periodic vortex structures close to the edge of the barrier.

Mesh-free simulation of two-phase fluid flow in porous media based on the shock-fitting method

The physical interactions between two immiscible fluids in porous media may lead to the formation of a shock front at the interface between the two fluid phases. Numerical methods for modeling fluid flow with shocks can be classified into two categories: the shock-capturing and the shock-fitting methods. The shock-fitting method can in general provide more accurate solution of the shock front than the shock-capturing method, but suffers from the extra complexity in accommodating the moving shock front inside the background mesh. In this study, we explore the possibility of integrating the mesh-free methods for solving partial-differential equations with the shock-fitting method. In this integrated method, which we call “mesh-free shock-fitting (MFSF)”, the nodes needed for mesh-free calculations of spatial derivatives are generated on the fly to adapt to the moving shock front. We demonstrate the implementation of MFSF using numerical experiments in one and two spatial dimensions. Preliminary results show that MFSF has the potential to provide more accurate solutions with lower computational cost than conventional shock-capturing methods and can simplify many of the operations in shock-fitting methods with unstructured meshes.

A new paradigm of dissipation-adjustable, multi-scale resolving schemes for compressible flows

The scale-resolving simulation of high speed compressible flow through direct numerical simulation (DNS) or large eddy simulation (LES) requires shock-capturing schemes to be more accurate for resolving broadband turbulence and robust for capturing strong shock waves. In this work, we develop a new paradigm of dissipation-adjustable, shock capturing scheme to resolve multi-scale flow structures in high speed compressible flow. The new scheme employs a polynomial of n-degree and non-polynomial THINC (Tangent of Hyperbola for INterface Capturing) functions of m-level steepness as reconstruction candidates. These reconstruction candidates are denoted as PnTm. From these candidates, the piecewise reconstruction function is selected through the boundary variation diminishing (BVD) algorithm. Unlike other shock-capturing techniques, the BVD algorithm effectively suppresses numerical oscillations without introducing excess numerical dissipation. Then, an adjustable dissipation (AD) algorithm is designed for scale-resolving simulations. This novel paradigm of shock-capturing scheme is named as PnTm−BVD−AD. The proposed PnTm−BVD−AD scheme has following desirable properties. First, it can capture large-scale discontinuous structures such as strong shock waves without obvious non-physical oscillations while resolving sharp contact, material interface and shear layer. Secondly, the numerical dissipation property of PnTm−BVD−AD can be effectively adjusted between n+1 order upwind-biased scheme and non-dissipative n+2 order central scheme through a simple tunable parameter λ. Thirdly, with λ=0.5 the scheme can recover to n+2 order non-dissipative central interpolation for smooth solution over all wavenumber, which is preferable for solving small-scale structures in DNS as well as resolvable-scale in explicit LES. Finally, the under-resolved small-scale can be solved with the dissipation adjustable algorithm through the so-called implicit LES (ILES) approach. Through simulating benchmark tests involving multi-scale flow structures and comparing with other central-upwind schemes, the superiority of the proposed scheme is evident. For instance, the simulation results of the supersonic planar jet show that PnTm−BVD−AD schemes can achieve competitive results as Pn+2Tm−BVD schemes which utilize a higher degree of reconstruction polynomial. Thus, in comparison with the previous work, the proposed PnTm−BVD−AD schemes have the benefit of a more compact stencil and lower cost. In summary, this work provides an alternative scheme for solving multi-scale problems in high speed compressible flows.

Protostellar discs subject to infall: a one-dimensional inviscid model and comparison with ALMA observations

ABSTRACT A new one-dimensional, inviscid, and vertically integrated disc model with prescribed infall is presented. The flow is computed using a second-order shock-capturing scheme. Included are vertical infall, radial infall at the outer radial boundary, radiative cooling, stellar irradiation, and heat addition at the disc-surface shock. Simulation parameters are chosen to target the L1527 IRS disc which has been observed using Atacama Large Millimeter Array (ALMA). The results give an outer envelope of radial infall and uϕ ∝ 1/r which encounters a radial shock at rshock ∼ 1.5 × the centrifugal radius (rc) across which the radial velocity is greatly reduced and the gas temperature rises from a pre-shock value of ≈25 to ≈180 K over a spatially thin region calculated using a separate shock structure code. At rc, the azimuthal velocity uϕ transitions from being ∝ 1/r to being nearly Keplerian. These results qualitatively agree with recent ALMA observations which indicate a radial shock where SO is sublimated as well as a transition from a uϕ ∼ 1/r region to a Keplerian inner disc. However, in one set of observations, the observed position-velocity map of cyclic-C3H2, together with a certain ballistic maximum velocity relation suggests that the radial shock coincides with a ballistic centrifugal barrier, which places the shock at rshock = 0.5rc, i.e. inward of rc, rather than outward as given by our simulations. It is argued that radial velocity plots from previous magnetic rotating-collapse simulations also indicate that the radial shock is located outward of rc. The discrepancy with observations is analysed and discussed, but remains unresolved.

Artificial Compressibility with Riemann Solvers: Convergence of Limiters on Unstructured Meshes

Free-surface flows and other variable density incompressible flows have numerous important applications in engineering.One way such flows can be modelled is to extend established numerical methods for compressible flows to incompressible flows using the method of artificial compressibility. Artificial compressibility introduces a pseudo-time derivative for pressure and, in each real-time step, the solution advances in pseudo-time until convergence to an incompressible limit - a fundamentally different approach than SIMPLE, PISO, and PIMPLE, the standard methods used in OpenFOAM. Although the artificial compressibility method is widespread in the literature, its application to free-surface flows is not. In this paper, we apply the method to variable density flows on 3D unstructured meshes for the first time, implementing a Godunov-type scheme with MUSCL reconstruction and Riemann solvers, where the free surface gets captured automatically by the contact wave in the Riemann solver. The critical problem in this implementation lies in the slope limiters used in the MUSCL reconstruction step. It is well-known that slope limiters can inhibit convergence to steady state on unstructured meshes; the problem is exacerbated here as convergence in pseudo-time is required not just once, but at every real-time step. We compare the limited gradient schemes included in OpenFOAM with an improved limiter from the literature, testing the solver against dam-break and hydrostatic pressure benchmarks. This work opens OpenFOAM up to the method of artificial compressibility, breaking the mould of PIMPLE and harnessing high-resolution shock-capturing schemes that are easier to parallelise.

No Evidence of Kinetic Screening in Simulations of Merging Binary Neutron Stars beyond General Relativity.

We have conducted fully relativistic simulations in a class of scalar-tensor theories with derivative self-interactions and screening of local scales. By using high-resolution shock-capturing methods and a nonvanishing shift vector, we have managed to avoid issues plaguing similar attempts in the past. We have first confirmed recent results by ourselves in spherical symmetry, obtained with an approximate approach and pointing at a partial breakdown of the screening in black-hole collapse. Then, we considered the late inspiral and merger of binary neutron stars. We found that screening tends to suppress the (subdominant) dipole scalar emission, but not the (dominant) quadrupole scalar mode. Our results point at quadrupole scalar signals as large as (or even larger than) in Fierz-Jordan-Brans-Dicke theories with the same conformal coupling, for strong-coupling scales in the MeV range that we can simulate.

An application of Gaussian process modeling for high-order accurate adaptive mesh refinement prolongation

We present a new polynomial-free prolongation scheme for Adaptive Mesh Refinement (AMR) simulations of compressible and incompressible computational fluid dynamics. The new method is constructed using a multi-dimensional kernel-based Gaussian Process (GP) prolongation model. The formulation for this scheme was inspired by the GP methods introduced by A. Reyes et al. (A New Class of High-Order Methods for Fluid Dynamics Simulation using Gaussian Process Modeling, Journal of Scientific Computing, 76 (2017), 443-480; A variable high-order shock-capturing finite difference method with GP-WENO, Journal of Computational Physics, 381 (2019), 189-217). In this paper, we extend the previous GP interpolations and reconstructions to a new GP-based AMR prolongation method that delivers a high-order accurate prolongation of data from coarse to fine grids on AMR grid hierarchies. In compressible flow simulations special care is necessary to handle shocks and discontinuities in a stable manner. To meet this, we utilize the shock handling strategy using the GP-based smoothness indicators developed in the previous GP work by A. Reyes et al. We demonstrate the efficacy of the GP-AMR method in a series of testsuite problems using the AMReX library, in which the GP-AMR method has been implemented.

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.

Dynamical evolution of dark matter admixed neutron stars

We dynamically evolve for the first time dark matter admixed neutron stars with fermionic dark matter. These systems are mixtures of the ordinary nuclear matter of a neutron star and dark matter. To perform our dynamical evolutions, we derive the equations of motion, in conservation form, for spherically symmetric systems with an arbitrary number of perfect fluids. Using finite volume and high-resolution shock-capturing methods, we dynamically evolve the two-fluid case, with the first fluid modeling ordinary matter and the second fluid modeling dark matter. We use our dynamical solutions to study nonlinear stability, radial oscillation frequencies, and a dynamical formation process.

On the Solution Accuracy Downstream of Shocks When Using Shock-Capturing Methods. I. Sources of Errors in One-Dimensional Problems

On the Solution Accuracy Downstream of Shocks When Using Shock-Capturing Methods. I. Sources of Errors in One-Dimensional Problems

Gp-Mood: A Positive-Preserving High-Order Finite Volume Method for Hyperbolic Conservation Laws

We present an a posteriori shock-capturing finite volume method algorithm called GP-MOOD that solves a compressible hyperbolic conservative system at high-order solution accuracy (e.g., third-, fifth-, and seventh-order) in multiple spatial dimensions. The GP-MOOD method combines two methodologies, the polynomial-free spatial reconstruction methods of GP (Gaussian Process) and the a posteriori detection algorithms of MOOD (Multidimensional Optimal Order Detection). The spatial approximation of our GP-MOOD method uses GP's unlimited spatial reconstruction that builds upon our previous studies on GP reported in Reyes et al., Journal of Scientific Computing, 76 (2017) and Journal of Computational Physics, 381 (2019). This paper focuses on extending GP's flexible variability of spatial accuracy to an a posteriori detection formalism based on the MOOD approach. We show that GP's polynomial-free reconstruction provides a seamless pathway to the MOOD's order cascading formalism by utilizing GP's novel property of variable (2R+1)th-order spatial accuracy on a multidimensional GP stencil defined by the GP radius R, whose size is smaller than that of the standard polynomial MOOD methods. The resulting GP-MOOD method is a positivity-preserving method. We examine the numerical stability and accuracy of GP-MOOD on smooth and discontinuous flows in multiple spatial dimensions without resorting to any conventional, computationally expensive a priori nonlinear limiting mechanism to maintain numerical stability.

A Priori Subcell Limiting Approach for the FR/CPR Method on Unstructured Meshes

A priori subcell limiting approach is developed for high-order flux reconstruction/correction procedure via reconstruction (FR/CPR) on two-dimensional unstructured quadrilateral meshes. Firstly, a modified indicator based on modal energy coefficients is proposed to detect troubled cells. Then, troubled cells are decomposed into nonuniform subcells and each subcell has one solution point. A second-order finite difference shock-capturing scheme based on nonuniform nonlinear weighted (NNW) interpolation is constructed to calculate troubled cells while smooth cells are calculated by the CPR method. Numerical investigations show that the subcell limiting strategy on unstructured quadrilateral meshes is robust in shock-capturing.

A novel high-order low-dissipation TENO-THINC scheme for hyperbolic conservation laws

The high-order low-dissipation shock-capturing scheme is important for compressible fluid simulations, in particular for cases where both shock and turbulence present. However, the competing requirements for resolving the small-scale turbulence structures with low dissipation and capturing the discontinuities sharply render it non-trivial for numerical development. In this work, based on a novel parameter-free discontinuity-detection criterion, a new shock-capturing framework is proposed by combining the standard TENO (targeted essentially non-oscillatory) scheme for smooth regions with the non-polynomial based THINC (tangent of hyperbola for INterface capturing) reconstruction for non-smooth discontinuities. The resulting hybrid scheme retains the low-dissipation property of TENO scheme for smooth flow scales and the discontinuity-resolving capability of THINC reconstruction for shock and contact waves. A set of benchmark simulations involving a wide range of length scales demonstrates that the new TENO-THINC scheme performs significantly better than the standard TENO scheme without the necessity of parameter tuning case by case, and thus is promising for more complex compressible flow predictions where the excessive numerical dissipation of existing schemes prevents the targeted flow structures from being properly resolved. Moreover, the present results serve as the state-of-the-art references for future numerical method development.

Symmetry-preserving enforcement of low-dissipation method based on boundary variation diminishing principle

A class of high-order shock-capturing schemes, PnTm-BVD (Deng et al., 2019; Deng et al., 2020) schemes, have been devised to solve the Euler equations with substantially reduced numerical dissipation, which enable high-resolution simulations to resolve flow structures of wider range scales. In such simulations with low dissipation, errors of round-off level might grow and contaminate the numerical solutions. A typical example of such problems is the loss of symmetry in the numerical solutions for physical problems of symmetric configurations even if the schemes are mathematically in line with the symmetry rules. In this study, the mechanisms of symmetry-breaking in a finite volume framework with the P4T2-BVD reconstruction scheme are thoroughly examined. Particular attention has been paid to remove the possible causes due to the lack of associativity in floating-point arithmetic which is associated with round-off errors. Modifications and new techniques are proposed to completely remove the possible causes for symmetry breaking in different components of the P4T2-BVD finite volume solver. Benchmark tests that have symmetric solution structures are used to verify the proposed methods. The numerical results demonstrate the perfect symmetric solution structures.

High-order numerical scheme for compressible multi-component real gas flows using an extension of the Roe approximate Riemann solver and specific Monotonicity-Preserving constraints

The purpose of this paper is to develop a high-order shock-capturing scheme capable of predicting flows where shock waves with high-temperature jumps interact with multi-component real gas mixtures, assuming a local thermodynamic equilibrium. We first propose a generalization of the Roe solver for distinct species with non-ideal thermodynamic properties that relies on the original method proposed by Vinokur & Montagné [1]. This method uses an approximation of compressibility factors to estimate a coherent value of the speed of sound at the Roe averaged state.This Roe averaged state is introduced in the One-Step Monotonicity-Preserving (OSMP) scheme, originally developed by Daru and Tenaud [2], to obtain an extension to the high-order with Lax-Wendroff procedure adequate for dealing with non-ideal gas flows. To avoid thermodynamic inconsistencies in the evolution of the Roe average state over a large stencil, we propose to reformulate the discrete total energy flux of the initial solver. This new formulation uses a combination of Riemann invariants related to the species mass fractions and avoids the influence of the independent values of the compressibility factors in the total energy flux computation. An additional M-P constraint on this new combination allows dealing with discontinuities. Based on the averaged speed of sound estimated by our proposed extension of the Vinokur & Montagné method, we demonstrate that this new formulation is equivalent to selecting a new combination of compressibility factors that completely fulfill the jump relationships of the Riemann problem.To properly capture discontinuities while optimizing the number of numerical cells, the new high-order OSMP scheme is combined with an Adaptive Multiresolution [3] procedure to automatically refine grid in regions where steep gradients occur and coarsen grid elsewhere. The order of the numerical method is evaluated on the convection of density and mass fraction waves. Its capability of capturing discontinuities is validated on a 1-D shock tube problem with a mixture of Nitrogen, Oxygen and dense refrigerant R22 gases. We show that smooth solutions, as well as discontinuities, are recovered with high accuracy. The 2-D interaction between a shock wave in Air with a cylindrical bubble initially filled with dense refrigerant R22 gas is also considered. Present results compare very well with both a recent fully resolved numerical solution of ideal gases and experimental results obtained with real gases. Compared to ideal gas solutions corresponding to calorically perfect gas, drastic changes are recorded on the predicted temperature and the bubble flow patterns that fully justify the use of relevant thermodynamics and the proposed numerical method to account for real gas properties.

Enhanced fifth order WENO shock-capturing schemes with deep learning

In this paper we enhance the well-known fifth order WENO shock-capturing scheme by using deep learning techniques. This fine-tuning of an existing algorithm is implemented by training a rather small neural network to modify the smoothness indicators of the WENO scheme in order to improve the numerical results especially at discontinuities. In our approach no further post-processing is needed to ensure the consistency of the method. Moreover, the formal order of accuracy of the resulting scheme can be proven.We demonstrate our findings with the inviscid Burgers’ equation, the Buckley–Leverett equation and the 1-D Euler equations of gas dynamics. Hereby we investigate the classical Sod problem and the Lax problem and show that our novel method outperforms the classical fifth order WENO schemes in simulations where the numerical solution is too diffusive or tends to overshoot at shocks. Finally, the straight-forward extension of the method to two-dimensional problems is included and illustrated using the 2D Burgers’ equation.

A numerical study on near-field pressure fluctuations of symmetrical and anti-symmetrical flapping modes of twin-jet using a high-resolution shock-capturing scheme

Screeching supersonic jets appears at off-design operating conditions and is perceived as an intense tonal noise. In a twin nozzle configuration, mutual interactions between the two jet plumes may occur with various coupling modes developing depending on the operating conditions and lateral distance between the jets. The investigation of the detailed flow behaviors and near-field pressure fluctuations with relevance to the twin jets system, the analysis of the developed instabilities, will enhance understanding of fundamental features associated with jets located close to each other.In the present study, the single jet is considered first to assess the large eddy simulation (LES) approach used and the near-field pressure fluctuation predictions. Based on the validated solver, twin jets are simulated. Two different twin-nozzle configurations having different separation distance or nozzle-to-nozzle centerline spacing are scrutinized for the same Mach number of 1.358. Notably, the twin jets are screeching by the coupling mode for both set-ups; however, the case of closer inter-nozzle distance presents a symmetrical dominant flapping mode, while the other case shows an anti-symmetrical flapping mode. The strength of the pressure fluctuation at the fundamental frequency changes depending on the location of the observer point (upstream or downstream) and the reference plane (twin-jet and normal to the twin-jet plane). The screech tones of the two cases, observable in the upstream region, are significantly different in the normal to the twin-jet plane direction because of the phase difference of fluctuating pressure. However, the first harmonic component remains strong, regardless of the flapping mode. It is also observed that, at the fundamental frequency, the amplitude of the pressure fluctuation at downstream locations is found to be strong in the normal to the twin-jet plane when the symmetrical flapping mode occurs. This feature is also observed in the twin-jet plane in the case of the opposite mode. By analyzing the developed vertical structures and performing correlation analyses of pressure fluctuations along jet shear layers, the periodicity of the flow in the downstream region with relevance to the fundamental frequency is revealed.

Gmunu: paralleled, grid-adaptive, general-relativistic magnetohydrodynamics in curvilinear geometries in dynamical space–times

ABSTRACT We present an update on the General-relativistic multigrid numerical (Gmunu) code, a parallelized, multidimensional curvilinear, general relativistic magnetohydrodynamics code with an efficient non-linear cell-centred multigrid elliptic solver, which is fully coupled with an efficient block-based adaptive mesh refinement module. To date, as described in this paper, Gmunu is able to solve the elliptic metric equations in the conformally flat condition approximation with the multigrid approach and the equations of ideal general-relativistic magnetohydrodynamics by means of high-resolution shock-capturing finite-volume method with reference metric formularised multidimensionally in Cartesian, cylindrical, or spherical geometries. To guarantee the absence of magnetic monopoles during the evolution, we have developed an elliptical divergence cleaning method by using the multigrid solver. In this paper, we present the methodology, full evolution equations and implementation details of Gmunu and its properties and performance in some benchmarking and challenging relativistic magnetohydrodynamics problems.