Time Step Restriction Research Articles

Conservative averaging-reconstruction techniques (Ring Average) for 3-D finite-volume MHD solvers with axis singularity

In cylindrical/spherical geometries, MHD solvers using explicit finite-volume methods face restrictive time steps imposed by the CFL condition due to the clustering of cells in the azimuthal direction near the pole/axis. We use a conservative averaging-reconstruction method (Ring Average) on structured cylindrical/spherical mesh to remove this severe time step restriction for multi-dimensional curvilinear finite-volume MHD solvers. The Ring Average technique is implemented as a post-processing step and thus requires no changes to the existing data structure, grid definition or numerical methods in the original MHD solver. This paper describes the Ring Average algorithm and presents simulation results using field loop advection and MHD blast waves in cylindrical/spherical geometries for validation. The algorithm is shown to be inexpensive and easily implemented in existing curvilinear finite-volume MHD solvers.

On the optimal CFL number of SSP methods for hyperbolic problems

We show that the theory for strong stability preserving (SSP) time stepping methods employed with the method of lines-type discretizations of hyperbolic conservation laws may result in overly stringent time step restrictions. We analyze a fully discrete finite volume method with slope reconstruction and a second order SSP Runge–Kutta time integrator to show that the maximum stable time step can be increased over the SSP limit. Numerical examples indicate that this result extends to two-dimensional problems on triangular meshes.

Hyperbolic phase field modeling of brittle fracture: Part I—Theory and simulations

We propose a new phase field model of dynamic brittle fracture, in which material damage evolves according to a hyperbolic partial differential equation. This model can be stably discretized using explicit time integration, without imposing crippling time step restrictions with refinement in space. The model is derived from microforce balance by including effects of microscopic inertia. Quantitative predictions of the proposed model differ from those of parabolic and elliptic phase field models in that there is a mild rate-toughening effect, but major qualitative solution features are essentially the same. We compute finite element approximations of solutions to several dynamic fracture scenarios to support these claims. Part II of this series will incorporate the proposed model into a hybrid isogeometric–meshfree framework for air-blast–structure interaction and provide further demonstrations of the model’s physical and numerical properties.

An in-cell reconstruction finite volume method for flows of compressible immiscible fluids

We present a hybrid capturing/tracking method for finite volume solvers of compressible flows involving two immiscible fluids that are described by equations of states for ideal gases. The approach is an extension of the two-dimensional level set based in-cell-reconstruction method originally proposed by Smiljanovski et al. [1] for deflagration waves to three-dimensional flows involving interfaces between two immiscible fluids with surface tension forces. The interface motion is captured by an extension to the conservative, un-split geometric volume-of-fluid technique of Owkes and Desjardins [2]. The resulting method is a sharp interface method that avoids time-step restrictions or merging/mixing rules due to cut-cells by using cell-face-aperture averaged waves to update the volume averaged states of cells containing the interface directly. The method furthermore avoids the use of any mixed states in discretization stencils, by reconstructing the pure fluid states in each mixed cell, explicitly enforcing the jump conditions across the interface. Simulations of compressible flows that involve the interaction of shocks with interfaces between immiscible fluids are performed to demonstrate the performance of the proposed method.

Stability analysis of MacCormack rapid solver method for evolutionary Stokes–Darcy problem

A numerical method based on the explicit–implicit scheme of MacCormack finite difference scheme was applied to the solution of Parabolized Navier–Stokes equations. Partitioned methods, which solve the coupled problem by successively solving the sub-physics problems, have recently been studied for the full evolutionary groundwater–surface water flow with convergence established over both bounded and long time intervals. A MacCormack rapid solver method based on interface approximation via temporal extrapolation is proposed for devising decoupled marching algorithms for the mixed model. Under a time step limitation of the form C1h≤Δt≤C2h (where h, Δt are mesh size and time step, respectively, and Cj, j∈{1,2}, are two physical parameters) we prove both uniform and asymptotic stabilities over long time intervals. Some numerical experiments are presented and discussed.

Numerical Preservation of Velocity Induced Invariant Regions for Reaction\u2013Diffusion Systems on Evolving Surfaces

We propose and analyse a finite element method with mass lumping (LESFEM) for the numerical approximation of reaction–diffusion systems (RDSs) on surfaces in {mathbb {R}}^3 that evolve under a given velocity field. A fully-discrete method based on the implicit–explicit (IMEX) Euler time-discretisation is formulated and dilation rates which act as indicators of the surface evolution are introduced. Under the assumption that the mesh preserves the Delaunay regularity under evolution, we prove a sufficient condition, that depends on the dilation rates, for the existence of invariant regions (i) at the spatially discrete level with no restriction on the mesh size and (ii) at the fully-discrete level under a timestep restriction that depends on the kinetics, only. In the specific case of the linear heat equation, we prove a semi- and a fully-discrete maximum principle. For the well-known activator-depleted and Thomas reaction–diffusion models we prove the existence of a family of rectangles in the phase space that are invariant only under specific growth laws. Two numerical examples are provided to computationally demonstrate (i) the discrete maximum principle and optimal convergence for the heat equation on a linearly growing sphere and (ii) the existence of an invariant region for the LESFEM–IMEX Euler discretisation of a RDS on a logistically growing surface.

Sharp interface approaches and deep learning techniques for multiphase flows

We present a review on numerical methods for simulating multiphase and free surface flows. We focus in particular on numerical methods that seek to preserve the discontinuous nature of the solutions across the interface between phases. We provide a discussion on the Ghost-Fluid and Voronoi Interface methods, on the treatment of surface tension forces that avoid stringent time step restrictions, on adaptive grid refinement techniques for improved efficiency and on parallel computing approaches. We present the results of some simulations obtained with these treatments in two and three spatial dimensions. We also provide a discussion of Machine Learning and Deep Learning techniques in the context of multiphase flows and propose several future potential research thrusts for using deep learning to enhance the study and simulation of multiphase flows.

Application of implicit Roe-type scheme and Jacobian-Free Newton-Krylov method to two-phase flow problems

A new implicit second-order accurate spatial scheme for steady-state thermal-hydraulic simulations of the two-phase two-fluid six-equation model is proposed. The new scheme is based on a Roe-type numerical flux that is formulated with the help of analytical approximate eigenvalues and eigenvectors of the two-phase system. Approximate eigenvalues and eigenvectors are obtained with a structured Jacobian matrix that is general for arbitrary Equation of State (EOS). Three issues are solved in this article: (1) A fully-implicit scheme using the backward Euler method is proposed to improve the stability and avoid the time step limit of the originally explicit scheme. The fully implicit scheme is solved with the Jacobian-Free Newton-Krylov method. (2) A second-order accurate spatial scheme is proposed by converting the first-order Roe-type numerical flux to a second-order one. The conversion is made by extending an existing procedure for single-phase flows to two-phase flows. (3) Phase appearance issue is treated with a numerical procedure, which requires little modification to the schemes. The new scheme is verified and validated with two numerical tests: the faucet flow problem and the BFBT benchmark.

A Particle Module for the PLUTO Code. I. An Implementation of the MHD–PIC Equations

Abstract We describe an implementation of a particle physics module available for the PLUTO code appropriate for the dynamical evolution of a plasma consisting of a thermal fluid and a nonthermal component represented by relativistic charged particles or cosmic rays (CRs). While the fluid is approached using standard numerical schemes for magnetohydrodynamics, CR particles are treated kinetically using conventional Particle-In-Cell (PIC) techniques. The module can be used either to describe test-particle motion in the fluid electromagnetic field or to solve the fully coupled magnetohydrodynamics (MHD)–PIC system of equations with particle backreaction on the fluid as originally introduced by Bai et al. Particle backreaction on the fluid is included in the form of momentum–energy feedback and by introducing the CR-induced Hall term in Ohm’s law. The hybrid MHD–PIC module can be employed to study CR kinetic effects on scales larger than the (ion) skin depth provided that the Larmor gyration scale is properly resolved. When applicable, this formulation avoids resolving microscopic scales, offering substantial computational savings with respect to PIC simulations. We present a fully conservative formulation that is second-order accurate in time and space, and extends to either the Runge–Kutta (RK) or the corner transport upwind time-stepping schemes (for the fluid), while a standard Boris integrator is employed for the particles. For highly energetic relativistic CRs and in order to overcome the time-step restriction, a novel subcycling strategy that retains second-order accuracy in time is presented. Numerical benchmarks and applications including Bell instability, diffusive shock acceleration, and test-particle acceleration in reconnecting layers are discussed.

Dual-time schemes for reacting flows with stiff kinetics

Two operator-split schemes designed for handling stiff chemistry are proposed within a dual-time framework and tested using zero- and one-dimensional test problems. Both are shown to perform robustly without ad hoc time-step restrictions. The first scheme is a Physical Splitting Scheme that is a straightforward extension of the Strang splitting method. Like the Strang splitting method, it retains the issue of splitting errors degrading solution accuracy at large time steps. However, it performs well at small time-steps and remains competitive because the scheme involves only two ODE solutions per physical time step. The second scheme is the Consistent Splitting Scheme that is designed to eliminate the splitting errors, but is relatively more expensive because it requires two ODE solutions per pseudo-time step. This scheme also performs well at small to moderate time steps, but its accuracy also degrades at large time steps in part because of issues associated with finite precision. This study indicates the relative strengths and weaknesses of operator splitting approaches and points the way for continued research.

Long time stability and convergence rate of MacCormack rapid solver method for nonstationary Stokes–Darcy problem

We propose and study a combination of two second-order implicit–explicit (IMEX) methods for the coupled Stokes–Darcy system that governs flows in karst aquifers. The first is a second-order explicit two-step MacCormack scheme and the second is a second-order implicit Crank–Nicolson method. Both algorithms only require the solution of two decoupled problems at each time step, one Stokes and the other Darcy. This combination so called the MacCormack rapid solver method is very efficient (faster, at least of second order accuracy in time and space) and can be easily implemented using legacy codes. Under time step limitation of the form Δt≤Ch (where h,Δt are mesh size and time step, respectively, and C is a physical parameter) we prove both long time stability and the rate of convergence of the method. Some numerical experiments are presented and discussed.

A time-spectral approach to numerical weather prediction

Finite difference methods are traditionally used for modelling the time domain in numerical weather prediction (NWP). Time-spectral solution is an attractive alternative for reasons of accuracy and efficiency and because time step limitations associated with causal CFL-like criteria, typical for explicit finite difference methods, are avoided. In this work, the Lorenz 1984 chaotic equations are solved using the time-spectral algorithm GWRM (Generalized Weighted Residual Method). Comparisons of accuracy and efficiency are carried out for both explicit and implicit time-stepping algorithms. It is found that the efficiency of the GWRM compares well with these methods, in particular at high accuracy. For perturbative scenarios, the GWRM was found to be as much as four times faster than the finite difference methods. A primary reason is that the GWRM time intervals typically are two orders of magnitude larger than those of the finite difference methods. The GWRM has the additional advantage to produce analytical solutions in the form of Chebyshev series expansions. The results are encouraging for pursuing further studies, including spatial dependence, of the relevance of time-spectral methods for NWP modelling. Program summaryProgram Title: Time-adaptive GWRM Lorenz 1984Program Files doi:http://dx.doi.org/10.17632/4nxfyjj7nv.1Licensing provisions: MITProgramming language: MapleNature of problem: Ordinary differential equations with varying degrees of complexity are routinely solved with numerical methods. The set of ODEs pertaining to chaotic systems, for instance those related to numerical weather prediction (NWP) models, are highly sensitive to initial conditions and unwanted errors. To accurately solve ODEs such as the Lorenz equations (E. N. Lorenz, Tellus A 36 (1984) 98–110), small time steps are required by traditional time-stepping methods, which can be a limiting factor regarding the efficiency, accuracy, and stability of the computations.Solution method: The Generalized Weighted Residual Method, being a time-spectral algorithm, seeks to increase the time intervals in the computation without degrading the efficiency, accuracy, and stability. It does this by postulating a solution ansatz as a sum of weighted Chebyshev polynomials, in combination with the Galerkin method, to create a set of linear/non-linear algebraic equations. These algebraic equations are then solved iteratively using a Semi Implicit Root solver (SIR), which has been chosen due to its enhanced global convergence properties. Furthermore, to achieve a desired accuracy across the entire domain, a time-adaptive algorithm has been developed. By evaluating the magnitudes of the Chebyshev coefficients in the time dimension of the solution ansatz, the time interval can either be decreased or increased.

A quasi-Lagrangian finite element method for the Navier–Stokes equations in a time-dependent domain

The paper develops a finite element method for the Navier–Stokes equations of incompressible viscous fluid in a time-dependent domain. The method builds on a quasi-Lagrangian formulation of the problem. The paper provides stability and convergence analysis of the fully discrete (finite-difference in time and finite-element in space) method. The analysis does not assume any CFL time-step restriction, it rather needs mild conditions of the form Δt≤C, where C depends only on problem data, and h2mu+2≤cΔt, mu is polynomial degree of velocity finite element space. Both conditions result from a numerical treatment of practically important non-homogeneous boundary conditions. The theoretically predicted convergence rate is confirmed by a set of numerical experiments. Further we apply the method to simulate a flow in a simplified model of the left ventricle of a human heart, where the ventricle wall dynamics is reconstructed from a sequence of contrast enhanced computed tomography images.

Corrected implicit Monte Carlo

In this work we develop a set of nonlinear correction equations to enforce a consistent time-implicit emission temperature for the original semi-implicit IMC equations. We present two possible forms of correction equations: one results in a set of non-linear, zero-dimensional, non-negative, explicit correction equations, and the other results in a non-linear, non-negative, Boltzman transport correction equation. The zero-dimensional correction equations adheres to the maximum principle for the material temperature, regardless of frequency-dependence, but does not prevent maximum principle violation in the photon intensity, eventually leading to material overheating. The Boltzman transport correction guarantees adherence to the maximum principle for frequency-independent simulations, at the cost of evaluating a reduced source non-linear Boltzman equation. We present numerical evidence suggesting that the Boltzman transport correction, in its current form, significantly improves time step limitations but does not guarantee adherence to the maximum principle for frequency-dependent simulations.

Asymptotic-Preserving Monte Carlo Methods for Transport Equations in the Diffusive Limit

We develop a new Monte Carlo method that solves hyperbolic transport equations with stiff terms, characterized by a (small) scaling parameter. In particular, we focus on systems which lead to a reduced problem of parabolic type in the limit when the scaling parameter tends to zero. Classical Monte Carlo methods suffer severe time step limitations in these situations, due to the fact that the characteristic speeds go to infinity in the diffusion limit. This makes the problem a real challenge, since the scaling parameter may differ by several orders of magnitude in the domain. To circumvent these time step limitations, we construct a new, asymptotic-preserving Monte Carlo method that is stable independently of the scaling parameter and degenerates to a standard probabilistic approach for solving the limiting equation in the diffusion limit. The method uses an implicit time discretization to formulate a modified equation in which the characteristic speeds do not grow indefinitely when the scaling factor tends to zero. The resulting modified equation can readily be discretized by a Monte Carlo scheme, in which the particles combine a finite propagation speed with a time step dependent diffusion term. We show the performance of the method by comparing it with standard (deterministic) approaches in the literature.

A Crank--Nicolson Finite Element Method and the Optimal Error Estimates for the Modified Time-Dependent Maxwell--Schrödinger Equations

In this paper, we study a Crank--Nicolson finite element method for the modified Maxwell--Schrodinger equations, which are deduced from the time-dependent Maxwell--Schrodinger equations under some assumptions. We present an (almost) unconditionally optimal error estimate for the numerical algorithm which only requires the time step $\Delta t < 1$. The key to our analysis is twofold. For one thing, we use some tricks to tackle the tough nonlinear terms without applying the inverse inequality to estimate the $L^{\infty}$ norm of the numerical solutions and get rid of certain restrictions like $\Delta t \leq C \, h^{\alpha}$ on the time step. For another, we employ a new technique to weaken the time step restriction required by the application of the discrete Gronwall inequality. Numerical tests are carried out to verify the theoretical results.

An unconditionally stable technique for uncertainty assessment in random media based on the ADI scheme

Stochastic media problems present a severe impact on the properties of the electromagnetic fields that usually deterministic methods tend to overlook. In this paper, we present a finite-difference time-domain (FDTD) technique, which is free from time step restrictions and manages to efficiently calculate the statistics of the involved fields in one realization. The proposed algorithm is a variation of the Stochastic FDTD method and exhibits satisfactory outcomes in high density grids. Numerical results verify the validity of our scheme compared to other methods, such as the Stochastic FDTD and the Monte Carlo techniques.

A comparative study of discrete velocity methods for low-speed rarefied gas flows

In the study of rarefied gas dynamics, the discrete velocity method (DVM) has been widely employed to solve the gas kinetic equations. Although various versions of DVM have been developed, their performance, in terms of modeling accuracy and computational efficiency, is yet to be comprehensively studied in all the flow regimes. Here, the traditional third-order time-implicit Godunov DVM (GDVM) and the recently developed discrete unified gas-kinetic scheme (DUGKS) are analysed in finding steady-state solutions of the low-speed force-driven Poiseuille and lid-driven cavity flows. With the molecular collision and free streaming being treated simultaneously, the DUGKS preserves the second-order accuracy in the spatial and temporal discretizations in all flow regimes. Towards the hydrodynamic flow regime, not only is the DUGKS faster than the GDVM when using the same spatial mesh, but also requires less spatial resolution than that of the GDVM to achieve the same numerical accuracy. From the slip to free molecular flow regimes, however, the DUGKS is slower than the GDVM, due to the complicated flux evaluation and the restrictive time step which is smaller than the maximum effective time step of the GDVM. Therefore, the DUGKS is preferable for problems involving different flow regimes, particularly when the hydrodynamic flow regime is dominant. For highly rarefied gas flows, if the steady-state solution is mainly concerned, the implicit GDVM, which can boost the convergence significantly, is a better choice.

A stabilized finite volume method for the evolutionary Stokes–Darcy system

In this paper, we propose a stabilized fully discrete finite volume method based on two local Gauss integrals for a non-stationary Stokes–Darcy problem. This stabilized method is free of stabilized parameters and uses the lowest equal-order finite element triples P1–P1–P1 for approximating the velocity, pressure and hydraulic head of the Stokes–Darcy model. Under a modest time step restriction in relation to physical parameters, we give the stability analysis and the error estimates for the stabilized finite volume scheme by means of a relationship between finite volume and finite element approximations with the lower order elements. Finally, a series of numerical experiments are provided to demonstrate the validity of the theoretical results.

On the stability of projection methods for the incompressible Navier–Stokes equations based on high-order discontinuous Galerkin discretizations

The present paper deals with the numerical solution of the incompressible Navier–Stokes equations using high-order discontinuous Galerkin (DG) methods for discretization in space. For DG methods applied to the dual splitting projection method, instabilities have recently been reported that occur for small time step sizes. Since the critical time step size depends on the viscosity and the spatial resolution, these instabilities limit the robustness of the Navier–Stokes solver in case of complex engineering applications characterized by coarse spatial resolutions and small viscosities. By means of numerical investigation we give evidence that these instabilities are related to the discontinuous Galerkin formulation of the velocity divergence term and the pressure gradient term that couple velocity and pressure. Integration by parts of these terms with a suitable definition of boundary conditions is required in order to obtain a stable and robust method. Since the intermediate velocity field does not fulfill the boundary conditions prescribed for the velocity, a consistent boundary condition is derived from the convective step of the dual splitting scheme to ensure high-order accuracy with respect to the temporal discretization. This new formulation is stable in the limit of small time steps for both equal-order and mixed-order polynomial approximations. Although the dual splitting scheme itself includes inf–sup stabilizing contributions, we demonstrate that spurious pressure oscillations appear for equal-order polynomials and small time steps highlighting the necessity to consider inf–sup stability explicitly.