Large Time Step Size Research Articles

Dynamically regularized Lagrange multiplier schemes with energy dissipation for the incompressible Navier-Stokes equations

In this paper, we present efficient numerical schemes based on the Lagrange multiplier approach for the Navier-Stokes equations. By introducing a dynamic equation (involving the kinetic energy, the Lagrange multiplier, and a regularization parameter), we form a new system which incorporates the energy evolution process but is still equivalent to the original equations. Such nonlinear system is then discretized in time based on the backward differentiation formulas, resulting in a dynamically regularized Lagrange multiplier (DRLM) method. First- and second-order DRLM schemes are derived and shown to be unconditionally energy stable with respect to the original variables. The proposed schemes require only the solutions of two linear Stokes systems and a scalar quadratic equation at each time step. Moreover, with the introduction of the regularization parameter, the Lagrange multiplier can be uniquely determined from the quadratic equation, even with large time step sizes, without affecting accuracy and stability of the numerical solutions. Fully discrete energy stability is also proved with the Marker-and-Cell (MAC) discretization in space. Various numerical experiments in two and three dimensions verify the convergence and energy dissipation as well as demonstrate the accuracy and robustness of the proposed DRLM schemes.

Exponential time difference methods with spatial exponential approximations for solving boundary layer problems

AbstractIn this work, an efficient and stable exponential time difference method is presented for solving boundary layer problems. By combining exponential time difference schemes with spatial direct discontinuous Galerkin discretization based on exponential boundary layer approximations, the proposed algorithm not only may admit large time step sizes but also could provide good spatial approximations even if on rather coarse spatial grids in the boundary layer. Some energy stabilities of the numerical scheme are rigorously derived. Numerical examples illustrate the accuracy, stability and efficiency of the algorithm.

An Integral-like Numerical Approach for Solving Burgers’ Equation

The Burgers’ equation, commonly appeared in the study of turbulence, fluid dynamics, shock waves, and continuum mechanics, is a crucial part of the dynamical core of any numerical weather model, influencing simulated meteorological phenomena. While past studies have suggested several robust numerical approaches for solving the equation, many are too complicated for practical adaptation and too computationally expensive for operational deployment. This paper introduces an unconventional approach based on spline polynomial interpolations and the Hopf-Cole transformation. Using Taylor expansion to approximate the exponential term in the Hopf-Cole transformation, the analytical solution of the simplified equation is discretized to form our proposed scheme. The scheme is explicit and adaptable for parallel computing, although certain types of boundary conditions need to be employed implicitly. Three distinct test cases were utilized to evaluate its accuracy, parallel scalability, and numerical stability. In the aspect of accuracy, the schemes employed cubic and quintic spline interpolation perform equally well, managing to reduce the <i>ӏ</i><sub>1</sub>, <i>ӏ</i><sub>2</sub>, and <i>ӏ</i><sub>∞</sub> error norms down to the order of 10<sup>−4</sup>. Parallel scalability observed in the weak-scaling experiment depends on time step size but is generally as good as any explicit scheme. The stability condition is <i>ν</i>∆<i>t</i>/∆<i>x</i><sup>2</sup> > 0.02, including the viscosity coefficient <i>ν</i> due to the Hopf-Cole transformation step. From the stability condition, the schemes can run at a large time step size ∆<i>t</i> even when using a small grid spacing ∆<i>x</i>, emphasizing its suitability for practical applications such as numerical weather prediction.

Precise integration boundary element method for solving nonlinear transient heat conduction problems with temperature dependent thermal conductivity and heat capacity

The precise integration boundary element method (PIBEM) is used to solve nonlinear transient heat conduction problems with temperature dependent thermal conductivity and heat capacity in this article. At first, a boundary-domain integral equation can be obtained by using the fundamental solution for steady linear heat conduction problems. Secondly, the domain integrals are converted into boundary integrals by using the radial integration method to get a pure boundary integral equation, which can retain the advantage of dimension reduction of the boundary element method. Then, after discretization, a first-order nonlinear differential equation system in time is obtained, and the precise integration method with predictor-corrector technique is used to solve the system of nonlinear equations. In the end, several numerical examples are presented to verify the accuracy and stability of the presented method. The results obtained by the presented method are less affected by the time step size and satisfactory results can be obtained even for a large time step size.

An analysis of high order FEM and IGA for explicit dynamics: Mass lumping and immersed boundaries

SummaryWe investigate the behavior of different shape functions for the discretization of hyperbolic problems. In particular, we consider classical Lagrange polynomials and B‐splines. The studies focus on the performance of the these functions as a spatial discretization approach combined with an explicit time marching scheme. In this regard, a major concern is the maximum eigenvalue that imposes restrictions on the critical time step size and suitable lumping techniques that yield a diagonal mass matrix. The accuracy of the discretization methods is assessed in an asymptotic manner in terms of the convergence of eigenvalues and eigenvectors. Further, the global accuracy is investigated in terms of the full spectrum. The results show that B‐spline discretization with a consistent mass matrix are more accurate than those based on Lagrange shape functions, which holds true in the boundary‐fitted as well as in the immersed setting. On the other hand, Lagrange shape functions are more robust with respect to standard lumping techniques, which cannot be directly applied for B‐splines without loss of accuracy. In general, we observe that none of the standard lumping schemes yields optimal results for B‐splines, even in the boundary‐fitted setting. For the immersed setting, also Lagrange shape functions show a drop in accuracy which depends on the position of the boundary that cuts the element. Several remedies are considered in order to overcome these issues, including interpolatory B‐spline bases as well as eigenvalue stabilization methods. While accuracy and stability can be improved using these remedies, we conclude from our study that it is still an open question, how to design a discretization method that achieves large critical time step sizes in combination with a diagonal mass matrix and high accuracy in the immersed setting. We note that these considerations primarily relate to linear structural dynamics applications, such as for example, structural acoustics. In nonlinear problems, such as automotive crash dynamics, other considerations predominate. An example of a one‐dimensional elastic‐plastic bar impacting a rigid wall is illustrative.

Simulation of nuclear fuel clad high-temperature ballooning under loss-of-coolant accident conditions considering anisotropic creep

Zirconium (Zr) alloys are standard clad material in commercial light water reactors. However, it is vulnerable to ballooning and burst under high-temperature conditions such as LOCA (Loss-Of-Coolant Accident), due to large internal/external pressure difference. Considering that clad ballooning is detrimental in terms of nuclear reactor safety especially for high-burnup fuel rod, high-fidelity numerical simulation of Zircaloy clad ballooning is desirable. This paper presented a finite-element based implementation for the simulation of Zircaloy clad ballooning. The developed simulation capability predicts progressive clad ballooning up to failure given empirical burst limits. The material models coupled together are phase transformation and high-temperature creep. Especially, anisotropic creep was considered for α phase Zircaloy. The implemented algorithm for stress update is implicit in time integration, which allows relatively large time step size to efficiently simulate the process that can last for hundreds of seconds. An analytic benchmark, the PUZRY ballooning and burst test, and the Hardy stress-rupture test were utilized to verify and validate the simulation capability. The anisotropic creep calculation conforms to the analytical solution. The comparison with PUZRY and Hardy test data shows that considering anisotropic creep in the α phase region is essential. 3D ballooning was also simulated by considering circumferential temperature difference for both conventional and novel annular clad designs. Interestingly, it is the first time to note that the risk of buckling failure is predicted for the double-cooled annular nuclear fuel clad during the accident condition, only after considering the anisotropic creep algorithm.

Higher‐order inverse mass matrices for the explicit transient analysis of heterogeneous solids

AbstractNew methods are presented for the direct computation of higher‐order inverse mass matrices (also called reciprocal mass matrices) that are used for explicit transient finite element analysis. The motivation of this work lies in the need of having appropriate sparse inverse mass matrices, which present the same structure as the consistent mass matrix, preserve the total mass, predict suitable frequency spectrum and dictate sufficiently large critical time step sizes. For an efficient evaluation of the reciprocal mass matrix, the projection matrix should be diagonal. This condition can be satisfied by adopting dual shape functions for the momentum field, generated from the same shape functions used for the displacement field. A theoretically consistent derivation of the inverse mass matrix is based on the three‐field Hamilton principle and requires the projection matrix to be evaluated from the integral of these shape functions. Unfortunately, for higher‐order FE shape functions and serendipity FE elements, the projection matrix is not positive definitive and can not be employed. Therefore, we study several lumping procedures for higher order reciprocal mass matrices considering their effect on total‐mass preserving, frequency spectra and accuracy in explicit transient simulations. The article closes with several numerical examples showing suitability of the direct inverse mass matrix in dynamics.

A collision-based hybrid method for the BGK equation

We apply the collision-based hybrid method introduced by Hauck and McClarren [1] to the Boltzmann equation with the BGK operator and a hyperbolic scaling. An implicit treatment of the source term is used to handle stiffness associated with the BGK operator. Although it helps the numerical scheme become stable with a large time step size, it is still not obvious to achieve the desired order of accuracy due to the relationship between the size of the spatial cell and the mean free path. Without asymptotic preserving property, a very restricted grid size is required to resolve the mean free path, which is not practical. Our approaches are based on the noncollision-collision decomposition of the BGK equation. We introduce the arbitrary order of nodal discontinuous Galerkin (DG) discretization in space with a semi-implicit time-stepping method; we employ the backward Euler time integration for the uncollided equation and the 2nd order predictor-corrector scheme for the collided equation, i.e., both source terms in uncollided and collided equations are treated implicitly and only streaming term in the collided equation is solved explicitly. This improves the computational efficiency without the complexity of the numerical implementation. Numerical results are presented for various Knudsen numbers to present the effectiveness and accuracy of our hybrid method. Also, we compare the solutions of the hybrid and non-hybrid schemes.

An unconditionally energy-stable and orthonormality-preserving iterative scheme for the Kohn-Sham gradient flow based model

We propose an unconditionally energy-stable, orthonormality-preserving, component-wise splitting iterative scheme for the Kohn-Sham gradient flow based model in the electronic structure calculation. We first study the scheme discretized in time but still continuous in space. The component-wise splitting iterative scheme changes one wave function at a time, similar to the Gauss-Seidel iteration for solving a linear equation system. At the time step n, the orthogonality of the wave function being updated to other wave functions is preserved by projecting the gradient of the Kohn-Sham energy onto the subspace orthogonal to all other wave functions known at the current time, while the normalization of this wave function is preserved by projecting the gradient of the Kohn-Sham energy onto the subspace orthogonal to this wave function at tn+1/2. The unconditional energy stability is nontrivial, and it comes from a subtle treatment of the two-electron integral as well as a consistent treatment of the two projections. Rigorous mathematical derivations are presented to show our proposed scheme indeed satisfies the desired properties. We then study the fully-discretized scheme, where the space is further approximated by a conforming finite element subspace. For the fully-discretized scheme, not only the preservation of orthogonality and normalization (together we called orthonormalization) can be quickly shown using the same idea as for the semi-discretized scheme, but also the highlight property of the scheme, i.e., the unconditional energy stability can be rigorously proven. The scheme allows us to use large time step sizes and deal with small systems involving only a single wave function during each iteration step. Several numerical experiments are performed to verify the theoretical analysis, where the number of iterations is indeed greatly reduced as compared to similar examples solved by the Kohn-Sham gradient flow based model in the literature.

Second-order efficient algorithm for coupled nonlinear model of groundwater transport system

In this paper, an efficient modified two-grid block-centered finite difference algorithm with second-order backward difference formula is presented for coupled nonlinear reaction-diffusion equations of the solutes in the groundwater flow. The main steps of the modified two-grid algorithm include dealing with a small-scale nonlinear problem on the coarse grid with a large time step size and then dealing with a large-scale linear problem on the fine grid with a smaller time step size base on the second-order backward difference formula. Moreover, the theoretical analysis illustrates that the modified two-grid algorithm is able to maintain second convergence order in space-time domain. Several numerical experiments including a practical example are carried out to support the theoretical results and verify that the modified two-grid block-centered finite difference with second-order backward difference formula algorithm can achieve high efficiency while maintaining the second-order accuracy, meanwhile, has a certain applicable value in actual problem.

High-order composite implicit time integration schemes based on rational approximations for elastodynamics

A new approach for developing implicit composite time integration schemes is established starting with rational approximations of the matrix exponential in the solution of the equations of motion. The rational approximations are designed to have the same effective stiffness matrix in all sub-steps. An efficient algorithm is devised so that the implicit equation for a sub-step is in the same form as that of the trapezoidal method. The proposed M-schemes are Mth order accurate using M sub-steps. The amount of numerical dissipation is controlled by a user-specified parameter ρ∞. The (M+1)-schemes that are (M+1)th order accurate using M sub-steps can also be constructed, but the amount of numerical dissipation is built in the schemes and not adjustable. The order of accuracy of all the schemes is not affected by external forces and physical damping. Numerical examples demonstrate that the proposed high-order composite schemes are effective for introducing numerical dissipation and allow the use of large time step sizes in wave propagation problems. The source code written in MATLAB is available for download at: https://github.com/ChongminSong/CompsiteTimeIntegration.

Development of a Composite Implicit Time Integration Scheme for Three-Dimensional Discontinuous Deformation Analysis

Discontinuous deformation analysis (DDA) is a discontinuum-based and implicit method for investigating the deformational behavior of block systems. The constant acceleration integration (CAI) scheme characterized by unconditional stability is employed in the traditional DDA. In this study, the problems of the CAI scheme regarding the time step in DDA are pointed out. A too large or too small time step size adopted in the CAI scheme will have adverse effects on the DDA computation. To overcome the weaknesses, an alternative composite implicit time integration (CITI) scheme, which is a combination of the trapezoidal rule and the three-point backward Euler method, is implemented in the three-dimensional (3D) DDA method. Verification examples and slope numerical simulations are presented to illustrate the accuracy and effectiveness of the proposed methodology. The results showed that the CITI scheme can overcome the numerical error caused by the large time step size, and the algorithm damping is closely related to the time step size and the selected splitting parameter.

A semi-implicit semi-Lagrangian method for simulating immersed boundary motion under high inertia and elasticity

In this paper, we present an efficient and stable fractional-step 2D immersed boundary (IB) method for solving the interaction problems between bulk fluid and elastic interface, in particular, when the fluid inertia and the interfacial elasticity are the significant factors affecting its dynamics. In myriads of real-world applications, the effects of high inertia and elasticity are dominant. So the complex fluid dynamics under such harsh conditions is an important topic in computational physics and is inherently challenging due to high computational complexity. In turn, it requires to solve the governing equations of elastic interfacial motion in an implicit manner so that more stable simulations can be performed by relaxing the Courant-Friedrichs-Lewy (CFL) condition. The contributions of the proposed approach are three folds. First, an iteration-free semi-Lagrangian method is employed in Navier-Stokes (NS) equations. Second, the elastic force acting along the interface is treated semi-implicitly in IB formulations. Both approaches improve the numerical stability associated with the high fluidic inertia and interfacial elasticity. Finally, to solve the resulting linear system, our novel idea is to transform the original 3-by-3 block matrix system into a reduced 2-by-2 block matrix system using a discrete projection operator in a staggered grid, and then explicitly represent the exact solution via the Schur complement of the Helmholtz operator. Owing to this feature, we refer to this proposed approach as reduced immersed boundary method (rIBM). We show that the two systems are equivalent in theory, whereas the conventional immersed boundary projection method (IBPM) modifies the discrete momentum equation in the original system. A series of numerical tests is conducted to confirm the stability of the rIBM using relatively larger time-step sizes, specifically with Reynolds number and inverse capillary number equal to or larger than approximately 1000. By estimating the computational time, the numerical efficiency of the proposed method is further verified in comparison with the conventional IBPM and the Crank-Nicolson scheme-based IB method. In conclusion, the proposed approach not only improves the numerical stability, but also increases the computational speed, suitable for solving more realistic problems.

Scale-Resolving Hybrid RANS-LES Simulation of a Model Kaplan Turbine on a 400-Million-Element Mesh

Double-regulated Kaplan turbines with adjustable guide vanes and runner blades offer a high degree of flexibility and good efficiency for a wide range of operating points. However, this also leads to a complex geometry and flow guidance with, for example, vortices of different sizes and strengths. The flow in a draft tube is especially challenging to simulate mainly due to flow phenomena, like swirl, separation and strong adverse pressure gradients, and a strong dependency on the upstream flow conditions. Standard simulation approaches with RANS turbulence models, a coarse mesh and large time step size often fail to correctly predict performance and can even lead to wrong tendencies in the overall behavior. To reveal occurring flow phenomena and physical effects, a scale-resolving hybrid RANS-LES simulation on a block structured mesh of about 400 million hexahedral elements of a double-regulated five-blade model Kaplan turbine is carried out. In this paper, first, the results of the ongoing simulation are presented. The major part of the simulation domain is running in LES mode and seems to be properly resolved. The validation of the simulation results with the experimental data shows mean deviations of less than 0.8% in the global results, i.e., total head and power, and a good visual agreement with the three-dimensional PIV measurements of the velocity in the cone and both diffuser channels of the draft tube. In particular, the trend of total head and the results for the draft tube differ significantly between the scale-resolving simulation and a standard RANS simulation. The standard RANS simulation exhibits a highly unsteady behavior of flow, which is not observed in the experiments or scale-resolving simulation.

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.

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.

CDFIB: A convolutional direct forcing immersed boundary method for solving incompressible flows with time-varying geometries

When incorporating the immersed boundary method to model solid obstacles in incompressible flows, the conventional direct forcing immersed boundary (DFIB) method may face the issue of spurious force oscillations (SFO) due to the mass imbalance when integrating drag/lift forces around the solid obstacle. Although the SFO can be controlled with large time step size when employing conventional DFIB method, the method remains prone to large SFO when the time step size is small. When we deal with flow-structure interaction (FSI) scenarios, the solid bodies are driven by hydrodynamic forces that are produced by flow fluids. Unphysical acceleration/deceleration may happen during the simulations due to SFO, thus making the simulations unstable or even diverge. To alleviate the SFO, we thus propose a novel method for modeling the solid bodies in incompressible fluid flows by performing the convolution kernel, dubbed as cDFIB method. cDFIB method inherits the merits of DFIB method, while effectively suppressing SFO for both small and large time step size. From the numerical results of the SFO oscillating cylinder problem, the SFO is proportional to O(Δx) and O(Δt) for present cDFIB method. This property allows cDFIB method to solve challenging problems where employing large time step size becomes infeasible.

Hybrid asynchronous isogeometric Perfectly Matched Layer for transient elastodynamics

This paper presents an isogeometric formulation (IGA) of the Perfectly Matched Layer (PML) for two-dimensional time-domain elastodynamics in unbounded media. The displacement-based unsplit-field PML formulation is first discretized in space using B-spline basis functions of increasing polynomial order p, from 2 to 10, able to deal with the strong solution gradients in the thickness of the PML thanks to the high order continuity provided by the B-spline basis functions. In order to keep the high accuracy of IGA, implicit time integration is chosen in the PML subdomains, whereas the physical domain is integrated with an explicit time integration scheme. The domain decomposition is based on a dual formulation and Heterogeneous Asynchronous Time Integrator (HATI) approach. The coupling between non-conforming IGA-patches for the PML subdomains and the physical domain is carried out using a standard mortar method. It enables us to select for the IGA-PML subdomains both the appropriate space discretization and the appropriate time discretization with an implicit time integrator and a larger time step size. The hybrid multi-time step IGA-PML performance in terms of computation time, is assessed for two-dimensional elastodynamic problems in single and multi-layer media depending on the chosen discretization parameters in space and time.

Carleman linearization approach for chemical kinetics integration toward quantum computation

The Harrow, Hassidim, Lloyd (HHL) algorithm, known as the pioneering algorithm for solving linear equations in quantum computers, is expected to accelerate solving large-scale linear ordinary differential equations (ODEs). To efficiently combine classical and quantum computers for high-cost chemical problems, non-linear ODEs (e.g., chemical reactions) must be linearized to the highest possible accuracy. However, the linearization approach has not been fully established yet. In this study, Carleman linearization was examined to transform nonlinear first-order ODEs of chemical reactions into linear ODEs. Although this linearization theoretically requires the generation of an infinite matrix, the original nonlinear equations can be reconstructed. For the practical use, the linearized system should be truncated with finite size and the extent of the truncation determines analysis precision. Matrix should be sufficiently large so that the precision is satisfied because quantum computers can treat such huge matrix. Our method was applied to a one-variable nonlinear dot{y}=-{y}^{2} system to investigate the effect of truncation orders and time step sizes on the computational error. Subsequently, two zero-dimensional homogeneous ignition problems for H2–air and CH4–air gas mixtures were solved. The results revealed that the proposed method could accurately reproduce reference data. Furthermore, an increase in the truncation order improved accuracy with large time-step sizes. Thus, our approach can provide accurate numerical simulations rapidly for complex combustion systems.