Large Time Steps Research Articles

In-Timestep Remeshing for Contacting Elastodynamics

We propose In-Timestep Remeshing, a fully coupled, adaptive meshing algorithm for contacting elastodynamics where remeshing steps are tightly integrated, implicitly, within the timestep solve. Our algorithm refines and coarsens the domain automatically by measuring physical energy changes within each ongoing timestep solve. This provides consistent, degree-of-freedom-efficient, productive remeshing that, by construction, is physics-aware and so avoids the errors, over-refinements, artifacts, per-example hand-tuning, and instabilities commonly encountered when remeshing with timestepping methods. Our in-timestep computation then ensures that each simulation step's output is both a converged stable solution on the updated mesh and a temporally consistent trajectory with respect to the model and solution of the last timestep. At the same time, the output is guaranteed safe (intersection- and inversion-free) across all operations. We demonstrate applications across a wide range of extreme stress tests with challenging contacts, sharp geometries, extreme compressions, large timesteps, and wide material stiffness ranges - all scenarios well-appreciated to challenge existing remeshing methods.

RAM: Rapid Advection Algorithm on Arbitrary Meshes

The study of many astrophysical flows requires computational algorithms that can capture high Mach number flows, while resolving a large dynamic range in spatial and density scales. In this paper we present a novel method, RAM: Rapid Advection Algorithm on Arbitrary Meshes. RAM is a time-explicit method to solve the advection equation in problems with large bulk velocity on arbitrary computational grids. In comparison with standard upwind algorithms, RAM enables advection with larger time steps and lower truncation errors. Our method is based on the operator splitting technique and conservative interpolation. Depending on the bulk velocity and resolution, RAM can decrease the numerical cost of hydrodynamics by more than one order of magnitude. To quantify the truncation errors and speed-up with RAM, we perform one- and two-dimensional hydrodynamics tests. We find that the order of our method is given by the order of the conservative interpolation and that the effective speed-up is in agreement with the relative increment in time step. RAM will be especially useful for numerical studies of disk−satellite interaction, characterized by high bulk orbital velocities and nontrivial geometries. Our method dramatically lowers the computational cost of simulations that simultaneously resolve the global disk and potential well inside the Hill radius of the secondary companion.

Multi-type dynamic load identification algorithm in continuous system: A numerical and experimental study based on SSM-Newmark-β

Dynamic load identification based on structural responses is an important problem in the field of engineering and plays an important role in the condition assessment of mechanical structures. Current popular load identification methods such as the state-space method (SSM) and the Green's kernel function method (GKFM) are implemented on discrete systems with outstanding performance, but when dealing with complex continuous systems, there are some limitations such as low efficiency and inaccuracy. In this paper, a continuous system load identification algorithm based on SSM and Newmark-β is proposed for the first time. Using modal coordinate transformation and modal truncation methods, the number of infinite vibration differential equations of a continuous system in physical space are converted to a finite number of vibration differential equations in modal space, where the modal truncation order and the optimal layout of the sensors are combined by modal strain energy. The Newmark-β method in modal space is derived and thus combined with SSM to obtain a load identification model Y=HF called SSM-Newmark-β method, which reduces the size of the transfer matrix H. The solution time of the proposed algorithm is 0.38 s for continuous systems, which is rather shorter than that for discrete systems. Furthermore, the simulations show that it has better accuracy and noise immunity than SSM and GKFM in the identification of sinusoidal load, impact load, and random load. The effect of time step is also discussed which reveals that the larger time step has less effect on the proposed algorithm. An experimental study is carried out on a cantilever beam system. The result verifies that the SSM-Newmark-β algorithm has better accuracy in load identification for the continuous system. This research provides a new sight for the real-time identification of dynamic loads for complex structures.

Lattice Boltzmann solidification modeling of forced convection internal flows applied to Gen-IV nuclear reactor coolants

Generation-IV nuclear reactors involve novel coolants such as molten salts and lead. One of the issues that must be assessed with these coolants is the risk of solidification during reactor operation owing to their high melting temperature. However, the numerical simulation of solidification problems is computationally demanding. In this work, we assess the applicability of Lattice Boltzmann Method (LBM) to nuclear coolants solidification problems under laminar internal forced flow conditions. The double distribution Total Enthalpy-Partially Saturated Method is implemented in the open-source code OpenLB and compared against more standard Finite-Volume Computational Fluid Dynamics (FV-CFD) results, which uses the enthalpy-porosity Eulerian Multiphase method, in commercial code Star-CCM+. The test cases comprise the solidification of two fluids with significantly different Prandtl numbers, molten salt and lead, in a 3D pipe geometry in laminar flow conditions with an initial Reynolds number of 100. The solid interface advances from the cooled wall towards the bulk of the pipe reducing the area of flow passage as the transient evolves. Results of LBM velocity fields, solid fraction profiles, temperature fields, and pressure fields are in good agreement with FV-CFD results for both coolants. Furthermore, LBM stability and accuracy limits are independently explored for molten salts and lead using Two-Relaxation Time (TRT) and Bhatnagar-Gross-Krook (BGK) collision operators. Simulation time between FV-CFD and LBM is also contrasted for both coolants. In this work, LBM solvers use explicit time schemes with a regular grid, whereas FV-CFD solutions are able to include implicit time schemes and near-wall mesh refinement. For molten salt flows (high Prandtl number) FV-CFD allows for considerably larger timesteps than LBM, which is limited by the maximum lattice Mach number limitation and stability considerations. For lead flows (low Prandtl number), the implicit scheme in FV-CFD is limited in terms of accuracy for large timesteps due to the limiting Fourier-number. On the other hand, LBM does not present notable stability issues for lead and the lattice Mach number limit is not as relevant as for the molten salt case for this Reynolds number, allowing for faster simulation times than FV-CFD.

Hybridizing FDM and FVM scheme of high-precision interface fast capture for mixed free-surface-pressurized flow in large cascade water delivery system

In recent years, China has been building several inter-basin water conveyance projects across mountains and deep valleys, leading to the extensive use of long tunnels and inverted siphons. The dynamics of mixed free-surface-pressurized flows are critical for tunnel design and operational safety. However, traditional numerical computation schemes cannot precisely and efficiently capture the pressure interface because the tunnel may be more than 70 km long. This study aims to develop a hybrid scheme that is as fast as the large-time-step Preissmann four-point scheme (FDM) and has the approximate interface accuracy of the finite volume method (FVM) by dynamic grid meshing. A dynamic mesh domain model (DM) is proposed by adopting an FVM mesh with a small time step to dynamically capture the interface and applying an FDM mesh with a large time step to improve computational efficiency. The results show that the method can simulate flow patterns, such as transcritical and pressurized flows, by ignoring the acceleration convection terms and capturing mixed free-surface-pressurized flows conveniently, accurately and efficiently. Furthermore, it can accelerate the computational speed of the transient mixed flow by a factor of approximately 100 when the target tunnel length exceeds 1200 m. The proposed scheme cannot capture the water hammer pressure because of the large time step. However, it can be effectively utilized in large cascade water delivery systems where the flow changes gradually.

Time‐parallel integration and phase averaging for the nonlinear shallow‐water equations on the sphere

We describe a proof‐of‐concept development and application of a phase‐averaging technique to the nonlinear rotating shallow‐water equations on the sphere, discretised using compatible finite‐element methods. Phase averaging consists of averaging the nonlinearity over phase shifts in the exponential of the linear wave operator. Phase averaging aims to capture the slow dynamics in a solution that is smoother in time (in transformed variables), so that larger timesteps may be taken. We overcome the two key technical challenges that stand in the way of studying the phase averaging and advancing its implementation: (1) we have developed a stable matrix exponential specific to finite elements and (2) we have developed a parallel finite averaging procedure. Following recent studies, we consider finite‐width phase‐averaging windows, since the equations have a finite timescale separation. In our numerical implementation, the averaging integral is replaced by a Riemann sum, where each term can be evaluated in parallel. This creates an opportunity for parallelism in the timestepping method, which we use here to compute our solutions. Here, we focus on the stability and accuracy of the numerical solution. We confirm that there is an optimal averaging window, in agreement with theory. Critically, we observe that the combined time discretisation and averaging error is much smaller than the time discretisation error in a semi‐implicit method applied to the same spatial discretisation. An evaluation of the parallel aspects will follow in later work.

A compact ADI finite difference method for 2D reaction–diffusion equations with variable diffusion coefficients

Reaction–diffusion systems on a spatially heterogeneous domain have been widely used to model various biological applications. However, solving such partial differential equations (PDEs) analytically is rarely possible. Therefore, efficient and accurate numerical methods for solving such PDEs are desired. In this paper, we apply the well-known Padé approximation-based operator splitting techniques and develop a fourth-order compact alternative directional implicit (ADI) scheme. The new scheme is compact and fourth-order accurate in space. Combined with the Richardson extrapolation, the method can be improved to fourth-order accuracy in time. Stability analysis shows that the method is unconditionally stable; thus, a large time step can be used to improve the overall computational efficiency. Numerical examples have also demonstrated the new scheme’s high efficiency and high order accuracy.

An optimal spatial-filtering method derived from eigenvalue perturbation for extending the Courant-Friedrichs-Lewy stability limit

Explicit time-marching schemes are widely used in numerical simulations. However, their maximum time step is constrained by the Courant-Friedrichs-Lewy (CFL) stability limit. Two methods have been developed to extend this limit: (1) the eigenvalue perturbation method, which exhibits high numerical accuracy but incurs unaffordable memory demand and computational costs, even for middle-scale models, and (2) the spatial-filtering method, which can be implemented easily but results in significant numerical errors under large time steps. However, the intrinsic relationship between these two methods remains unknown. We reveal the intrinsic relationship between these two methods by considering the eigenvalue perturbation method for a homogeneous model. Using this relationship, we derive an analytical spatial filter for the homogeneous model and develop an optimal spatial filter for heterogeneous models. Compared to the classical spatial-filtering method, which removes all high-wavenumber components, our method retains all high-wavenumber components that contribute to wave propagation while eliminating other high-wavenumber components that cause instability. For the same numerical accuracy, the maximum time step allowed by the proposed method is approximately twice that of the classical spatial-filtering method. Compared to the eigenvalue perturbation method, our method can be used in large-scale models without additional memory consumption and computational cost. This can significantly accelerate the pseudospectral method with significantly larger time steps beyond the CFL stability limit, which is particularly promising for large-scale models with a fine grid.

A monolithic, finite element-based strategy for solving fluid structure interaction problems coupled with electrostatics

In this work, we present a monolithic finite-element-based strategy for problems involving the deformation of a hyperelastic solid, an incompressible fluid and electrostatics. We use a two-field hybrid Lagrangian formulation for the structure, and a velocity-based ALE mixed formulation for the fluid with appropriately chosen interpolations for the various field variables to ensure stability of the resulting numerical procedure. The equations of electrostatics are solved on the reference configuration over both the solid and fluid domains, with voltage and electric displacement continuity imposed at the interface. Keeping in view that the thickness of typical MEMS structures is small, a stress-based hybrid formulation is used to prevent locking. The use of a monolithic strategy provides a robust and stable algorithm, and allows the user to take large time steps. The consistent linearization ensures a quadratic rate of convergence of the non-linear iterations at each time step. Detailed expressions for the tangent stiffness and associated matrices for the three-way coupled problem are provided. The robustness and accuracy of the proposed method are demonstrated by solving several benchmark problems from the literature.

Accelerating high order discontinuous Galerkin solvers using neural networks: 3D compressible Navier-Stokes equations

We propose to accelerate a high order discontinuous Galerkin solver using neural networks. We include a corrective forcing to a low polynomial order simulation to enhance its accuracy. The forcing is obtained by training a deep fully connected neural network, using a high polynomial order simulation but only for a short time frame. With this corrective forcing, we can run the low polynomial order simulation faster (with large time steps and low cost per time step) while improving its accuracy.We explored this idea for a 1D Burgers' equation in [1], and we have extended this work to the 3D Navier-Stokes equations, with and without a Large Eddy Simulation closure model. We test the methodology with the turbulent Taylor Green Vortex case and for various Reynolds numbers (30, 200 and 1600). In addition, the Taylor Green Vortex evolves with time and covers laminar, transitional, and turbulent regimes, as time progresses.The proposed methodology proves to be applicable to a variety of flows and regimes. The results show that the corrective forcing is effective in all Reynolds numbers and time frames (excluding the initial flow development). We can train the corrective forcing with a polynomial order of 8, to increase the accuracy of simulations from a polynomial order 3 to 6, when correcting outside the training time frame. The low order corrected solution is 4 to 5 times faster than a simulation with comparable accuracy (polynomial order 6).Additionally, we explore changes in the hyperparameters and use transfer learning to speed up the training. We observe that it is not useful to train a corrective forcing using a different flow condition. However, an already trained corrective forcing can be used to initialise a new training (at the correct flow conditions) to obtain an effective forcing with only a few training iterations.

Analysis and evaluation of steady-state and non-steady-state preserving operator splitting schemes for reaction-diffusion(-advection) problems

Reaction-diffusion(-advection) problems are well-known in chemical engineering and computational fluid dynamics. A common feature of these systems is a linear (or non-stiff) transport sub-system and a non-linear (or stiff), highly coupled chemistry sub-system. The expected numerical effort often prohibits a fully coupled solution of the system. Therefore, the system is split into a transport and a reaction sub-system, and each sub-system is solved using specialized solvers. Operator splitting schemes are required to reconstruct the solution of the initial system from the sub-system solutions. Steady-state preserving splitting schemes are particularly essential for steady-state calculations since local time stepping (LTS) or other methods based on fictional time rely on large time steps to be efficient. This work formally analyzes common splitting schemes for reaction-diffusion problems by stability analysis and checking the steady-state preservation of a representative linear scalar problem. Balanced and Simpler Balanced splitting are the only steady-state conservative schemes analyzed. Three new steady-state preserving splitting schemes are proposed based on the findings of the formal analysis. To achieve steady-state preservation, the new schemes use splitting constants based on either the mixing derivative or the chemistry derivative. The formal analysis is accompanied by a dimensionless perfectly stirred reactor (PSR) case and a hydrogen combustion case, both known to be challenging for operator splitting schemes. The test case results are in line with the theoretical results but indicate that the scalar linear analysis is nonviable to capture the full effects of the chemical sub-system. The Simpler and the newly proposed Consistent Staggered splitting schemes give significantly better results than the remaining ones while being second-order and first-order accurate, respectively. If temporal accuracy is irrelevant, e.g., for steady-state solvers, the proposed Consistent Adaptive splitting scheme is promising since it preserves the steady-state solution first-order time accurate with less function evaluations.

An unconditionally energy-stable scheme for the convective heat transfer equation

PurposeThis paper aims to present an unconditionally energy-stable scheme for approximating the convective heat transfer equation.Design/methodology/approachThe scheme stems from the generalized positive auxiliary variable (gPAV) idea and exploits a special treatment for the convection term. The original convection term is replaced by its linear approximation plus a correction term, which is under the control of an auxiliary variable. The scheme entails the computation of two temperature fields within each time step, and the linear algebraic system resulting from the discretization involves a coefficient matrix that is updated periodically. This auxiliary variable is given by a well-defined explicit formula that guarantees the positivity of its computed value.FindingsCompared with the semi-implicit scheme and the gPAV-based scheme without the treatment on the convection term, the current scheme can provide an expanded accuracy range and achieve more accurate simulations at large (or fairly large) time step sizes. Extensive numerical experiments have been presented to demonstrate the accuracy and stability performance of the scheme developed herein.Originality/valueThis study shows the unconditional discrete energy stability property of the current scheme, irrespective of the time step sizes.

An energy stable incompressible SPH method with consistent solid boundary treatment

We propose two novel solid wall boundary treatments on the energy-stable Smoothed Particle Hydrodynamics (SPH) method for the Navier–Stokes system which models single-phase incompressible flows. In this work, the solid wall is discretized with dummy particles, and an external force is introduced to restrict the motion of solid particles, whose physical significance and numerical necessity are analyzed. We construct consistent fluid–solid coupling schemes based on the approaches so that the boundary condition is involved in the same equations system with particular constraints. Furthermore, the schemes are unconditionally energy-stable, which allows large time step sizes and leads to efficient computation. A variety of examples are demonstrated to verify the corresponding numerical analysis result, which shows that the approaches can easily handle no-slip and moving boundaries. Finally, we highlight that the energy-stable ISPH with consistent boundary treatments greatly benefits solving long-term problems with steady states.

Computing enhancement of the nonlinear SPN approximations of radiative heat transfer in participating material

Anisotropic mesh adaptation is an efficient procedure for controlling the output error of finite element simulations, particularly when used for three-dimensional problems. In this paper, we present an enhanced computational algorithm based on an anisotropic mesh adaptation for nonlinear SPN approximations of radiative heat transfer in both two- and three-dimensional enclosures. Using an asymptotic analysis for the optical scale in the radiative transfer, the integro-differential equation is replaced by a series of partial differential equations of elliptic type. The nonlinear coupling between the heat transfer and radiation in the SPN equations does not depend on the direction ordinates. In an anisotropic participating media, internal boundary layers with different magnitudes occur for each direction and developing an efficient numerical algorithm to accurately resolve them is a challenging task. In the present study, we propose an adaptive finite element method using an efficient hierarchical error estimator. A second-order scheme is used for the time integration and a Newton-type solver is implemented for the fully coupled nonlinear system. The proposed method has the potential to use large timesteps in the simulations for radiative heat transfer in anisotropic media at high temperatures. To demonstrate the viability of the method, several examples are presented for two- and three-dimensional problems. The numerical results confirm the capability of the proposed method to efficiently solve the nonlinear SPN approximations of radiative transfer in anisotropic media.

High order asymptotic preserving finite difference WENO schemes with constrained transport for MHD equations in all sonic Mach numbers

In this paper, a high-order semi-implicit (SI) asymptotic preserving (AP) and divergence-free finite difference weighted essentially nonoscillatory (WENO) scheme is proposed for magnetohydrodynamic (MHD) equations. We consider the sonic Mach number ε ranging from 0 to O(1). High-order accuracy in time is obtained by SI implicit-explicit Runge–Kutta (IMEX-RK) time discretization. High-order accuracy in space is achieved by finite difference WENO schemes with characteristic-wise reconstructions. A constrained transport method is applied to maintain a discrete divergence-free condition. We formally prove that the scheme is AP. Asymptotic accuracy (AA) in the incompressible MHD limit is obtained if the implicit part of the SI IMEX-RK scheme is stiffly accurate. Numerical experiments are provided to validate the AP, AA, and divergence-free properties of our proposed approach. Besides, the scheme can well capture discontinuities such as shocks in an essentially non-oscillatory fashion in the compressible regime, while it is also a good incompressible solver with uniform large-time step conditions in the low sonic Mach limit.

A high efficiency Lie derivative algorithm for the nonautonomous nonlinear systems

Numerical procedure plays a key role in tackling the solutions of nonlinear dynamical systems. With the advent of the age of big data and high-power computing, developing efficient and fast numerical algorithms is an urgent task. This paper extends the Lie derivative discretization algorithm to the nonautonomous nonlinear systems and investigates the numerical solutions of the systems. The periodic solutions of three different classical nonlinear systems are calculated, and the results are compared to those values calculated from the Runge–Kutta fourth-order algorithm, which demonstrated that the Lie derivative algorithm has the advantages of large time step and short computation time.

A modified and efficient phase field model for the biological transport network

This paper aims to establish the biological transport network based on the phase field model. In order to ensure that the topological shape is formed under the guidance of electrical conductivity, we generate the biological networks with sufficient information based on the distributions of the venation of the leaf represented by the reaction-diffusion model. We modify the original energy of the network generating model by considering the auxin gradient property. By applying the gradient flow method to minimize the modified energy, we derive the Poisson type equation for pressure, the reaction-diffusion type equation for the network conductance, and the Allen-Cahn type equation for the phase field. The proposed model is significant on the investigation of phase transitions by considering the gradient properties on the boundaries. We have innovatively added conductivity and phase-field coupling terms that inhibit perpendicular transport of nutrients, making it easy to generate thin branches from the trunk through this model. In order to obtain the second-order temporal accuracy, we take the Crank–Nicolson method for the governing system. To obtain the second-order spatial accuracy, we discretize the coupling system with the central finite difference method and linearize the nonlinear terms semi-explicitly to form a linear system at each time step. The discrete energy dissipation is provably preserved and we can use a larger time step. We apply the preconditioned conjugate gradient method with the multigrid method as a preconditioner to implement a practical algorithm with only linear algebraic complexity. The proposed algorithm is easy to implement and achieves a fast convergence. Various numerical tests are demonstrated to verify the efficiency, stability, and robustness of the proposed method.

On the Application of Non-Gaussian Noise in Stochastic Langevin Simulations

In light of recent advances in time-step independent stochastic integrators for Langevin equations, we revisit the considerations for using non-Gaussian distributions for the thermal noise term in discrete-time thermostats. We find that the desirable time-step invariance of the modern methods is rooted in the Gaussian noise, and that deviations from this distribution will distort the Boltzmann statistics arising from the fluctuation-dissipation balance of the integrators. We use the GJ stochastic Verlet methods as the focus of our investigation since these methods are the ones that contain the most accurate thermodynamic measures of existing methods. Within this set of methods we find that any distribution of applied noise, which satisfies the two first moments given by the fluctuation-dissipation theorem, will result in correct, time-step independent results that are generated by the first two moments of the system coordinates. However, if non-Gaussian noise is applied, undesired deviations in higher moments of the system coordinates will appear to the detriment of several important thermodynamic measures that depend especially on the fourth moments. The deviations, induced by non-Gaussian noise, become significant with the one-time-step velocity attenuation, thereby inhibiting the benefits of the methods. Thus, we conclude that the application of Gaussian noise is necessary for reliable thermodynamic results when using modern stochastic thermostats with large time steps.

A Practical Approach to Wave Function Propagation, Hopping Probabilities, and Time Steps in Surface Hopping Calculations.

We compare several established approaches for propagating wave functions and calculating hopping probabilities within the fewest switches surface hopping (FSSH) algorithm for difficult cases with many electronic states and many trivial crossings. If only a single time step (Δtc) is employed, we find that no published approach can accurately capture the dynamics correctly unless Δtc → 0 (which is not computationally feasible). If multiple time steps are employed, for a fixed classical time step (Δtc), a robust scheme can be found for dynamically choosing quantum time steps (δtq1 and δtq2) and calculating hopping probabilities so that one can systematically reduce all errors and achieve maximally efficient accuracy; scattering calculations confirm that one can choose a fairly large classical time step. The robust scheme presented here uses both the "local diabatic" and adiabatic interpolation and thus borrows elements from both the Granucci/Persico and Meek/Levine algorithms. Our findings should be broadly applicable in the future.

A Novel Energy and Momentum Consistent Mixed Formulation for Coupled Nonlinear Electro‐Thermo‐Elastodynamics

AbstractIn the present contribution, a novel mixed formulation aimed at the energy‐momentum consistent simulation of coupled nonlinear electro‐thermo‐elastodynamic systems, in particular dielectric elastomer actuators, is proposed. It is essentially based on a mixed framework for elastodynamics in the case of polyconvex stored energy functions. In accordance with this framework, the properties of the rediscovered tensor cross product are exploited in a first step to derive a mixed formulation via a Hu‐Washizu type extension of the stored energy function. Afterwards, the corresponding strong form is derived and supplemented with the (initial) boundary value problems of thermodynamics and electrostatics. By additionally choosing an appropriate polyconvexity‐inspired energy density function, this procedure leads to a fully coupled electro‐thermo‐elastodynamic formulation that benefits from the properties of the underlying mixed framework.Furthermore, an energy‐momentum consistent time integration scheme is proposed for the novel framework, where discrete derivatives in the sense of Gonzalez are employed. The formulation is second‐order accurate and stable even for large time step sizes. Eventually, the performance of the novel formulation is illustrated in a numerical example.