Stiff Source Terms Research Articles

Exponential Runge-Kutta Methods for the Multispecies Boltzmann Equation

AbstractThis paper generalizes the exponential Runge-Kutta asymptotic preserving (AP) method developed in [G. Dimarco and L. Pareschi,SIAM Numer. Anal., 49 (2011), pp. 2057-2077] to compute the multi-species Boltzmann equation. Compared to the single species Boltzmann equation that the method was originally applied on, this set of equation presents a new difficulty that comes from the lack of local conservation laws due to the interaction between different species. Hence extra stiff nonlinear source terms need to be treated properly to maintain the accuracy and the AP property. The method we propose does not contain any nonlinear nonlocal implicit solver, and can capture the hydrodynamic limit with time step and mesh size independent of the Knudsen number. We prove the positivity and strong AP properties of the scheme, which are verified by two numerical examples.

Computational challenges for simulations related to the NASA electric arc shock tube (EAST) experiments

The goal of this study is to gain some physical insights and an understanding of the computational challenges for the simulations related to the hypersonic nonequilibrium multi-species and multi-reaction experiments on the NASA Electric Arc Shock Tube (EAST). While experimental measurement does not provide any information about the radial structure of this type of flow, accurate and reliable numerical simulations can provide more insight into the physical structure of the flow to aid the design of atmospheric entry spacecrafts. The paper focuses on the spurious numerics which take place in numerical simulations of the subject physics containing stiff source terms and discontinuities. This paper is based on the knowledge gained from Yee et al. on simple reacting test cases (Yee et al. 2013, [9]) as a guide to reveal the computational challenges involved for such an extreme flow type. The results of the 1D and 2D EAST viscous and inviscid simulations using a simplified physical model are presented. The computation reveals, for the first time, that the 2D viscous model which contains both shocks and shears exhibits Tollmien–Schlichting-like instability complex patterns at the boundary layer. In addition to exhibiting spurious numerical behavior of wrong propagation speed of discontinuities by typical shock-capturing methods, there is improved understanding on the cause of numerical difficulties by previous investigators. One example is that the relative distance between the shocks and shear/contact is different from one grid spacing to another for each considered high order shock-capturing scheme. The results presented can provide insight on the numerical instability observed by previous investigations and future algorithm development for this type of extreme flow.

Hyperbolic reformulation of a 1D viscoelastic blood flow model and ADER finite volume schemes

The applicability of ADER finite volume methods to solve hyperbolic balance laws with stiff source terms in the context of well-balanced and non-conservative schemes is extended to solve a one-dimensional blood flow model for viscoelastic vessels, reformulated as a hyperbolic system, via a relaxation time. A criterion for selecting relaxation times is found and an empirical convergence rate assessment is carried out to support this result. The proposed methodology is validated by applying it to a network of viscoelastic vessels for which experimental and numerical results are available. The agreement between the results obtained in the present paper and those available in the literature is satisfactory. Key features of the present formulation and numerical methodologies, such as accuracy, efficiency and robustness, are fully discussed in the paper.

The equilibrium state method for hyperbolic conservation laws with stiff reaction terms

A new fractional-step method is proposed for numerical simulations of hyperbolic conservation laws with stiff source terms arising from chemically reactive flows. In stiff reaction problems, a well-known spurious numerical phenomenon, the incorrect propagation speed of discontinuities, may occur in general fractional-step algorithm due to the underresolved numerical solution in both space and time. The basic idea of the present proposed scheme is to replace the cell average representation with a two-equilibrium states reconstruction during the reaction step, which allows us to obtain the correct propagation of discontinuities for stiff reaction problems in an underresolved mesh. Because the definition of these two-equilibrium states for each transition cell is independent of its neighboring cells, the proposed method can be extended to multi-dimensional problems directly. In addition, this method is promising to deal with more complicated real-world problems after being extended to multi-species/multi-reactions system. Extensive numerical examples for one- and two-dimensional scalar and Euler system demonstrate the reliability and robustness of this novel method.

A strong conservative Riemann solver for the solution of the coupled Maxwell and Navier–Stokes equations

The coupled system of the Navier–Stokes and Maxwell equations are recast into a strong conservative form, which allows the fluid coupling to the Maxwell system to be written in terms of flux divergence rather than explicit source terms. This effectively removes source terms from the Navier–Stokes equations, although retaining an exact coupling to the electromagnetics. While this relieves the stiff source terms and potentially stabilizes the system, it introduces a much more complicated eigenstructure to the governing equations. The flux Jacobian and eigenvectors for this strong conservative system are presented in the current paper for the first time. An approximate Riemann solver based upon these eigenvectors is then introduced and tested. The solver is implemented in a preconditioned, dual-time implicit form. Validations for classic one- and two-dimensional problems are presented, and the performances of the new formulation and the traditional source-coupled formulation are compared.

A TWO-MOMENT RADIATION HYDRODYNAMICS MODULE IN ATHENA USING A TIME-EXPLICIT GODUNOV METHOD

We describe a module for the Athena code that solves the gray equations of radiation hydrodynamics (RHD), based on the first two moments of the radiative transfer equation. We use a combination of explicit Godunov methods to advance the gas and radiation variables including the non-stiff source terms, and a local implicit method to integrate the stiff source terms. We adopt the M1 closure relation and include all leading source terms. We employ the reduced speed of light approximation (RSLA) with subcycling of the radiation variables in order to reduce computational costs. Our code is dimensionally unsplit in one, two, and three space dimensions and is parallelized using MPI. The streaming and diffusion limits are well-described by the M1 closure model, and our implementation shows excellent behavior for a problem with a concentrated radiation source containing both regimes simultaneously. Our operator-split method is ideally suited for problems with a slowly varying radiation field and dynamical gas flows, in which the effect of the RSLA is minimal. We present an analysis of the dispersion relation of RHD linear waves highlighting the conditions of applicability for the RSLA. To demonstrate the accuracy of our method, we utilize a suite of radiation and RHD tests covering a broad range of regimes, including RHD waves, shocks, and equilibria, which show second-order convergence in most cases. As an application, we investigate radiation-driven ejection of a dusty, optically thick shell in the interstellar medium (ISM). Finally, we compare the timing of our method with other well-known iterative schemes for the RHD equations. Our code implementation, Hyperion, is suitable for a wide variety of astrophysical applications and will be made freely available on the Web.

Additive Runge-Kutta methods for H2/O2/Ar detonation with a detailed elementary chemical reaction model

We report here the additive Runge-Kutta methods for computing reactive Euler equations with a stiff source term, and in particular, their applications in gaseous detonation simulations. The source term in gaseous detonation is stiff due to the presence of wide range of time scales during thermal-chemical non-equilibrium reactive processes and some of these time scales are much smaller than that of hydrodynamic flow. The high order, L-stable, additive Runge-Kutta methods proposed in this paper resolved the stiff source term into the stiff part and non-stiff part, in which the stiff part was solved implicitly while the non-stiff part was handled explicitly. The proposed method was successfully applied to simulating the gaseous detonation in a stoichiometric H2/O2/Ar mixture based on a detailed elementary chemical reaction model comprised of 9 species and 19 elementary reactions. The results showed that the stiffly accurate additive Runge-Kutta methods can capture the discontinuity well, and describe the detonation complex wave configurations accurately such as the triple wave structure and cellular pattern.

Spurious behavior of shock-capturing methods by the fractional step approach: Problems containing stiff source terms and discontinuities

The goal of this paper is to relate numerical dissipations that are inherited in high order shock-capturing schemes with the onset of wrong propagation speed of discontinuities. For pointwise evaluation of the source term, previous studies indicated that the phenomenon of wrong propagation speed of discontinuities is connected with the smearing of the discontinuity caused by the discretization of the advection term. The present study focuses only on solving the reactive system by the fractional step method using the Strang splitting. Studies shows that the degree of wrong propagation speed of discontinuities is highly dependent on the accuracy of the numerical method. The manner in which the smearing of discontinuities is contained by the numerical method and the overall amount of numerical dissipation being employed play major roles. Depending on the numerical method, time step and grid spacing, the numerical simulation may lead to (a) the correct solution (within the truncation error of the scheme), (b) a divergent solution, (c) a wrong propagation speed of discontinuities solution or (d) other spurious solutions that are solutions of the discretized counterparts but are not solutions of the governing equations. The findings might shed some light on the reported difficulties in numerical combustion and problems with stiff nonlinear (homogeneous) source terms and discontinuities in general.

Prospects of UDMH Organic Gel Droplet Evaporation and Combustion Model

The organic gel propellants have the advantages of both solid and liquid propellant, but the catching fire and burning process of organic gel droplet appears obvious vapor ejection and flame disturbance. This unsteady combustion phe- nomenon of organic gel droplet is different from the quasi-steady combustion of conventional liquid droplet and produces great influence on the spray combustion flow field of gel liquid rocket engine and its optimization design. Based on the un- steady evaporation characteristics inside UDMH (unsymmetrical dimethylhydrazine) organic gel droplet and the micro molecule structure of organic gel propellant, three kinds of liquid phase evaporation model, the unsteady multicomponent evaporation model, the interface tracking multiphase flow evaporation model and the discrete element multiphase flow evaporation model, are proposed. The unsteady multicomponent evaporation model regards the gel droplet as a binary multi- conponent droplet, in which the gellant inside the gel droplet is been equivalent to a species with low volatility and high boiling temperature. In order to depict the shape and location variety of vapor bubble inside the gel droplet more precisely, the interface tracking multiphase flow evaporation model regards the gel droplet and the vapor bubble inside it as a single fluid, and uses the Level Set function tracking the interface location of vapor bubble. Based on the micro molecule structure of organic gel propellant, the discrete element multiphase flow evaporation model for the gel droplet regards the gellant and vapor bubble as the aggregation of micro-particles with different properties respectively, and uses the Lagrange method tracking the movement and aggregation of each kind micro-particle. The multicomponent evaporation model has been im- plemented in our foregoing work, so a summary of the model is presented. And next the main difficulties and the technique preparation about the last two multiphase flow models have been presented and analyzed in detail. In view of mixing and combustion process of gas phase surrounding the UDMH organic gel droplet, we propose that the finite rate chemical reac- tion model and the flamelet model should be used for the stagnation ambient and convective ambient respectively. The diffi- culties and numerical methods for handling the stiff chemical source term in finite rate chemical reaction model have been presented. And next the consistent flamelet model which takes consideration of species diffusion coefficient and non-unity Lewis number is presented. It would be more precise for the temperature and species interpolation of the gel droplet laminar non-premixed combustion flow field. Keywords organic gel; droplet; evaporation model; finite rate chemistry model; flamelet model; stiff source term

High-Resolution Finite Volume Modeling of Wave Propagation in Orthotropic Poroelastic Media

Poroelasticity theory models the dynamics of porous, fluid-saturated media. It was pioneered by Maurice Biot in the 1930s through the 1960s and has applications in several fields, including geophysics and modeling of in vivo bone. A wide variety of methods have been used to model poroelasticity, including finite difference, finite element, pseudospectral, and discontinuous Galerkin methods. In this work we use a Cartesian-grid high-resolution finite volume method to numerically solve Biot's equations in the time domain for orthotropic materials with the stiff relaxation source term in the equations incorporated using operator splitting. This class of finite volume method has several useful properties, including the ability to use wave limiters to reduce numerical artifacts in the solution, ease of incorporating material inhomogeneities, low memory overhead, and an explicit time-stepping approach. To the authors' knowledge, this is the first use of high-resolution finite volume methods to model poroelastic...

Asymptotic‐preserving Godunov‐type numerical schemes for hyperbolic systems with stiff and nonstiff relaxation terms

AbstractWe devise a new class of asymptotic‐preserving Godunov‐type numerical schemes for hyperbolic systems with stiff and nonstiff relaxation source terms governed by a relaxation time ε. As an alternative to classical operator‐splitting techniques, the objectives of these schemes are twofold: first, to give accurate numerical solutions for large, small, and in‐between values of ε and second, to make optional the choice of the numerical scheme in the asymptotic regime ε tends to zero. The latter property may be of particular interest to make easier and more efficient the coupling at a fixed spatial interface of two models involving very different values of ε. © 2012 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2013

GPU-based flow simulation with detailed chemical kinetics

The current paper reports on the implementation of a numerical solver on the Graphic Processing Units (GPUs) to model reactive gas mixtures with detailed chemical kinetics. The solver incorporates high-order finite volume methods for solving the fluid dynamical equations coupled with stiff source terms. The chemical kinetics are solved implicitly via an operator-splitting method. We explored different approaches in implementing a fast kinetics solver on the GPU. The detail of the implementation is discussed in the paper. The solver is tested with two high-order shock capturing schemes: MP5 (Suresh and Huynh, 1997) [9] and ADERWENO (Titarev and Toro, 2005) [10]. Considering only the fluid dynamics calculation, the speed-up factors obtained are 30 for the MP5 scheme and 55 for ADERWENO scheme. For the fully-coupled solver, the performance gain depended on the size of the reaction mechanism. Two different examples of chemistry were explored. The first mechanism consisted of 9 species and 38 reactions, resulting in a speed-up factor up to 35. The second, larger mechanism, consisted of 36 species and 308 reactions, resulting in a speed-up factor of up to 40.

General relativistic radiation hydrodynamics of accretion flows - II. Treating stiff source terms and exploring physical limitations

We present the implementation of an implicit-explicit (IMEX) Runge-Kutta numerical scheme for general relativistic hydrodynamics coupled to an optically thick radiation field in two existing GR-hydrodynamics codes. We argue that the necessity of such an improvement arises naturally in astrophysically relevant regimes where the optical thickness is high as the equations become stiff. By performing several 1D tests we verify the codes' new ability to deal with this stiffness and show consistency. Then, still in 1D, we compute a luminosity versus accretion rate diagram for the setup of spherical accretion onto a Schwarzschild black hole and find good agreement with previous work. Lastly, we revisit the supersonic Bondi Hoyle Lyttleton (BHL) accretion in 2D where we can now present simulations of realistic temperatures, down to T~10^6 K. Here we find that radiation pressure plays an important role, but also that these highly dynamical set-ups push our approximate treatment towards the limit of physical applicability. The main features of radiation hydrodynamics BHL flows manifest as (i) an effective adiabatic index approaching gamma_effective ~ 4/3; (ii) accretion rates two orders of magnitude lower than without radiation pressure; (iii) luminosity estimates around the Eddington limit, hence with an overall radiative efficiency as small as eta ~ 10^{-2}; (iv) strong departures from thermal equilibrium in shocked regions; (v) no appearance of the flip-flop instability. We conclude that the current optically thick approximation to the radiation transfer does give physically substantial improvements over the pure hydro also in set-ups departing from equilibrium, and, once accompanied by an optically thin treatment, is likely to provide a fundamental tool for investigating accretion flows in a large variety of astrophysical systems.

BSSN equations in spherical coordinates without regularization: Vacuum and nonvacuum spherically symmetric spacetimes

Brown [Phys. Rev. D 79, 104029 (2009)] has recently introduced a covariant formulation of the BSSN equations which is well suited for curvilinear coordinate systems. This is particularly desirable as many astrophysical phenomena are symmetric with respect to the rotation axis or are such that curvilinear coordinates adapt better to their geometry. However, the singularities associated with such coordinate systems are known to lead to numerical instabilities unless special care is taken (e.g., regularization at the origin). Cordero-Carri\'on will present a rigorous derivation of partially implicit Runge-Kutta methods in forthcoming papers, with the aim of treating numerically the stiff source terms in wavelike equations that may appear as a result of the choice of the coordinate system. We have developed a numerical code solving the BSSN equations in spherical symmetry and the general relativistic hydrodynamic equations written in flux-conservative form. A key feature of the code is that it uses a second-order partially implicit Runge-Kutta method to integrate the evolution equations. We perform and discuss a number of tests to assess the accuracy and expected convergence of the code---namely a pure gauge wave, the evolution of a single black hole, the evolution of a spherical relativistic star in equilibrium, and the gravitational collapse of a spherical relativistic star leading to the formation of a black hole. We obtain stable evolutions of regular spacetimes without the need for any regularization algorithm at the origin.

On stability of difference schemes. Central schemes for hyperbolic conservation laws with source terms

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.

Discontinuous Galerkin Approximation of Relaxation Models for Linear and Nonlinear Diffusion Equations

In this work we present finite element approximations of relaxed systems for nonlinear diffusion problems, which can also tackle the cases of degenerate and strongly degenerate diffusion equations. Relaxation schemes take advantage of the replacement of the original partial differential equation (PDE) with a semilinear hyperbolic system of equations, with a stiff source term, tuned by a relaxation parameter $\epsilon$. When $\epsilon \rightarrow 0^+$, the system relaxes onto the original PDE: in this way, a consistent discretization of the relaxation system for vanishing $\epsilon$ yields a consistent discretization of the original PDE. The numerical schemes obtained with this procedure do not require solving implicit nonlinear problems and possess the robustness of upwind discretizations. The proposed approximations are based on a discontinuous Galerkin method in space and on suitable implicit-explicit integration in time. Then, in principle, we can achieve any order of accuracy and obtain stable solutions, even when the diffusion equation becomes degenerate and solution singularities develop. Moreover, when needed, we can easily incorporate slope limiters within our schemes in order to handle spurious oscillatory phenomena. Some preliminary theoretical results are given, along with several numerical tests in one and two space dimensions, both for linear and nonlinear diffusion problems, including a degenerate diffusion equation, that provide numerical evidence of the properties of the presented approach.

THE PLUTO CODE FOR ADAPTIVE MESH COMPUTATIONS IN ASTROPHYSICAL FLUID DYNAMICS

We present a description of the adaptive mesh refinement (AMR) implementation of the PLUTO code for solving the equations of classical and special relativistic magnetohydrodynamics (MHD and RMHD). The current release exploits, in addition to the static grid version of the code, the distributed infrastructure of the CHOMBO library for multidimensional parallel computations over block-structured, adaptively refined grids. We employ a conservative finite-volume approach where primary flow quantities are discretized at the cell-center in a dimensionally unsplit fashion using the Corner Transport Upwind (CTU) method. Time stepping relies on a characteristic tracing step where piecewise parabolic method (PPM), weighted essentially non-oscillatory (WENO) or slope-limited linear interpolation schemes can be handily adopted. A characteristic decomposition-free version of the scheme is also illustrated. The solenoidal condition of the magnetic field is enforced by augmenting the equations with a generalized Lagrange multiplier (GLM) providing propagation and damping of divergence errors through a mixed hyperbolic/parabolic explicit cleaning step. Among the novel features, we describe an extension of the scheme to include non-ideal dissipative processes such as viscosity, resistivity and anisotropic thermal conduction without operator splitting. Finally, we illustrate an efficient treatment of point-local, potentially stiff source terms over hierarchical nested grids by taking advantage of the adaptivity in time. Several multidimensional benchmarks and applications to problems of astrophysical relevance assess the potentiality of the AMR version of PLUTO in resolving flow features separated by large spatial and temporal disparities.

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.

Relaxation method for unsteady convection–diffusion equations

We propose and implement a relaxation method for solving unsteady linear and nonlinear convection–diffusion equations with continuous or discontinuity-like initial conditions. The method transforms a convection–diffusion equation into a relaxation system, which contains a stiff source term. The resulting relaxation system is then solved by a third-order accurate implicit–explicit (IMEX) Runge–Kutta method in time and a fifth-order finite difference WENO scheme in space. Numerical results show that the method can be used to effectively solve convection–diffusion equations with both smooth structures and discontinuities.

Relaxation model for extended magnetohydrodynamics: Comparison to magnetohydrodynamics for dense Z-pinches

It is shown that the two-fluid model under a generalized Ohm’s law formulation and the resistive magnetohydrodynamics (MHD) can both be described as relaxation systems. In the relaxation model, the under-resolved stiff source terms constrain the dynamics of a set of hyperbolic equations to give the correct asymptotic solution. When applied to the collisional two-fluid model, the relaxation of fast time scales associated with displacement current and finite electron mass allows for a natural transition from a system where Ohm’s law determines the current density to a system where Ohm’s law determines the electric field. This result is used to derive novel algorithms, which allow for multiscale simulation of low and high frequency extended-MHD physics. This relaxation formulation offers an efficient way to implicitly advance the Hall term and naturally simulate a plasma-vacuum interface without invoking phenomenological models. The relaxation model is implemented as an extended-MHD code, which is used to analyze pulsed power loads such as wire arrays and ablating foils. Two-dimensional simulations of pulsed power loads are compared for extended-MHD and MHD. For these simulations, it is also shown that the relaxation model properly recovers the resistive-MHD limit.