Stiff Source Terms Research Articles

Explicit one-step time discretizations for discontinuous Galerkin and finite volume schemes based on local predictors

We consider a family of explicit one-step time discretizations for finite volume and discontinuous Galerkin schemes, which is based on a predictor-corrector formulation. The predictor remains local taking into account the time evolution of the data only within the grid cell. Based on a space–time Taylor expansion, this idea is already inherent in the MUSCL finite volume scheme to get second order accuracy in time and was generalized in the context of higher order ENO finite volume schemes. We interpret the space–time Taylor expansion used in this approach as a local predictor and conclude that other space–time approximate solutions of the local Cauchy problem in the grid cell may be applied. Three possibilities are considered in this paper: (1) the classical space–time Taylor expansion, in which time derivatives are obtained from known space-derivatives by the Cauchy–Kovalewsky procedure; (2) a local continuous extension Runge–Kutta scheme and (3) a local space–time Galerkin predictor with a version suitable for stiff source terms. The advantage of the predictor–corrector formulation is that the time evolution is done in one step which establishes optimal locality during the whole time step. This time discretization scheme can be used within all schemes which are based on a piecewise continuous approximation as finite volume schemes, discontinuous Galerkin schemes or the recently proposed reconstructed discontinuous Galerkin or P N P M schemes. The implementation of these approaches is described, advantages and disadvantages of different predictors are discussed and numerical results are shown.

ADER Schemes for Nonlinear Systems of Stiff Advection–Diffusion–Reaction Equations

In this article we extend the high order ADER finite volume schemes introduced for stiff hyperbolic balance laws by Dumbser, Enaux and Toro (J. Comput. Phys. 227:3971---4001, 2008) to nonlinear systems of advection---diffusion---reaction equations with stiff algebraic source terms. We derive a new efficient formulation of the local space-time discontinuous Galerkin predictor using a nodal approach whose interpolation points are tensor-products of Gauss---Legendre quadrature points. Furthermore, we propose a new simple and efficient strategy to compute the initial guess of the locally implicit space-time DG scheme: the Gauss---Legendre points are initialized sequentially in time by a second order accurate MUSCL-type approach for the flux term combined with a Crank---Nicholson method for the stiff source terms. We provide numerical evidence that when starting with this initial guess, the final iterative scheme for the solution of the nonlinear algebraic equations of the local space-time DG predictor method becomes more efficient. We apply our new numerical method to some systems of advection---diffusion---reaction equations with particular emphasis on the asymptotic preserving property for linear model systems and the compressible Navier---Stokes equations with chemical reactions.

Implicit–explicit schemes for flow equations with stiff source terms

In this paper, we design stable and accurate numerical schemes for conservation laws with stiff source terms. A prime example and the main motivation for our study is the reactive Euler equations of gas dynamics. Furthermore, we consider widely studied scalar model equations. We device one-step IMEX (implicit–explicit) schemes for these equations that treats the convection terms explicitly and the source terms implicitly. For the non-linear scalar equation, we use a novel choice of initial data for the resulting Newton solver and obtain correct propagation speeds, even in the difficult case of rarefaction initial data. For the reactive Euler equations, we choose the numerical diffusion suitably in order to obtain correct wave speeds on under-resolved meshes. We prove that our implicit–explicit scheme converges in the scalar case and present a large number of numerical experiments to validate our scheme in both the scalar case as well as the case of reactive Euler equations. Furthermore, we discuss fundamental differences between the reactive Euler equations and the scalar model equation that must be accounted for when designing a scheme.

A hybrid Godunov method for radiation hydrodynamics

From a mathematical perspective, radiation hydrodynamics can be thought of as a system of hyperbolic balance laws with dual multiscale behavior (multiscale behavior associated with the hyperbolic wave speeds as well as multiscale behavior associated with source term relaxation). With this outlook in mind, this paper presents a hybrid Godunov method for one-dimensional radiation hydrodynamics that is uniformly well behaved from the photon free streaming (hyperbolic) limit through the weak equilibrium diffusion (parabolic) limit and to the strong equilibrium diffusion (hyperbolic) limit. Moreover, one finds that the technique preserves certain asymptotic limits. The method incorporates a backward Euler upwinding scheme for the radiation energy density E r and flux F r as well as a modified Godunov scheme for the material density ρ, momentum density m, and energy density E. The backward Euler upwinding scheme is first-order accurate and uses an implicit HLLE flux function to temporally advance the radiation components according to the material flow scale. The modified Godunov scheme is second-order accurate and directly couples stiff source term effects to the hyperbolic structure of the system of balance laws. This Godunov technique is composed of a predictor step that is based on Duhamel’s principle and a corrector step that is based on Picard iteration. The Godunov scheme is explicit on the material flow scale but is unsplit and fully couples matter and radiation without invoking a diffusion-type approximation for radiation hydrodynamics. This technique derives from earlier work by Miniati and Colella (2007) [41]. Numerical tests demonstrate that the method is stable, robust, and accurate across various parameter regimes.

An implementation of the microphysics in full general relativity: a general relativistic neutrino leakage scheme

Performing fully general relativistic simulations taking account of microphysical processes (e.g. weak interactions and neutrino cooling) is one of the long-standing problems in numerical relativity. One of main difficulties in implementation of weak interactions in the general relativistic framework lies in the fact that the characteristic timescale of weak interaction processes (the WP timescale, ) in hot dense matters is much shorter than the dynamical timescale (tdyn). Numerically this means that stiff source terms appear in the equations so that an implicit scheme is in general necessary to stably solve the relevant equations. Otherwise a very short timestep (Δt < twp ≪ tdyn) will be required to solve them explicitly, which is unrealistic in the present computational resources. Furthermore, in the relativistic framework, the Lorentz factor is coupled with the rest mass density and the energy density. The specific enthalpy is also coupled with the momentum. Due to these couplings, it is very complicated to recover the primitive variables and the Lorentz factor from conserved quantities. Consequently, it is very difficult to solve the equations implicitly in the fully general relativistic framework. At the current status, no implicit procedure has been proposed except for the case of the spherical symmetry. Therefore, an approximate explicit procedure is developed in the fully general relativistic framework in this paper as a first implementation of the microphysics toward a more realistic sophisticated model. The procedure is based on the so-called neutrino leakage schemes which are based on the property that the characteristic timescale in which neutrinos leak out of the system (the leakage timescale, tleak) is much longer than the WP timescale. In the previous leakage schemes, however, the problems of the stiff source terms are avoided in an artificial manner. In this paper, I present a detailed neutrino leakage scheme and a simple and stable method for solving the equations explicitly in the fully general relativistic framework. The drawback of the artificial treatment of the stiff source terms is improved. I also perform a test simulation to check the validity of the present method, showing that it works fairly well.

On Stability, Monotonicity, and Construction of Difference Schemes I: Theory

The stability of nonlinear explicit difference schemes with not, in general, open domains of the scheme operators are studied. For the case of path-connected, bounded, and Lipschitz domains, we establish the notion that a multi-level nonlinear explicit scheme is stable iff (if and only if) the corresponding scheme in variations is stable. A new modification of the central Lax-Friedrichs (LxF) scheme is developed to be of the second order accuracy. The modified scheme is based on nonstaggered grids. A monotone piecewise cubic interpolation is used in the central scheme to give an accurate approximation for the model in question. The stability of the modified scheme is investigated. Some versions of the modified scheme are tested on several conservation laws, and the scheme is found to be accurate and robust. As applied to hyperbolic conservation laws with, in general, stiff source terms, it is constructed a second order nonstaggered central scheme based on operator-splitting techniques.

A Micro-Macro Decomposition-Based Asymptotic-Preserving Scheme for the Multispecies Boltzmann Equation

In this paper we extend the micro-macro decomposition-based asymptotic-preserving scheme developed by Bennoune, Lemou, and Mieussens [J. Comput. Phys., 227 (2008), pp. 3781-3803] for the single species Boltzmann equation to the multispecies problems. An asymptotic-preserving scheme for the kinetic equation is very efficient in the fluid regime where the Knudsen number is small and the collision term becomes stiff. It allows a coarse (independent of the Knudsen number) mesh size and a large time step in the fluid regime. The difficulty associated with multispecies problems is that there are no local conservation laws for each species, resulting in extra stiff nonlinear source terms that need to be discretized properly in order to (1) avoid Newton-type solvers for nonlinear algebraic systems, and (2) to be asymptotic-preserving. We show that these extra nonlinear source terms can be solved using only linear system solvers, and the scheme preserves the correct Euler and Navier-Stokes limits. Numerical examples are used to demonstrate the efficiency and applicability of the schemes for both Euler and Navier-Stokes regimes.

A higher-order Godunov method for radiation hydrodynamics: Radiation subsystem

A higher order Godunov method for the radiation subsystem of radiation hydrodynamics is presented. A key ingredient of the method is the direct coupling of stiff source term effects to the hyperbolic structure of the system of conservation laws; it is composed of a predictor step that is based on Duhamel's principle and a corrector step that is based on Picard iteration. The method is second order accurate in both time and space, unsplit, asymptotically preserving, and uniformly well behaved from the photon free streaming (hyperbolic) limit through the weak equilibrium diffusion (parabolic) limit and to the strong equilibrium diffusion (hyperbolic) limit. Numerical tests demonstrate second order convergence across various parameter regimes.

Dispersive nonlinear shallow-water equations: some preliminary numerical results

The dispersive nonlinear shallow-water equations of Antuono et al. (Stud Appl Math 122:1–28, 2008) are solved by means of an explicit arbitrary high-order accurate finite-volume scheme for nonlinear hyperbolic systems with stiff source terms. Tests against typical benchmark solutions are used to illustrate the robustness and accuracy of the solver while typical solutions for the propagation of solitary waves on a slope highlight the solution value in reproducing nearshore flows.

ADER schemes on unstructured meshes for nonconservative hyperbolic systems: Applications to geophysical flows

We develop a new family of well-balanced path-conservative quadrature-free one-step ADER finite volume and discontinuous Galerkin finite element schemes on unstructured meshes for the solution of hyperbolic partial differential equations with non-conservative products and stiff source terms. The fully discrete formulation is derived using the recently developed framework of explicit one-step P N P M schemes of arbitrary high order of accuracy in space and time for conservative hyperbolic systems [Dumbser M, Balsara D, Toro EF, Munz CD. A unified framework for the construction of one-step finite-volume and discontinuous Galerkin schemes. J Comput Phys 2008;227:8209–53]. The two key ingredients of our high order approach are: first, the high order accurate P N P M reconstruction operator on unstructured meshes, using the WENO strategy presented in [Dumbser M, Käser M, Titarev VA Toro EF. Quadrature-free non-oscillatory finite volume schemes on unstructured meshes for nonlinear hyperbolic systems. J Comput Phys 2007;226:204–43] to ensure monotonicity at discontinuities, and second, a local space–time Galerkin scheme to predict the evolution of the reconstructed polynomial data inside each element during one time step to obtain a high order accurate one-step time discretization. This approach is also able to deal with stiff source terms as shown in [Dumbser M, Enaux C, Toro EF. Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws. J Comput Phys 2008;227:3971–4001]. These two key ingredients are combined with the recently developed path-conservative methods of Parés [Parés C. Numerical methods for nonconservative hyperbolic systems: a theoretical framework. SIAM J Numer Anal 2006;44:300–21] and Castro et al. [Castro MJ, Gallardo JM, Parés C. High-order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. Applications to shallow-water systems. Math Comput 2006;75:1103–34] to treat the non-conservative products properly. We show applications of our method to the two-layer shallow water equations as well as applications to the recently published depth-averaged two-fluid flow model of Pitman and Le [Pitman EB, Le L. A two-fluid model for avalanche and debris flows. Philos Trans Roy Soc A 2005;363:1573–601].

An implicit box scheme for subsonic compressible flow with dissipative source term

We investigate the stability and convergence of an implicit box scheme for subsonic flows modelled by scalar conservation laws with dissipative and possibly stiff source terms. The scheme is proposed for solving transient gas flow problems in pipeline networks. Such networks are operated in the subsonic flow region and are characterized by pressure losses due to dissipative friction terms. We verify the properties stated by Kružkov’s theorem (Kružkov, Math. USSR-Sb. 10:217–243, 1970) for the approximate solution and prove its convergence to the entropy solution.

A second-order TVD implicit–explicit finite volume method for time-dependent convection-reaction equations

A class of conservative methods is developed in the more general framework of cell-centered upwind differences to approximate numerically the solution of one-dimensional non-linear conservation laws with (possibly) stiff reaction source terms. These methods are based on a non-oscillatory piecewise linear polynomial representation of the discrete solution within any mesh interval to compute pointwise solution values. The piecewise linear approximate solution is obtained by approximating the cell average of the analytical solution and the solution slope in every mesh cell. These two quantities are evolved in time by solving a set of discrete equations that are suitably designed to ensure formal second-order consistency. Several numerical tests which are taken from literature illustrate the performance of the method in solving non-stiff and stiff convection-reaction equations in conservative form.

Conditional Moment Closure (CMC) applied to autoignition of high pressure methane jets in a shock tube

Autoignition of non-premixed methane–air mixtures is investigated using first-order Conditional Moment Closure (CMC). Turbulent velocity and mixing fields simulations are decoupled from the CMC calculations due to low temperature changes until ignition occurs. The CMC equations are cross-stream averaged and finite differences are applied to discretize the equations. A three-step fractional method is implemented to treat separately the stiff chemical source term. Two detailed chemical kinetics mechanisms are tested as well as two mixing models. The present results show good agreement with published experimental measurements for the magnitude of both ignition delay and kernel location. The slope of the predicted ignition delay is overpredicted and possible sources of discrepancy are identified. Both scalar dissipation rate models produce comparable results due to the turbulent flow homogeneity assumption. Further, ignition always occurs at low scalar dissipation rates, much lower than the flamelet critical value of ignition. Ignition is found to take place in lean mixtures for a value of mixture fraction around 0.02. The conditional species concentrations are in qualitative agreement with previous research. Homogeneous and inhomogeneous CMC calculations are also performed in order to investigate the role of physical transport in the present autoignition study. It is found that spatial transport is small at ignition time. Predicted ignition delays are shown to be sensitive to the chemical kinetics. Reasonable agreement with previous simulations is found. Improved formulations for the mixing model based on non-homogeneous turbulence are expected to have an impact.

Modelling of autoignition for methane-based fuel blends using Conditional Moment Closure

First-order Conditional Moment Closure (CMC) is applied to investigate the autoignition of high pressure CH4, CH4/C2H6, CH4/H2, and CH4/N2 jets injected into lower pressure heated air over a wide range of pre-combustion temperatures. Detailed and optimized chemical kinetics mechanisms are used for each fuel blend. The turbulent flow field and CMC calculations are decoupled based on the frozen mixing assumption. The standard k–ε turbulence model is used to perform the velocity and mixing fields calculations. The CMC equations are cross-stream averaged due to the small radial dependence of the conditional scalars in shear flows, and are solved elliptically. A three-step fractional method is implemented to treat the stiff chemical source term separately. Two formulations are implemented to model the conditional scalar dissipation term. The numerical results are compared to shock tube experiments. Predicted ignition delays and locations are found to be in reasonable agreement with the experiments. Further analysis is provided.

Construction of a multirate RODAS method for stiff ODEs

Multirate time stepping is a numerical technique for efficiently solving large-scale ordinary differential equations (ODEs) with widely different time scales localized over the components. This technique enables one to use large time steps for slowly varying components, and small steps for rapidly varying ones. Multirate methods found in the literature are normally of low order, one or two. Focusing on stiff ODEs, in this paper we discuss the construction of a multirate method based on the fourth-order RODAS method. Special attention is paid to the treatment of the refinement interfaces with regard to the choice of the interpolant and the occurrence of order reduction. For stiff, linear systems containing a stiff source term, we propose modifications for the treatment of the source term which overcome order reduction originating from such terms and which we can implement in our multirate method.

Godunov’s Method for Simulatinons of Fluid-Structure Interaction in Piping Systems

Constant coefficient one-dimensional linear hyperbolic systems of partial differential equations (PDEs) are often used for description of fluid-structure interaction (FSI) phenomena during transient conditions in piping systems. In the past, these systems of equations have been numerically solved with method of characteristics (MOCs). The MOC method is actually the most efficient and accurate method for description of the single-phase transient in the cold liquid where the constant coefficient mathematical model describes phenomenon with sufficient accuracy. In energy production systems where hot pressurized liquid is used for heat transfer between the heat source and the steam generator, more complex and nonlinear mathematical models are needed to describe transient flow and these models cannot be solved with MOC method because the models are not constant. In addition, the MOC method can be used for pipes having discontinuities like elbows, geometrical changes, material properties changes, etc., but only with some extra numerical modeling. An interesting alternative is explicit characteristic upwind numerical method, known as Godunov’s method that is frequently used for nonlinear systems or systems where properties change with position. In the present study, applicability of the Godunov’s method for the FSI analyses is tested with eight first order PDEs mathematical model. The conventional linear mathematical model is improved with convective term that makes the system nonlinear and additional terms that enable simulations of the FSI in arbitrarily shaped piping systems located in a plane. Two PDEs describe pressure waves in the single-phase fluid and six PDEs describe axial, lateral, and rotational stress waves in the pipe. The applied system of equations has stiff source terms. This numerical problem is solved introducing implicit iterations. The proposed model is verified with a rod impact experiment that is carried out on single-elbow pipe hanging on wires. Godunov’s method is found as a very promising numerical method for simulations of the FSI problems.

A computational study of inviscid hypersonic flows using energy relaxation method

Reasonable analysis of hypersonic flows requires a thermodynamic non-equilibrium model to properly simulate strong shock waves or high pressure and temperature states in the flow field. The energy relaxation method (ERM) has been used to model such a non-equilibrium effect which is generally expressed as a hyperbolic system of equations with a stiff relaxation source term. Relaxation time that is multiplied with source terms is responsible for non-equilibrium in the system. In the present study, a numerical analysis has been carried out with varying values of relaxation time for several hypersonic flows with AUSM (advection upstream splitting method) as a numerical scheme. Vibration modes of thermodynamic nonequilibrium effects are considered. The results obtained showed that, as the relaxation time reduces to zero, the solution marches toward equilibrium, while it shows non-equilibrium effects, as the relaxation time increases. The present computations predicted the experiment results of hypersonic flows with good accuracy. The work carried out suggests that the present energy relaxation method can be robust for analysis of hypersonic flows.

Numerical investigation of ethylene flame bubble instability induced by shock waves

In this paper, the ethylene/oxygen/nitrogen premixed flame instabilities induced by incident and reflected shock wave were investigated numerically. The effects of grid resolutions and chemical mechanisms on the flame bubble deformation process are evaluated. In the computational frame, the 2D multi-component Navier–Stokes equations with second-order flux-difference splitting scheme were used; the stiff chemical source term was integrated using an implicit ordinary differential equations (ODEs) solver. The two ethylene/oxygen/nitrogen chemical mechanisms, namely 3-step reduced mechanism and 35-step elementary skeletal mechanism, were used to examine the reliability of chemistry. On the other hand, the different grid sizes, Δx × Δy = 0.25 × 0.5mm and Δx × Δy = 0.15 × 0.2mm, were implemented to examine the accuracy of the grid resolution. The computational results were qualitatively validated with experimental results of Thomas et al. (Combust Theory Model 5:573–594, 2001). Two chemical mechanisms and two grid resolutions used in present study can qualitatively reproduce the ethylene spherical flame instability process generated by an incident shock wave of Mach number 1.7. For the case of interaction between the flame and reflected shock waves, the 35-steps mechanism qualitatively predicts the physical process and is somewhat independent on the grid resolutions, while the 3-steps mechanism fails to reproduce the instability of ethylene flame for the two selected grid resolutions. It is concluded that the detailed chemical mechanism, which includes the chain elementary reactions of fuel combustion, describes the flame instability induced by shock wave, in spite of the fact that the flame thickness (reaction zone) is represented by 1–2 grids only.

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.

Numerical solution of under-resolved detonations

A new fractional-step method is proposed for the numerical solution of high speed reacting flows, where the chemical time scales are often much smaller than the fluid dynamical time scales. When the problem is stiff, because of insufficient spatial/temporal resolution, a well-known spurious numerical phenomenon occurs in standard finite volume schemes: the incorrect calculation of the speed of propagation of discontinuities. The new method is first illustrated considering a one-dimensional scalar hyperbolic advection/reaction equation with stiff source term, which may be considered as a model problem to under-resolved detonations. During the reaction step, the proposed scheme replaces the cell average representation with a two-value reconstruction, which allows us to locate the discontinuity position inside the cell during the computation of the source term. This results in the correct propagation of discontinuities even in the stiff case. The method is proved to be second-order accurate for smooth solutions of scalar equations and is applied successfully to the solution of the one-dimensional reactive Euler equations for Chapman–Jouguet detonations.