Semi-implicit Time Integration Scheme Research Articles

A Smoothing Scheme for Seismic Wave Propagation Simulation with a Mixed FEM of Diffusionized Wave Equation

In this paper, we propose a smoothing scheme for seismic wave propagation simulation. The proposed scheme is based on a diffusionized wave equation with the fourth-order spatial derivative term. So, the solution requires higher regularity in the usual weak formulation. Reducing the diffusionized wave equation to a coupled system of diffusion equations yields a mixed FEM to ease the regularity. We mathematically explain how our scheme works for smoothing. We construct a semi-implicit time integration scheme and apply it to the wave equation. This numerical experiment reveals that our scheme is effective for filtering short wavelength components in seismic wave propagation simulation.

Semi-Implicit Time Integration with Hessian Eigenvalue Corrections for a Larger Time Step in Molecular Dynamics Simulations.

In molecular dynamics simulations, the limited time step size has been a barrier to simulating long-time behaviors. Implicit time integration methods allow markedly larger time steps than the standard explicit time method, although they have major drawbacks such as overheads solving linear systems and instability of Newton iterations. To overcome these issues, we propose a semi-implicit time integration scheme, the semi-implicit Hessian correction (SimHec) scheme, for overdamped Langevin dynamics. The method focuses on the Hessian matrices of bonded and nonbonded interactions, where components with large negative Hessian eigenvalues are cut off in the linear approximation of momentum equations to avoid instability. The narrow band Hessian matrix enables an efficient parallelized linear solution with an overlapping approximation. We tested SimHec for the interdomain fluctuations in adenylate kinase and the powerstroke transition of myosin II using a coarse-grained protein model. SimHec reproduced the same dynamics as the explicit method, although the transition dynamics tended to be accelerated and fluctuations in bonded potentials were slightly reduced. These deviations were corrected using a hybrid method, SimHec-H, which adds explicit time steps after the semi-implicit time step. The proposed scheme allowed us to use time steps 50-200 times larger than those in explicit time integration, which resulted in a speedup factor of 7-30 taking the overhead into account.

Discontinuous Galerkin discretization for two-equation turbulence closure model

Accurate representation of vertical turbulent fluxes is crucial for numerical ocean modeling, both in global and coastal applications. The state-of-the-art approach is to use two-equation turbulence closure models which introduces two dynamic equations to the system. Solving these equations numerically, however, is challenging due to the strict requirement of positivity of the turbulent quantities (e.g., turbulence kinetic energy and its dissipation rate), and the non-linear source terms that may render the numerical system unstable. In this paper, we present a Discontinuous Galerkin (DG) finite element discretization of the Generic Length Scale (GLS) equations designed to be incorporated in a DG coastal ocean model, Thetis. To ensure numerical stability, the function space for turbulent quantities must be chosen carefully. In this work, we propose to use zeroth degree elements for the turbulent quantities and linear discontinuous elements for the tracers and velocity. The spatial discretization is completed with a positivity preserving semi-implicit time integration scheme. We validate the implementation with standard turbulence closure model benchmarks and an idealized estuary simulation. Finally, we use the full three-dimensional model to simulate the Columbia River plume. The results confirm that the coupled model generates realistic vertical mixing, and remains stable under strongly stratified conditions and strong tidal forcing. River plume characteristics are well captured.

A non-isothermal thermodynamically consistent phase field model for damage, fracture and fatigue evolutions in elasto-plastic materials

This paper presents a general thermodynamically consistent non-isothermal phase field framework to model effects of damage, fracture and fatigue evolutions in elasto-plastic materials under the hypothesis of small strains. The proposed methodology is obtained from the principle of virtual power, energy balance and second law of thermodynamics in the form of a generalized Clausius–Duhem inequality for entropy. Aspects of the energy degradation functions are thoroughly investigated for an isotropic elasto-plastic material with viscous dissipation and constant specific heat. A new combined degradation function is proposed to degrade both elastic and plastic energy densities as an alternative to the classical degradation functions. Our conceptual framework leads to thermodynamically consistent models that may include contributions not usually considered in the literature, as temperature and inertia effects as well as time-rate dependent processes. The governing nonlinear transient equations obtained are solved adopting a semi-implicit time integration scheme combined with the classical Newton–Raphson iterative procedure. Results for an I-shaped specimen made of 7075-T7351 aluminum alloy under different conditions are presented. The proposed model is able to reproduce qualitative and quantitatively the ductile fracture and the fatigue phenomena. In particular for fatigue, both the S–N experimental data and the Paris crack growth curve over cycles are recovered. In addition, the cycle jump strategy reduces four times the CPU time for fatigue simulations.

A finite-element coarse-grid projection method for incompressible flow simulations

Coarse grid projection (CGP) methodology is a novel multigrid method for systems involving decoupled nonlinear evolution equations and linear elliptic Poisson equations. The nonlinear equations are solved on a fine grid and the linear equations are solved on a corresponding coarsened grid. Mapping operators execute data transfer between the grids. The CGP framework is constructed upon spatial and temporal discretization schemes. This framework has been established for finite volume/difference discretizations as well as explicit time integration methods. In this article we present for the first time a version of CGP for finite element discretizations, which uses a semi-implicit time integration scheme. The mapping functions correspond to the finite-element shape functions. With the novel data structure introduced, the mapping computational cost becomes insignificant. We apply CGP to pressure-correction schemes used for the incompressible Navier-Stokes flow computations. This version is validated on standard test cases with realistic boundary conditions using unstructured triangular meshes. We also pioneer investigations of the effects of CGP on the accuracy of the pressure field. It is found that although CGP reduces the pressure field accuracy, it preserves the accuracy of the pressure gradient and thus the velocity field, while achieving speedup factors ranging from approximately 2 to 30. The minimum speedup occurs for velocity Dirichlet boundary conditions, while the maximum speedup occurs for open boundary conditions.

A Parallel Finite Element Method for 3D Two-Phase Moving Contact Line Problems in Complex Domains

Moving contact line problem plays an important role in fluid-fluid interface motion on solid surfaces. The problem can be described by a phase-field model consisting of the coupled Cahn–Hilliard and Navier–Stokes equations with the generalized Navier boundary condition (GNBC). Accurate simulation of the interface and contact line motion requires very fine meshes, and the computation in 3D is even more challenging. Thus, the use of high performance computers and scalable parallel algorithms are indispensable. In this paper, we generalize the GNBC to surfaces with complex geometry and introduce a finite element method on unstructured 3D meshes with a semi-implicit time integration scheme. A highly parallel solution strategy using different solvers for different components of the discretization is presented. More precisely, we apply a restricted additive Schwarz preconditioned GMRES method to solve the systems arising from implicit discretization of the Cahn–Hilliard equation and the velocity equation, and an algebraic multigrid preconditioned CG method to solve the pressure Poisson system. Numerical experiments show that the strategy is efficient and scalable for 3D problems with complex geometry and on a supercomputer with a large number of processors.

Implementation of a Semi-implicit Time Integration Scheme in Non-Hydrostatic Euler Equations

A semi-implicit time integration scheme is implemented in a non-hydrostatic Euler problem on the cubed-sphere grid. The semi-implicit time integration scheme is a 3-stage additive RungeKutta method, which is an Implicit-Explicit (IMEX) multi-stage single-step scheme. The system treats acoustic and gravity waves implicitly and advection explicitly. In the implicit part, we compute the linear system by defining a proper linear operator. Then, the numerical results are presented to compare the semi-implicit time integration scheme with the explicit RungeKutta time integration scheme in non-hydrostatic Euler equations. In terms of accuracy, we demonstrate that the proposed method performs better than other schemes.

Constitutive relations and their time integration for anisotropic elasto-plastic porous materials

Finite deformation constitutive relations are developed for a class of plastically anisotropic porous solids with an underlying evolving microstructure. They are based on a model obtained by homogenization for rigid-perfectly plastic materials containing non-spherical voids. To facilitate numerical implementation, heuristic extensions are proposed to incorporate weak elasticity, strain hardening and accurate void shape evolution. A semi-implicit time integration scheme is used along with the Newton–Raphson method to solve the system of equations resulting from the discretization of the constitutive equations. The procedure to calculate the consistent tangent matrix, which is needed to solve the global force–displacement matrix equation, is summarized. The framework is used to illustrate the predictive capabilities of the model, first under conditions previously assessed against finite element cell model calculations, then under conditions heretofore not examined. The latter include situations of initial anisotropy as well as situations involving significant void distortions, not only in terms of void enlargement or shape change, but also in terms of void rotations. In particular, various combinations of stress triaxiality, initial void shape, void orientation, matrix orthotropy properties and loading directions are simulated. In addition, the finite element implementation of the model is addressed and illustrated for simple cases.

Accurately simulating anisotropic thermal conduction on a moving mesh

We present a novel implementation of an extremum preserving anisotropic diffusion solver for thermal conduction on the unstructured moving Voronoi mesh of the Arepo code. The method relies on splitting the one-sided facet fluxes into normal and oblique components, with the oblique fluxes being limited such that the total flux is both locally conservative and extremum preserving. The approach makes use of harmonic averaging points and a simple, robust interpolation scheme that works well for strong heterogeneous and anisotropic diffusion problems. Moreover, the required discretization stencil is small. Efficient fully implicit and semi-implicit time integration schemes are also implemented. We perform several numerical tests that evaluate the stability and accuracy of the scheme, including applications such as point explosions with heat conduction and calculations of convective instabilities in conducting plasmas. The new implementation is suitable for studying important astrophysical phenomena, such as the conductive heat transport in galaxy clusters, the evolution of supernova remnants, or the distribution of heat from black hole-driven jets into the intracluster medium.

A Semi-Implicit Version of the MPAS-Atmosphere Dynamical Core

Abstract An important question for atmospheric modeling is the viability of semi-implicit time integration schemes on massively parallel computing architectures. Semi-implicit schemes can provide increased stability and accuracy. However, they require the solution of an elliptic problem at each time step, creating concerns about their parallel efficiency and scalability. Here, a semi-implicit (SI) version of the Model for Prediction Across Scales (MPAS) is developed and compared with the original model version, which uses a split Runge–Kutta (SRK3) time integration scheme. The SI scheme is based on a quasi-Newton iteration toward a Crank–Nicolson scheme. Each Newton iteration requires the solution of a Helmholtz problem; here, the Helmholtz problem is derived, and its solution using a geometric multigrid method is described. On two standard test cases, a midlatitude baroclinic wave and a small-planet nonhydrostatic gravity wave, the SI and SRK3 versions produce almost identical results. On the baroclinic wave test, the SI version can use somewhat larger time steps (about 60%) than the SRK3 version before losing stability. The SI version costs 10%–20% more per step than the SRK3 version, and the weak and strong scalability characteristics of the two versions are very similar for the processor configurations the authors have been able to test (up to 1920 processors). Because of the spatial discretization of the pressure gradient in the lowest model layer, the SI version becomes unstable in the presence of realistic orography. Some further work will be needed to demonstrate the viability of the SI scheme in this case.

Highly Stable Explicit Temporal Integration for Discrete Element Computations

AbstractThis paper presents the adaptation of a semi-implicit time integration scheme that has been reported in the literature and its implementation for use in discrete element methods (DEM). The computational efficiency of DEM methods is primarily associated with the selection of time increment sizes, as dictated by stability requirements, and other factors that depend on the particulars of the DEM implementation and the problems solved. The proposed time integration scheme and the associated DEM are developed for problems pertaining to rigid-particle interaction and interaction of elastic bodies that are modeled as a cluster of rigid interconnected particles. Verification studies that consider nonlinear problems demonstrate that the proposed algorithm is unconditionally stable and accurate even for large time step sizes and when applicable does not require any inversions of the system matrices. Assessment studies on the accuracy, stability, and computational efficiency of the method have been conducted...

Finite element method implementation of a nonlocal viscoplastic model with consideration of temperature dependence

This article presents development of a damage-coupled viscoplastic constitutive model with temperature consideration. The model is subjected to an explicit nonlocal treatment within the characteristic length that is not limited to one local element as conventional ones. Temperature dependence of material behavior is incorporated into the model to account for the material property degradation at elevated temperature. A test program to determine the correlation between material parameters and temperature is also presented. The nonlocal damage-coupled viscoplastic material model is implemented in a commercial finite element program ABAQUS through its user-defined material subroutine UMAT using a semi-implicit time integration scheme. The model is applied to predict the behavior of 63Sn37Pb soldered structure. The mesh sensitivity of the model is discussed, as well as the efficiency of deduced consistent tangent modulus.

Interactive real-time physics: An intuitive approach to form-finding and structural analysis for design and education

Real-time physics simulation has been extensively used in computer games, but its potential has yet to be fully realized in design and education. We present an interactive 3D physics engine with a wide variety of applications.In common with traditional FEM, the use of a local element stiffness matrix is retained. However, unlike typical non-linear FEM routines elements forces, moments and inertia are appropriately lumped at nodes following the dynamic relaxation method. A semi-implicit time integration scheme updates linear and angular momentum, and subsequently the local coordinate frames of the nodes. A co-rotational approach is used to compute the resultant field of displacements in global coordinates including the effect of large deformations. The results obtained compare well against established commercial software.We demonstrate that the method presented allows the making of interactive structural models that can be used in teaching to develop an intuitive understanding of structural behaviour. We also show that the same interactive physics framework allows real-time optimization that can be used for geometric and structural design applications.

An efficient fully-discrete local discontinuous Galerkin method for the Cahn–Hilliard–Hele–Shaw system

In this paper, we develop an efficient and energy stable fully-discrete local discontinuous Galerkin (LDG) method for the Cahn–Hilliard–Hele–Shaw (CHHS) system. The semi-discrete energy stability of the LDG method is proved firstly. Due to the strict time step restriction (Δt=O(Δx4)) of explicit time discretization methods for stability, we introduce a semi-implicit time integration scheme which is based on a convex splitting of the discrete Cahn–Hilliard energy. The unconditional energy stability has been proved for this fully-discrete LDG scheme. The fully-discrete equations at the implicit time level are nonlinear. Thus, the nonlinear Full Approximation Scheme (FAS) multigrid method has been applied to solve this system of algebraic equations, which has been shown the nearly optimal complexity numerically. Numerical results are also given to illustrate the accuracy and capability of the LDG method coupled with the multigrid solver.

PPBLAB: A New Multivariate Population Balance Environment for Particulate System Modelling and Simulation

Abstract This work presents a new windows-based MATLAB program, which is called PPBLAB (Particulate Population Balance Laboratory) for modelling and numerical simulation of particulate systems. As a first step, liquid-liquid extraction columns are modelled using the population balance equation as a mathematical framework. Up to date population balance models and solvers are incorporated. The discretization of the spatial domain is based on the finite volume method with flux vector splitting. A strongly stable semi-implicit first order time integration scheme is used to resolve such a large and stiff ODE system. The MATLAB GUI is used to make PPBLAB a user friendly program, which allows the user to define and simulate liquid extraction columns. A thermodynamics package TEA-COCO, which obtained from CAPE-OPEN, is linked to PPBLAB. Therefore, a special interface is designed for the purpose of data exchange between PPBLAB and CAPE-OPEN TEA tool. The solute distribution coefficient in ternary systems is predicted using the UNIQUAQ model with a special optimization tool to estimate the binary interaction parameters based on infinite conditions. A pilot plant Kuhni extraction column is simulated and tested. Full analysis and performance of the PPBLAB are carried out and validated against experimental data.

High-order accurate solution of the incompressible Navier–Stokes equations on massively parallel computers

The emergence of “petascale” supercomputers requires us to develop today’s simulation codes for (incompressible) flows by codes which are using numerical schemes and methods that are better able to exploit the offered computational power. In that spirit, we present a massively parallel high-order Navier–Stokes solver for large incompressible flow problems in three dimensions. The governing equations are discretized with finite differences in space and a semi-implicit time integration scheme. This discretization leads to a large linear system of equations which is solved with a cascade of iterative solvers. The iterative solver for the pressure uses a highly efficient commutation-based preconditioner which is robust with respect to grid stretching. The efficiency of the implementation is further enhanced by carefully setting the (adaptive) termination criteria for the different iterative solvers. The computational work is distributed to different processing units by a geometric data decomposition in all three dimensions. This decomposition scheme ensures a low communication overhead and excellent scaling capabilities. The discretization is thoroughly validated. First, we verify the convergence orders of the spatial and temporal discretizations for a forced channel flow. Second, we analyze the iterative solution technique by investigating the absolute accuracy of the implementation with respect to the different termination criteria. Third, Orr–Sommerfeld and Squire eigenmodes for plane Poiseuille flow are simulated and compared to analytical results. Fourth, the practical applicability of the implementation is tested for transitional and turbulent channel flow. The results are compared to solutions from a pseudospectral solver. Subsequently, the performance of the commutation-based preconditioner for the pressure iteration is demonstrated. Finally, the excellent parallel scalability of the proposed method is demonstrated with a weak and a strong scaling test on up to O ( 10 4 ) processing units and O ( 10 11 ) grid points.

On finite strain micromorphic elastoplasticity

In the micromorphic continuum theory of Eringen, it was proposed that microstructure of materials could be represented in a continuum framework using a micro-deformation tensor governing micro-element deformation, in addition to the deformation gradient governing macro-element deformation. The paper formulates finite strain micromorphic elastoplasticity based on micromorphic continuum mechanics in the sense of Eringen. Multiplicative decomposition into elastic and plastic parts of the deformation gradient and micro-deformation are assumed, and the Clausius–Duhem inequality is formulated in the intermediate configuration B ¯ to analyze what stresses, elastic deformation measures, and plastic deformation rates are used/defined in the constitutive equations. The resulting forms of plastic and internal state variable evolution equations can be viewed as phenomenological at their various scales (i.e., micro-continuum and macro-continuum). The phenomenology of inelastic mechanical material response at the various scales can be different, but for demonstration purposes, J 2 flow plasticity is assumed for each of three levels of plastic evolution equations identified, with different stress, internal state variables, and material parameters. All evolution equations and a semi-implicit time integration scheme are formulated in the intermediate configuration for future coupled Lagrangian finite element implementation. A simpler two-dimensional model for anti-plane shear kinematics is formulated to demonstrate more clearly how such model equations simplify for future finite element implementation.

Equation-Free/Galerkin-Free Reduced-Order Modeling of the Shallow Water Equations Based on Proper Orthogonal Decomposition

In this paper, two categories of reduced-order modeling (ROM) of the shallow water equations (SWEs) based on the proper orthogonal decomposition (POD) are presented. First, the traditional Galerkin-projection POD/ROM is applied to the one-dimensional (1D) SWEs. The result indicates that although the Galerkin-projection POD/ROM is suitable for describing the physical properties of flows (during the POD basis functions’ construction time), it cannot predict that the dynamics of the shallow water flows properly as it was expected, especially with complex initial conditions. Then, the study is extended to applying the equation-free/Galerkin-free POD/ROM to both 1D and 2D SWEs. In the equation-free/Galerkin-free framework, the numerical simulation switches between a fine-scale model, which provides data for construction of the POD basis functions, and a coarse-scale model, which is designed for the coarse-grained computational study of complex, multiscale problems like SWEs. In the present work, the Beam & Warming and semi-implicit time integration schemes are applied to the 1D and 2D SWEs, respectively, as fine-scale models and the coefficients of a few POD basis functions (reduced-order model) are considered as a coarse-scale model. Projective integration is applied to the coarse-scale model in an equation-free framework with a time step grater than the one used for a fine-scale model. It is demonstrated that equation-free/Galerkin-free POD/ROM can resolve the dynamics of the complex shallow water flows. Moreover, the computational cost of the approach is less than the one for a fine-scale model.

On a solenoidal Fourier–Chebyshev spectral method for stability analysis of the Hagen–Poiseuille flow

In this work we study the efficiency of a spectral Petrov–Galerkin method for the linear and nonlinear stability analysis of the pipe or Hagen–Poiseuille flow. We formulate the problem in solenoidal primitive variables for the velocity field and the pressure term is eliminated from the scheme suitably projecting the equations on another solenoidal subspace. The method is unusual in being based on Chebyshev polynomials of selected parity for the radial variable, avoiding clustering of the quadrature points near the origin, satisfying appropriate regularity conditions at the pole and allowing the use of a fast cosine transform if required. Besides, this procedure provides good conditions for the time marching schemes. For the time evolution, we use semi-implicit time integration schemes. Special attention is given to the explicit treatment and efficient evaluation of the nonlinear terms via pseudospectral partial summations. The method provides spectral accuracy and the linear and nonlinear results obtained are in very good agreement with previous works. The scheme presented can be applied to other flows in unbounded cylindrical geometries.

The NCAR Spectral Element Climate Dynamical Core: Semi-Implicit Eulerian Formulation

A prototype dynamical core for the Community Atmospheric Model (CAM) component of the Community Climate System Model (CCSM) is presented. The 3D governing primitive equations are specified in curvilinear coordinates on the cubed-sphere combined with a hybrid pressure ? vertical coordinate. The horizontal space discretisation is based on a $$\mathbb{P}_N - \mathbb{P}_N$$ spectral element variational formulation. A semi-implicit time integration scheme is derived in order to circumvent the severe time step restrictions associated with gravity waves. Eigen-mode decomposition of the vertical structure matrix results in a set of decoupled 2D modified Helmholtz problems which are solved using a preconditioned conjugate gradient iteration. An idealized climate simulation is presented, where the semi-implicit scheme permits a much larger time step