Time Step Constraint Research Articles

Large Eddy Simulation of a Free-Burning Arc Discharge in Argon with a Turbulent Cross Flow and External Field

In most of the numerical approaches proposed for modeling high-intensity plasma-arcs, the effects of turbulence on the arc structure are often excluded because of the intricate physics originating from the interaction of turbulent scales, high-temperature gas dynamics, magnetohydrodynamics (MHD) and chemical kinetics. The goal of this study is threefold: to develop a generic turbulent MHD model to simulate free-burning arc discharges, to validate the code with available experimental data, and to investigate the effect of an external field and turbulent cross flow on the free-burning arc configuration. The governing equations are solved in conservative form using a hybrid scheme that combines a high-order monotonic upwind scheme with a second-order central scheme. The fluid and MHD turbulence are resolved using a large eddy simulation (LES) approach with a recently developed sub-grid closure model. An implicit scheme is used to compute the magnetic diffusion term appearing in the magnetic induction equation to alleviate the severe time-step constraint. The comparison of the model prediction with experimental data for Argon arcs at different current intensities shows generally good agreement. When an external field is applied, the overall shape of the free-burning arc drastically changes. The straightening of the arc indicates the potential for stabilization of a free-burning arc by magnetic forces. Even though the turbulence is significantly attenuated as a result of the thermal expansion near the cathode, it adds an unsteady characteristic to the arc and, in general, has a negative impact on the stabilization of the electrical discharge.

Discrete gradient methods for preserving a first integral of an ordinary differential equation

In this paper we consider discrete gradient methods for approximating the solution and preserving a first integral (also called a constant of motion) of autonomous ordinary differential equations. We prove under mild conditions for a large class of discrete gradient methods that the numerical solution exists and is locally unique, and that for arbitrary $p\in \mathbb{N}$ we may construct a method that is of order $p$. In the proofs of these results we also show that the constants in the time step constraint and the error bounds may be chosen independently from the distance to critical points of the first integral. &nbsp In the case when the first integral is quadratic, for arbitrary $p \in \mathbb{N}$, we have devised a new method that is linearly implicit at each time step and of order $p$. A numerical example suggests that this new method has advantages in terms of efficiency.

A parallel local timestepping Runge–Kutta discontinuous Galerkin method with applications to coastal ocean modeling

Geophysical flows over complex domains often encompass both coarse and highly resolved regions. Approximating these flows using shock-capturing methods with explicit timestepping gives rise to a Courant–Friedrichs–Lewy (CFL) timestep constraint. This approach can result in small global timesteps often dictated by flows in small regions, vastly increasing computational effort over the whole domain. One approach for coping with this problem is to use locally varying timesteps. In previous work, we formulated a local timestepping (LTS) method within a Runge–Kutta discontinuous Galerkin framework and demonstrated the accuracy and efficiency of this method on serial machines for relatively small-scale shallow water applications. For more realistic models involving large domains and highly complex physics, the LTS method must be parallelized for multi-core parallel computers. Furthermore, additional physics such as strong wind forcing can effect the choice of local timesteps. In this paper, we describe a parallel LTS method, parallelized using domain decomposition and MPI. We demonstrate the method on tidal flows and hurricane storm surge applications in the coastal regions of the Western North Atlantic Ocean.

An asymptotic preserving scheme for the streamer simulation

In this work we describe the development of a robust simulator of the streamer phenomenon. A very common problem of computational plasma physics is the stability of the numerical scheme when the quasi-neutral regime is approached. Actually the stability condition may imply very tight time step constraints. We prove, both numerically and analytically, the good behavior of our approach when the system approaches the quasi-neutral regime. In particular, the semi-implicit method we employ requires mild time step restrictions to guarantee the stability. Some other properties such as the positivity of the scheme are investigated and a detailed review of the numerical implementation is included.

A conservative numerical method for the Cahn–Hilliard equation with Dirichlet boundary conditions in complex domains

In this paper we present a conservative numerical method for the Cahn–Hilliard equation with Dirichlet boundary conditions in complex domains. The method uses an unconditionally gradient stable nonlinear splitting numerical scheme to remove the high-order time-step stability constraints. The continuous problem has the conservation of mass and we prove the conservative property of the proposed discrete scheme in complex domains. We describe the implementation of the proposed numerical scheme in detail. The resulting system of discrete equations is solved by a nonlinear multigrid method. We demonstrate the accuracy and robustness of the proposed Dirichlet boundary formulation using various numerical experiments. We numerically show the total energy decrease and the unconditionally gradient stability. In particular, the numerical results indicate the potential usefulness of the proposed method for accurately calculating biological membrane dynamics in confined domains.

Implementation of Point-Implicit Diffusion Operator in a Moving-Particle Method (MPPM)

A point-implicit strategy to circumvent the time-step constraint for diffusion process in a moving-particle method (MPPM) is proposed in this article. Since the Eulerian mesh has been established to deal with all pressure-related operators, it is also employed in constructing an accurate and effective local mesh (LM) formulation to represent the diffusion operator on the Lagranian particles. Explicit evaluation of the diffusion process will become inefficient due to the limited time step with a refined particle distribution. On the other hand, a fully implicit method requires an iterative procedure which will inevitably increase the computational burden. In the proposed point-implicit strategy, only the solution value at the objective particle is handled implicitly, and those at neighboring ones are treated explicitly. The time-step constraint on the diffusion process can then be taken out. All these efforts will result in an effective solution procedure for incompressible flow simulations.

A comparison study of phase-field models for an immiscible binary mixture with surfactant

The Ginzburg-Landau free energy functional with two order parameters has been widely used to describe surfactant adsorption phenomena at the interface between two immiscible fluids such as oil and water. To model surfactant adsorption, additional surfactant related terms are added to the original free energy functional which models an immiscible binary mixture. In this paper, we present a detailed comparison of phase-field models for an immiscible binary mixture with surfactant. In particular, we investigate the effects of mathematical model parameters on equilibrium surfactant profile across the interface between the immiscible binary mixture. Most previous models have severe time-step constraints due to the nonlinear coupling of order parameters. To solve these stability problems, we propose a special case of these models which allows the use of a much larger time-step size. We also apply a type of unconditionally gradient stable scheme and a fast multigrid method to solve the proposed model efficiently and accurately.

A conservative orbital advection scheme for simulations of magnetized shear flows with the PLUTO code

Explicit numerical computations of super-fast differentially rotating disks are subject to the time-step constraint imposed by the Courant condition. When the bulk orbital velocity largely exceeds any other wave speed the time step is considerably reduced and a large number of steps may be necessary to complete the computation. We present a robust numerical scheme to overcome the Courant limitation by extending the algorithm previously known as FARGO (Fast Advection in Rotating Gaseous Objects) to the equations of magnetohydrodynamics (MHD). The proposed scheme conserves total angular momentum and energy to machine precision and works in Cartesian, cylindrical, or spherical coordinates. The algorithm is implemented in the PLUTO code for astrophysical gasdynamics and is suitable for local or global simulations of accretion or proto-planetary disk models. By decomposing the total velocity into an average azimuthal contribution and a residual term, the algorithm solves the MHD equations through a linear transport step in the orbital direction and a standard nonlinear solver applied to the MHD equations written in terms of the residual velocity. Since the former step is not subject to any stability restriction, the Courant condition is computed only in terms of the residual velocity, leading to substantially larger time steps. The magnetic field is advanced in time using the constrained transport method in order to preserve the divergence-free condition. Conservation of total energy and angular momentum is enforced at the discrete level by properly expressing the source terms in terms of upwind fluxes available during the standard solver. Our results show that applications of the proposed orbital-advection scheme to problems of astrophysical relevance provides, at reduced numerical cost, equally accurate and less dissipative results than standard time-marching schemes.

Local time-stepping in Runge–Kutta discontinuous Galerkin finite element methods applied to the shallow-water equations

Large geophysical flows often encompass both coarse and highly resolved regions. Approximating these flows using shock-capturing methods with explicit time stepping gives rise to a Courant–Friedrichs–Lewy (CFL) time step constraint. Even if the refined regions are sparse, they can restrict the global CFL condition to very small time steps, vastly increasing computational effort over the whole domain. One method to cope with this problem is to use locally varying time steps over the domain. These are also referred to as multi-rate methods in the ODE literature. Ideally, such methods must be conservative, accurate and easy to implement. In this study, we derive a second-order, local time stepping procedure within a Runge–Kutta discontinuous Galerkin (RKDG) framework to solve the shallow water equations. This procedure is based on previous first-order work of the second author and collaborator Kirby [1–3]. As we are interested in both coastal and overland flows due to, e.g., rainfall, wetting and drying is incorporated into the model. Numerical results are shown, which verify the accuracy and efficiency of the approach (compared to using a globally defined CFL time step), and the application of the method to rainfall–runoff scenarios.

An implicit free surface algorithm for geodynamical simulations

Identifying the dominant controls on Earth’s surface topography is of critical importance to understanding both the short- and long-term evolution of geological processes and past- and present-day dynamics of Earth’s coupled mantle–lithosphere system. The ability to simulate a stress free — or a so-called ‘free surface’ — boundary condition is required to examine such processes via numerical models. However, at present, geodynamical models incorporating a free surface are limited, as most underlying free surface algorithms place severe restrictions on the computational timestep. Consequently, the simulations are often intractable. In this study, we introduce a new approach for incorporating a free surface within geodynamical models: an algorithm, in which free surface elevation is treated as an independent variable and is solved for in conjunction with the momentum and continuity equation, using implicit time integration. We demonstrate that the method is straightforward to implement in existing models and, using a series of analytical and benchmark comparisons, we show that it does not suffer from the timestep constraints of previous schemes. Furthermore, the scheme can be made second order accurate in time, at no additional cost. The method therefore dramatically improves the computational efficiency of geodynamical simulations including a free surface, whilst maintaining solution accuracy.

Semi-Implicit Runge–Kutta Schemes: Development and Application to Compressible Channel Flow

Two semi-implicit six-stage Runge–Kutta algorithms are developed for the simulation of wall-bounded flows. Using these schemes, time integration is implicit in the wall-normal direction, and explicit in the other directions, to relax the time step constraint due to the fine mesh near the wall. The explicit subscheme is a six-stage fourth-order low-storage Runge–Kutta scheme. Based on analysis in Fourier space and results obtained for propagation test cases, the semi-implicit schemes are shown to be of order 3 and, for waves discretized by a number of points per period between 4 and 16, to be as accurate as, or more accurate than, the standard explicit fourth-order Runge–Kutta algorithm in terms of dissipation and dispersion. The large-eddy simulation of a compressible turbulent channel flow at a friction Reynolds number of 360 and a Mach number of 0.1 is then carried out with one of the proposed algorithms. The computational time is reduced by a factor 1.33 with respect to a large-eddy simulation using the explicit subscheme in all directions. Wall-pressure and velocity spectra from the large-eddy simulation are presented to give insights into the flow turbulent structures. In particular, wave number–frequency spectra are calculated. Acoustic components appear to be identified in these spectra.

A method for overcoming the surface tension time step constraint in multiphase flows II

SUMMARYA method for overcoming the surface tension time step constraint is presented. The algorithm presented in this work is an improvement on the work presented by Sussman and Ohta (SIAM J Sci Comput 2009). In this work, the method of Sussman and Ohta is extended in order to treat problems with contact angle dynamics. Furthermore, this work presents a more efficient method for computing volume‐preserving motion by mean curvature than the method presented previously. The new method is tested on the following four 2D problems: (1) 3D axisymmetric (r−z) surface tension driven zero gravity droplet oscillation, (2) measurement of the magnitude of parasitic currents for a droplet on a substrate initialized in static equilibrium, (3) relaxation of a 2D droplet on a substrate to static shape, and (3) relaxation of a 2D bubble on a substrate to static shape. Copyright © 2011 John Wiley & Sons, Ltd.

A High Order Sharp-Interface Method with Local Time Stepping for Compressible Multiphase Flows

AbstractIn this paper, a new sharp-interface approach to simulate compressible multiphase flows is proposed. The new scheme consists of a high order WENO finite volume scheme for solving the Euler equations coupled with a high order path-conservative discontinuous Galerkin finite element scheme to evolve an indicator function that tracks the material interface. At the interface our method applies ghost cells to compute the numerical flux, as the ghost fluid method. However, unlike the original ghost fluid scheme of Fedkiw et al. [15], the state of the ghost fluid is derived from an approximate-state Riemann solver, similar to the approach proposed in [25], but based on a much simpler formulation. Our formulation leads only to one single scalar nonlinear algebraic equation that has to be solved at the interface, instead of the system used in [25]. Away from the interface, we use the new general Osher-type flux recently proposed by Dumbser and Toro [13], which is a simple but complete Riemann solver, applicable to general hyperbolic conservation laws. The time integration is performed using a fully-discrete one-step scheme, based on the approaches recently proposed in [5,7]. This allows us to evolve the system also with time-accurate local time stepping. Due to the sub-cell resolution and the subsequent more restrictive time-step constraint of the DG scheme, a local evolution for the indicator function is applied, which is matched with the finite volume scheme for the solution of the Euler equations that runs with a larger time step. The use of a locally optimal time step avoids the introduction of excessive numerical diffusion in the finite volume scheme. Two different fluids have been used, namely an ideal gas and a weakly compressible fluid modeled by the Tait equation. Several tests have been computed to assess the accuracy and the performance of the new high order scheme. A verification of our algorithm has been carefully carried out using exact solutions as well as a comparison with other numerical reference solutions. The material interface is resolved sharply and accurately without spurious oscillations in the pressure field.

Compensated Monte Carlo Collision Model for Particle-in-Cell Simulation in High-Pressure Plasmas

The Monte Carlo collision (MCC) model is widely adopted to simulate discharge plasmas using the particle-in-cell (PIC) method; however, it has low efficiency in high-pressure plasmas because of the small time steps required due to the constraint of high collision frequency. To relax this time step constraint, a compensated Monte Carlo collision model (CMCC) is proposed which considers multiple collisions in a time step as a series of single collisions to compensate for the neglected collisions. The electron motion in a high-pressure He gas for various reduced electric fields E/N and the streamer formation process in a laser-triggered spark gap were simulated using the CMCC model. Simulation results showed that the CMCC model with long time step obtained reasonable electron velocity distribution, temperature, drift velocity, plasma density, and space-charge field. It was demonstrated that the CMCC model had high accuracy and high efficiency, particularly for PIC simulation in high-pressure plasmas.

A multiscale approach to mesh-based surface tension flows

We present an approach to simulate flows driven by surface tension based on triangle meshes. Our method consists of two simulation layers: the first layer is an Eulerian method for simulating surface tension forces that is free from typical strict time step constraints. The second simulation layer is a Lagrangian finite element method that simulates sub-grid scale wave details on the fluid surface. The surface wave simulation employs an unconditionally stable, symplectic time integration method that allows for a high propagation speed due to strong surface tension. Our approach can naturally separate the grid- and sub-grid scales based on a volume-preserving mean curvature flow. As our model for the sub-grid dynamics enforces a local conservation of mass, it leads to realistic pinch off and merging effects. In addition to this method for simulating dynamic surface tension effects, we also present an efficient non-oscillatory approximation for capturing damped surface tension behavior. These approaches allow us to efficiently simulate complex phenomena associated with strong surface tension, such as Rayleigh-Plateau instabilities and crown splashes, in a short amount of time.

IMPLEMENTATION OF THE SHEARING BOX APPROXIMATION IN ATHENA

We describe the implementation of the shearing box approximation for the study of the dynamics of accretion disks in the Athena magnetohydrodynamics (MHD) code. Second-order Crank-Nicholson time differencing is used for the Coriolis and tidal gravity source terms that appear in the momentum equation for accuracy and stability. We show this approach conserves energy for epicyclic oscillations in hydrodynamic flows to round-off error. In the energy equation, the tidal gravity source terms are differenced as the gradient of an effective potential in a way which guarantees that total energy (including the gravitational potential energy) is also conserved to round-off error. We introduce an orbital advection algorithm for MHD based on constrained transport to preserve the divergence-free constraint on the magnetic field. This algorithm removes the orbital velocity from the time step constraint, and makes the truncation error more uniform in radial position. Modifications to the shearing box boundary conditions applied at the radial boundaries are necessary to conserve the total vertical magnetic flux. In principle similar corrections are also required to conserve mass, momentum and energy, however in practice we find the orbital advection method conserves these quantities to better than 0.03% over hundreds of orbits. The algorithms have been applied to studies of the nonlinear regime of the MRI in very wide (up to 32 scale heights) horizontal domains.

A non-stiff boundary integral method for 3D porous media flow with surface tension

We present an efficient, non-stiff boundary integral method for 3D porous media flow with surface tension. Surface tension introduces high order (i.e., high derivative) terms into the evolution equations, and this leads to severe stability constraints for explicit time-integration methods. Furthermore, the high order terms appear in non-local operators, making the application of implicit methods difficult. Our method uses the fundamental coefficients of the surface as dynamical variables, and employs a special isothermal parameterization of the interface which enables efficient application of implicit or linear propagator time-integration methods via a small-scale decomposition. The method is tested by computing the relaxation of an interface to a flat surface under the action of surface tension. These calculations employ an approximate interface velocity to test the stiffness reduction of the method. The approximate velocity has the same mathematical form as the exact velocity, but avoids the numerically intensive computation of the full Birkhoff–Rott integral. The algorithm is found to be effective at eliminating the severe time-step constraint that plagues explicit time-integration methods.

Dynamic simulation of locally inextensible vesicles suspended in an arbitrary two-dimensional domain, a boundary integral method

We consider numerical algorithms for the simulation of hydrodynamics of two-dimensional vesicles suspended in a viscous Stokesian fluid. The motion of vesicles is governed by the interplay between hydrodynamic and elastic forces. Continuum models of vesicles use a two-phase fluid system with interfacial forces that include tension (to maintain local “surface” inextensibility) and bending. Vesicle flows are challenging to simulate. On the one hand, explicit time-stepping schemes suffer from a severe stability constraint due to the stiffness related to high-order spatial derivatives in the bending term. On the other hand, implicit time-stepping schemes can be expensive because they require the solution of a set of nonlinear equations at each time step. Our method is an extension of the work of Veerapaneni et al. [S.K. Veerapaneni, D. Gueyffier, D. Zorin, G. Biros, A boundary integral method for simulating the dynamics of inextensible vesicles suspended in a viscous fluid in 2D, Journal of Computational Physics 228(7) (2009) 2334–2353], in which a semi-implicit time-marching scheme based on a boundary integral formulation of the Stokes problem for vesicles in an unbounded medium was proposed. In this paper, we consider two important generalizations: (i) confined flows within arbitrary-shaped stationary/moving geometries; and (ii) flows in which the interior (to the vesicle) and exterior fluids have different viscosity. In the rest of the paper, we will refer to this as the “ viscosity contrast”. These two problems require solving additional integral equations and cause nontrivial modifications to the previous numerical scheme. Our method does not have severe time-step stability constraints and its computational cost-per-time-step is comparable to that of an explicit scheme. The discretization is pseudo-spectral in space, and multistep BDF in time. We conduct numerical experiments to investigate the stability, accuracy and the computational cost of the algorithm. Overall, our method achieves several orders of magnitude speed-up compared to standard explicit schemes. As a preliminary validation of our scheme, we study the dependence of the inclination angle of a single vesicle in shear flow on the viscosity contrast and the reduced area of the vesicle, the lateral migration of vesicles in shear flow, the dispersion of two vesicles, and the effective viscosity of a dilute suspension of vesicles.

An Investigation of Numerical Modeling of Transient Heat Conduction in a One-Dimensional Functionally Graded Material

This article presents a closed form analytical solution for one-dimensional transient heat conduction in a material where the thermal conductivity varies linearly through the thickness but the thermal diffusivity is held constant. This solution is used to validate the results from finite-difference and finite-element approximations that account for this variation at the element level. This was motivated by a suggested limitation on the minimum time step used in the commercial finite-element software code ABAQUS for quadratic elements. Good agreement was found between the analytical and numerical approximations, indicating that conventional numerical techniques may be sufficiently robust to analyze heat conduction problems in functionally graded materials without the use of special elements. The minimum time step constraint was found to be unnecessary for a convective boundary condition for the one-dimensional elements and property variation used in this study.

Parallel Galerkin domain decomposition procedures for wave equation

Parallel Galerkin domain decomposition procedures for wave equation are given. These procedures use implicit method in the sub-domains and simple explicit flux calculation on the inter-boundaries of sub-domains by integral mean method or extrapolation method. Thus, the parallelism can be achieved by these procedures. The explicit nature of the flux prediction induces a time step constraint that is necessary to preserve the stability. L 2 -norm error estimates are derived for these procedures. Experimental results are presented to confirm the theoretical results.