High order finite difference methods with subcell resolution for advection equations with stiff source terms

A new high order finite-difference method utilizing the idea of Harten ENO subcell resolution method is proposed for chemical reactive flows and combustion. In reaction problems, when the reaction time scale is very small, e.g., orders of magnitude smaller than the fluid dynamics time scales, the governing equations will become very stiff. Wrong propagation speed of discontinuity may occur due to the underresolved numerical solution in both space and time. The present proposed method is a modified fractional step method which solves the convection step and reaction step separately. In the convection step, any high order shock-capturing method can be used. In the reaction step, an ODE solver is applied but with the computed flow variables in the shock region modified by the Harten subcell resolution idea. For numerical experiments, a fifth-order finite-difference WENO scheme and its anti-diffusion WENO variant are considered. A wide range of 1D and 2D scalar and Euler system test cases are investigated. Studies indicate that for the considered test cases, the new method maintains high order accuracy in space for smooth flows, and for stiff source terms with discontinuities, it can capture the correct propagation speed of discontinuities in very coarse meshes with reasonable CFL numbers.     

Dynamics of Chapman-Jouguet pulsating detonations with Chain-Branching kinetics: Fickett’s analogue and Euler equations

The present study examines the spatiotemporal nonlinear dynamics of detonations over a wide range of reaction time scales away from the neutral stability region. This is addressed by one-dimensional numerical simulations with chain-branching kinetics. Fickett’s detonation analogue and Euler’s equations were used as evolution equations. A shock-fitting solver is used to reduce CPU time. Up to four thousand five hundred simulations have been carried out. Detailed bifurcation diagrams have been generated to explore the detonation dynamics. For long/intermediate reaction time scales, away from the neutral boundary, the traditional period-doubling cascade to chaos is seen. For square wave detonations, away from the neutral stability, almost periodic oscillations are recorded. This result might have implications for the existence of a characteristic length scale, the cell size, on typical cellular detonations which have a short reaction length.     

A geometric approach to bistable front propagation in scalar reaction–diffusion equations with cut-off

‘Cut-offs’ were introduced to model front propagation in reaction–diffusion systems in which the reaction is effectively deactivated at points where the concentration lies below some threshold. In this article, we investigate the effects of a cut-off on fronts propagating into metastable states in a class of bistable scalar equations. We apply the method of geometric desingularization from dynamical systems theory to calculate explicitly the change in front propagation speed that is induced by the cut-off. We prove that the asymptotics of this correction scales with fractional powers of the cut-off parameter, and we identify the source of these exponents, thus explaining the structure of the resulting expansion. In particular, we show geometrically that the speed of bistable fronts increases in the presence of a cut-off, in agreement with results obtained previously via a variational principle. We first discuss the classical Nagumo equation as a prototypical example of bistable front propagation. Then, we present corresponding results for the (equivalent) cut-off Schlögl equation. Finally, we extend our analysis to a general family of reaction–diffusion equations that support bistable fronts, and we show that knowledge of an explicit front solution to the associated problem without cut-off is necessary for the correction induced by the cut-off to be computable in closed form.     

一种基于特征理论模拟凝聚炸药爆轰的单元中心型Lagrange方法

提出一种数值模拟凝聚炸药爆轰问题的单元中心型Lagrange方法.利用有限体积离散爆轰反应流动方程组,基于双曲型偏微分方程组的特征理论获得离散网格节点的速度与压力,获得的网格节点速度与压力用于更新网格节点位置以及计算网格单元边的数值通量.以这种方式获得的网格节点解是一种"真正多维"的理论解,是一维Godunov格式在二维Riemann问题的推广.有限体积离散得到的爆轰反应流动的半离散系统使用一种显-隐Runge-Kutta格式来离散求解：显式格式处理对流项,隐式格式处理化学反应刚性源项.算例表明,提出的单元中心型Lagrange方法能够较好地模拟凝聚炸药的爆轰反应流动.     

DF激光器燃烧室反应流场的一种新的计算求解方法

 使用了一种新的方法来计算燃烧驱动DF激光器燃烧室多组分化学反应流场。首先引入了元素的分布方程,然后根据化学平衡的热力计算方法求解每一位置的组分及温度。这样,在计算低速流场时,就避免了处理化学反应源项引起的数值计算过程中的方程刚性问题,也减少了组分方程的数量。结果显示,这种燃烧室的注氦方式能够在一定程度上冷却燃烧室壁面。     

Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws

In this article, we propose a new class of finite volume schemes of arbitrary accuracy in space and time for systems of hyperbolic balance laws with stiff source terms. The new class of schemes is based on a three stage procedure. First a high-order WENO reconstruction procedure is applied to the cell averages at the current time level. Second, the temporal evolution of the reconstruction polynomials is computed locally inside each cell using the governing equations. In the original ENO scheme of Harten et al. and in the ADER schemes of Titarev and Toro, this time evolution is achieved via a Taylor series expansion where the time derivatives are computed by repeated differentiation of the governing PDE with respect to space and time, i.e. by applying the so-called Cauchy–Kovalewski procedure. However, this approach is not able to handle stiff source terms. Therefore, we present a new strategy that only replaces the Cauchy–Kovalewski procedure compared to the previously mentioned schemes. For the time-evolution part of the algorithm, we introduce a local space–time discontinuous Galerkin (DG) finite element scheme that is able to handle also stiff source terms. This step is the only part of the algorithm which is locally implicit. The third and last step of the proposed ADER finite volume schemes consists of the standard explicit space–time integration over each control volume, using the local space–time DG solutions at the Gaussian integration points for the intercell fluxes and for the space–time integral over the source term. We will show numerical convergence studies for nonlinear systems in one space dimension with both non-stiff and with very stiff source terms up to sixth order of accuracy in space and time. The application of the new method to a large set of different test cases is shown, in particular the stiff scalar model problem of LeVeque and Yee [R.J. LeVeque, H.C. Yee, A study of numerical methods for hyperbolic conservation laws with stiff source terms, Journal of Computational Physics 86 (1) (1990) 187–210], the relaxation system of Jin and Xin [S. Jin, Z. Xin, The relaxation schemes for systems of conservation laws in arbitrary space dimensions, Communications on Pure and Applied Mathematics 48 (1995) 235–277] and the full compressible Euler equations with stiff friction source terms.     

Collocation Methods for Hyperbolic Partial Differential Equations with Singular Sources

A numerical study is given on the spectral methods and the high order WENO finite difference scheme for the solution of linear and nonlinear hyperbolic partial differential equations with stationary and non-stationary singular sources. The singular source term is represented by the $δ$-function. For the approximation of the $δ$-function, the direct projection method is used that was proposed in [6]. The $δ$-function is constructed in a consistent way to the derivative operator. Nonlinear sine-Gordon equation with a stationary singular source was solved with the Chebyshev collocation method. The $δ$-function with the spectral method is highly oscillatory but yields good results with small number of collocation points. The results are compared with those computed by the second order finite difference method. In modeling general hyperbolic equations with a non-stationary singular source, however, the solution of the linear scalar wave equation with the non-stationary singular source using the direct projection method yields non-physical oscillations for both the spectral method and the WENO scheme. The numerical artifacts arising when the non-stationary singular source term is considered on the discrete grids are explained.     

Pulsating one-dimensional detonations in hydrogen–air mixtures

The propagation of one-dimensional detonations in hydrogen–air mixtures is investigated numerically by solving the one-dimensional Euler equations with detailed finite-rate chemistry. The numerical method is based on a second-order spatially accurate total-variation-diminishing scheme and a point implicit time marching algorithm. The hydrogen–air combustion is modelled with a 9-species, 19-step reaction mechanism. A multi-level, dynamically adaptive grid is utilized, in order to resolve the structure of the detonation. Parametric studies for an equivalence ratio range of 0.4–2.0, initial pressure range of 0.2–0.8 bar and different degrees of detonation overdrive demonstrate that the detonation is unstable for low degrees of overdrive, but the dynamics of wave propagation varies with fuel–air equivalence ratio and pressure. For equivalence ratios less than approximately 1.2 and for all pressures, the detonation exhibits a short-period oscillatory mode, characterized by high-frequency, low-amplitude waves. Richer mixtures exhibit a period-doubled bifurcation that depends on the initial pressure. Parametric studies over a degree of overdrive range of 1.0–1.2 for stoichiometric mixtures at 0.42 bar initial pressure indicate that stable detonation wave propagation is obtained at the high end of this range. For degrees of overdrive close to one, the detonation wave exhibits a low-frequency mode characterized by large fluctuations in the detonation wave speed. The McVey–Toong short-period wave-interaction theory is in qualitative agreement with the numerical simulations; however, the frequencies obtained from their theory are much higher, especially for near-stoichiometric mixtures at high pressure. Modification of this theory to account for the finite heat-release time significantly improves agreement with the numerically computed frequency over the entire equivalence ratio and pressure ranges.     

Evaluation of deconvolution modelling applied to numerical combustion

A possible modelling approach in the large eddy simulation (LES) of reactive flows is to deconvolve resolved scalars. Indeed, by inverting the LES filter, scalars such as mass fractions are reconstructed. This information can be used to close budget terms of filtered species balance equations, such as the filtered reaction rate. Being ill-posed in the mathematical sense, the problem is very sensitive to any numerical perturbation. The objective of the present study is to assess the ability of this kind of methodology to capture the chemical structure of premixed flames. For that purpose, three deconvolution methods are tested on a one-dimensional filtered laminar premixed flame configuration: the approximate deconvolution method based on Van Cittert iterative deconvolution, a Taylor decomposition-based method, and the regularised deconvolution method based on the minimisation of a quadratic criterion. These methods are then extended to the reconstruction of subgrid scale profiles. Two methodologies are proposed: the first one relies on subgrid scale interpolation of deconvolved profiles and the second uses parametric functions to describe small scales. Conducted tests analyse the ability of the method to capture the chemical filtered flame structure and front propagation speed. Results show that the deconvolution model should include information about small scales in order to regularise the filter inversion. a priori and a posteriori tests showed that the filtered flame propagation speed and structure cannot be captured if the filter size is too large.     

Very high order PNPM schemes on unstructured meshes for the resistive relativistic MHD equations

In this paper we propose the first better than second order accurate method in space and time for the numerical solution of the resistive relativistic magnetohydrodynamics (RRMHD) equations on unstructured meshes in multiple space dimensions. The nonlinear system under consideration is purely hyperbolic and contains a source term, the one for the evolution of the electric field, that becomes stiff for low values of the resistivity.     

Weak solutions for partial differential equations with source terms: Application to the shallow water equations

Weak solutions of problems with m equations with source terms are proposed using an augmented Riemann solver defined by m + 1 states instead of increasing the number of involved equations. These weak solutions use propagating jump discontinuities connecting the m + 1 states to approximate the Riemann solution. The average of the propagated waves in the computational cell leads to a reinterpretation of the Roe’s approach and in the upwind treatment of the source term of Vázquez-Cendón. It is derived that the numerical scheme can not be formulated evaluating the physical flux function at the position of the initial discontinuities, as usually done in the homogeneous case. Positivity requirements over the values of the intermediate states are the only way to control the global stability of the method. Also it is shown that the definition of well-balanced equilibrium in trivial cases is not sufficient to provide correct results: it is necessary to provide discrete evaluations of the source term that ensure energy dissipating solutions when demanded. The one and two dimensional shallow water equations with source terms due to the bottom topography and friction are presented as case study. The stability region is shown to differ from the one defined for the case without source terms, and it can be derived that the appearance of negative values of the thickness of the water layer in the proximity of the wet/dry front is a particular case, of the wet/wet fronts. The consequence is a severe reduction in the magnitude of the allowable time step size if compared with the one obtained for the homogeneous case. Starting from this result, 1D and 2D numerical schemes are developed for both quadrilateral and triangular grids, enforcing conservation and positivity over the solution, allowing computationally efficient simulations by means of a reconstruction technique for the inner states of the weak solution that allows a recovery of the time step size.     

Augmented Riemann solvers for the shallow water equations over variable topography with steady states and inundation

We present a class of augmented approximate Riemann solvers for the shallow water equations in the presence of a variable bottom surface. These belong to the class of simple approximate solvers that use a set of propagating jump discontinuities, or waves, to approximate the true Riemann solution. Typically, a simple solver for a system of m conservation laws uses m such discontinuities. We present a four wave solver for use with the the shallow water equations—a system of two equations in one dimension. The solver is based on a decomposition of an augmented solution vector—the depth, momentum as well as momentum flux and bottom surface. By decomposing these four variables into four waves the solver is endowed with several desirable properties simultaneously. This solver is well-balanced: it maintains a large class of steady states by the use of a properly defined steady state wave—a stationary jump discontinuity in the Riemann solution that acts as a source term. The form of this wave is introduced and described in detail. The solver also maintains depth non-negativity and extends naturally to Riemann problems with an initial dry state. These are important properties for applications with steady states and inundation, such as tsunami and flood modeling. Implementing the solver with LeVeque’s wave propagation algorithm [R.J. LeVeque, Wave propagation algorithms for multi-dimensional hyperbolic systems, J. Comput. Phys. 131 (1997) 327–335] is also described. Several numerical simulations are shown, including a test problem for tsunami modeling.     

A high-order discontinuous Galerkin method for wave propagation through coupled elastic–acoustic media

We introduce a high-order discontinuous Galerkin (dG) scheme for the numerical solution of three-dimensional (3D) wave propagation problems in coupled elastic–acoustic media. A velocity–strain formulation is used, which allows for the solution of the acoustic and elastic wave equations within the same unified framework. Careful attention is directed at the derivation of a numerical flux that preserves high-order accuracy in the presence of material discontinuities, including elastic–acoustic interfaces. Explicit expressions for the 3D upwind numerical flux, derived as an exact solution for the relevant Riemann problem, are provided. The method supports h-non-conforming meshes, which are particularly effective at allowing local adaptation of the mesh size to resolve strong contrasts in the local wavelength, as well as dynamic adaptivity to track solution features. The use of high-order elements controls numerical dispersion, enabling propagation over many wave periods. We prove consistency and stability of the proposed dG scheme. To study the numerical accuracy and convergence of the proposed method, we compare against analytical solutions for wave propagation problems with interfaces, including Rayleigh, Lamb, Scholte, and Stoneley waves as well as plane waves impinging on an elastic–acoustic interface. Spectral rates of convergence are demonstrated for these problems, which include a non-conforming mesh case. Finally, we present scalability results for a parallel implementation of the proposed high-order dG scheme for large-scale seismic wave propagation in a simplified earth model, demonstrating high parallel efficiency for strong scaling to the full size of the Jaguar Cray XT5 supercomputer.     

Depolarization of the electromagnetic field scattered by a three-dimensional large-scale irregularity of the lower ionosphere

This paper develops analytical and numerical methods for the solution of three-dimensional problems of radio wave propagation. We consider a three-dimensional vector problem for the electromagnetic field of a vertical electric dipole in a planar Earth-ionosphere waveguide with a largescale local irregularity of negative characteristics at the anisotropic ionospheric boundary. The field components at the boundary surfaces obey the Leontovich boundary conditions. The problem is reduced to a system of two-dimensional integral equations taking into account the overexcitation and depolarization of the field scattered by the irregularity. Using asymptotic (with respect to the parameter kr≫1, where r is the distance from the source or receiver to the nearest point of the irregularity, k=2π/λ, and λ is the radio wavelength) integration over the direction perpendicular to the ray path, we transform this system to one-dimensional integral equations where integration contours represent the geometric contour of the irregularity. The system is numerically solved in the diagonal approximation, combining direct inversion of the Volterra integral operator and subsequent iterations. The proposed numerical algorithm reduces the computer time required for the solution of this problem and is applicable for studying both small-scale and large-scale irregularities. We obtained novel estimates for the field components that are not excited by the source but result entirely from scattering by the sample three-dimensional ionospheric irregularity.     

Operator Splitting Implicit Integration Factor Methods for Stiff Reaction-Diffusion-Advection Systems

For reaction-diffusion-advection equations, the stiffness from the reaction and diffusion terms often requires very restricted time step size, while the nonlinear advection term may lead to a sharp gradient in localized spatial regions. It is challenging to design numerical methods that can efficiently handle both difficulties. For reaction-diffusion systems with both stiff reaction and diffusion terms, implicit integration factor (IIF) method and its higher dimensional analog compact IIF (cIIF) serve as an efficient class of time-stepping methods, and their second order version is linearly unconditionally stable. For nonlinear hyperbolic equations, weighted essentially non-oscillatory (WENO) methods are a class of schemes with a uniformly high-order of accuracy in smooth regions of the solution, which can also resolve the sharp gradient in an accurate and essentially non-oscillatory fashion. In this paper, we couple IIF/cIIF with WENO methods using the operator splitting approach to solve reaction-diffusion-advection equations. In particular, we apply the IIF/cIIF method to the stiff reaction and diffusion terms and the WENO method to the advection term in two different splitting sequences. Calculation of local truncation error and direct numerical simulations for both splitting approaches show the second order accuracy of the splitting method, and linear stability analysis and direct comparison with other approaches reveals excellent efficiency and stability properties. Applications of the splitting approach to two biological systems demonstrate that the overall method is accurate and efficient, and the splitting sequence consisting of two reaction-diffusion steps is more desirable than the one consisting of two advection steps, because CWC exhibits better accuracy and stability.     

守恒律方程组的大时间步长叠波格式

大时间步长叠波格式最初思想为LeVeque提出的大时间步长Godunov格式,通过叠加间断分解发出的强波来构造数值格式.原方法只给出了间断强波的穿越叠加方法,文章对其进行了完善,并推广到多维.针对膨胀波提出了一种网格单元分解法可以自动满足熵条件,避免出现非物理解.给出了格式的具体计算公式,并用单个守恒律方程、一维/多维Euler方程组进行了数值计算.计算结果表明,新格式除了可以采用大时间步长的优点外,在一定范围内随CFL数增加其耗散反而更低,因而对激波接触间断膨胀波的分辨率更高.     

A large time step 1D upwind explicit scheme (CFL > 1): Application to shallow water equations

It is possible to relax the Courant–Friedrichs–Lewy condition over the time step when using explicit schemes. This method, proposed by Leveque, provides accurate and correct solutions of non-sonic shocks. Rarefactions need some adjustments which are explored in the present work with scalar equation and systems of equations. The non-conservative terms that appear in systems of conservation laws introduce an extra difficulty in practical application. The way to deal with source terms is incorporated into the proposed procedure. The boundary treatment is analysed and a reflection wave technique is considered. In presence of strong discontinuities or important source terms, a strategy is proposed to control the stability of the method allowing the largest time step possible. The performance of the above scheme is evaluated to solve the homogeneous shallow water equations and the shallow water equations with source terms.     

A review of direct numerical simulations of astrophysical detonations and their implications

Multi-dimensional direct numerical simulations (DNS) of astrophysical detonations in degenerate matter have revealed that the nuclear burning is typically characterized by cellular structure caused by transverse instabilities in the detonation front. Type Ia supernova modelers often use onedimensional DNS of detonations as inputs or constraints for their whole star simulations.While these one-dimensional studies are useful tools, the true nature of the detonation is multi-dimensional. The multi-dimensional structure of the burning influences the speed, stability, and the composition of the detonation and its burning products, and therefore, could have an impact on the spectra of Type Ia supernovae. Considerable effort has been expended modeling Type Ia supernovae at densities above 1×10⁷ g·cm^-3 where the complexities of turbulent burning dominate the flame propagation. However, most full star models turn the nuclear burning schemes off when the density falls below 1×10⁷ g·cm^-3 and distributed burning begins. The deflagration to detonation transition (DDT) is believed to occur at just these densities and consequently they are the densities important for studying the properties of the subsequent detonation. This work will review the status of DNS studies of detonations and their possible implications for Type Ia supernova models. It will cover the development of Detonation theory from the first simple Chapman–Jouguet (CJ) detonation models to the current models based on the time-dependent, compressible, reactive flow Euler equations of fluid dynamics.     

迎风紧致格式求解Hamilton-Jacobi方程

基于Hamilton-Jacobi(H-J)方程和双曲型守恒律之间的关系,将三阶和五阶迎风紧致格式推广应用于求解H-J方程,建立了高精度的H-J方程求解方法.给出了一维和二维典型数值算例的计算结果,其中包括一个平面激波作用下的Richtmyer Meshkov界面不稳定性问题.数值试验表明,在解的光滑区域该方法具有高精度,而在导数不连续的不光滑区域也获得了比较好的分辨效果.相比于同阶精度的WENO格式,本方法具有更小的数值耗散,从而有利于多尺度复杂流动的模拟中H-J方程的求解.     

Method of adaptive artificial viscosity

A new finite-difference method for the numerical solution of gas dynamics equations is proposed. This method is a uniform monotonous finite-difference scheme of second-order approximation on time and space outside of domains of shock and compression waves. This method is based on inputting adaptive artificial viscosity (AAV) into gas dynamics equations. In this paper, this method is analyzed for 2D geometry. The testing computations of the movement of contact discontinuities and shock waves and the breakup of discontinuities are demonstrated.