Implicit Time-stepping Scheme Research Articles

Dynamic modeling of a sliding ring on an elastic rod with incremental potential formulation

Abstract Mechanical interactions between rigid rings and flexible cables find broad application in both daily life (hanging clothes) and engineering system (closing a tether-net). A reduced-order method for the dynamic analysis of sliding rings on a deformable one-dimensional (1D) rod-like object is proposed. In contrast to the conventional approach of discretising joint rings into multiple nodes and edges for contact detection and numerical simulation, a single point is used to reduce the order of the model. To ensure that the sliding ring and flexible rod do not deviate from their desired positions, a new barrier function is formulated using the incremental potential theory. Subsequently, the interaction between tangent frictional forces is obtained through a delayed dissipative approach. The proposed barrier functional and the associated frictional functional are C2 continuous, hence the nonlinear elastodynamic system can be solved variationally by an implicit time-stepping scheme. The numerical framework is initially applied to simple examples where the analytical solutions are available for validation. Then, multiple complex practical engineering examples are considered to showcase the effectiveness of the proposed method. The simplified ring-to-rod interaction model has the capacity to enhance the realism of visual effects in image animations, while simultaneously facilitating the optimisation of designs for space debris removal systems.

A general positivity-preserving algorithm for implicit high-order finite volume schemes solving the Euler and Navier-Stokes equations

This paper presents a general positivity-preserving algorithm for implicit high-order finite volume schemes that solve compressible Euler and Navier-Stokes equations to ensure the positivity of density and internal energy (or pressure). Previous positivity-preserving algorithms are mainly based on the slope limiting or flux limiting technique, which rely on the existence of low-order positivity-preserving schemes. This dependency poses serious restrictions on extending these algorithms to temporally implicit schemes since it is difficult to know if a low-order implicit scheme is positivity-preserving. In the present paper, a new positivity-preserving algorithm is proposed in terms of the flux correction technique. And the factors of the flux correction are determined by a residual correction procedure. For a finite volume scheme that is capable of achieving a converged solution, we show that the correction factors are in the order of unity with additional high-order terms corresponding to the spatial and temporal rates of convergence. Therefore, the proposed positivity-preserving algorithm is accuracy-reserving and asymptotically consistent. The notable advantage of this method is that it does not rely on the existence of low-order positivity-preserving baseline schemes. Therefore, it can be applied to the implicit schemes solving Euler and especially Navier-Stokes equations. In the present paper, the proposed technique is applied to an implicit dual time-stepping finite volume scheme with temporal second-order and spatial high-order accuracy. The present positivity-preserving algorithm is implemented in an iterative manner to ensure that the dual time-stepping iteration will converge to the positivity-preserving solution. Another similar correction technique is also proposed to ensure that the solution remains positivity-preserving at each sub-iteration. Numerical results demonstrate that the proposed algorithm preserves positive density and internal energy in all test cases and significantly improves the robustness of the numerical schemes.

Stabilized state-based peridynamics for elasticity emanating from constrained Lagrangian

State-based peridynamics for elasticity paves the way to bridge the local theory with the inherent capabilities of general elasticity with non-local effects. However, the original non-ordinary state-based peridynamics with correspondence materiel models suffers from high-frequency oscillations due to zero-energy mode instability. To address this issue, we start from the fundamental Lagrangian with the introduction of a constraint potential energy that minimizes the numerical error and suppresses the zero-energy mode oscillations. A stabilized non-ordinary state-based peridynamics is then derived from this constrained Lagrangian that forms an integral–differential–algebraic system. In specifics, the equation of motion for particles is described by integral-differential equations while the constraint suppressing the zero-energy mode instability is expressed as algebraic equations (zero-error constraint). An energy–momentum conserving implicit time-stepping scheme for this integral–differential–algebraic system is then derived. End-to-end 2D and 3D simulations, including scalar wave propagation, crack path, and large deformation problems, are presented to show that the proposed method is promising for these applications in terms of accurately capturing the crack path and large deformation, efficiently eliminating the numeral oscillations, conserving the total energy for long duration simulations, and preserving the zero-error constraint.

A compatible finite element discretisation for the nonhydrostatic vertical slice equations

We present a compatible finite element discretisation for the vertical slice compressible Euler equations, at next-to-lowest order (i.e., the pressure space is bilinear discontinuous functions). The equations are numerically integrated in time using a fully implicit timestepping scheme which is solved using monolithic GMRES preconditioned by a linesmoother. The linesmoother only involves local operations and is thus suitable for domain decomposition in parallel. It allows for arbitrarily large timesteps but with iteration counts scaling linearly with Courant number in the limit of large Courant number. This solver approach is implemented using Firedrake, and the additive Schwarz preconditioner framework of PETSc. We demonstrate the robustness of the scheme using a standard set of testcases that may be compared with other approaches.

Sparse grid time-discontinuous Galerkin method with streamline diffusion for transport equations

High-dimensional transport equations frequently occur in science and engineering. Computing their numerical solution, however, is challenging due to its high dimensionality. In this work we develop an algorithm to efficiently solve the transport equation in moderately complex geometrical domains using a Galerkin method stabilized by streamline diffusion. The ansatz spaces are a tensor product of a sparse grid in space and discontinuous piecewise polynomials in time. Here, the sparse grid is constructed upon nested multilevel finite element spaces to provide geometric flexibility. This results in an implicit time-stepping scheme which we prove to be stable and convergent. If the solution has additional mixed regularity, the convergence of a 2d-dimensional problem equals that of a d-dimensional one up to logarithmic factors. For the implementation, we rely on the representation of sparse grids as a sum of anisotropic full grid spaces. This enables us to store the functions and to carry out the computations on a sequence regular full grids exploiting the tensor product structure of the ansatz spaces. In this way existing finite element libraries and GPU acceleration can be used. The combination technique is used as a preconditioner for an iterative scheme to solve the transport equation on the sequence of time strips. Numerical tests show that the method works well for problems in up to six dimensions. Finally, the method is also used as a building block to solve nonlinear Vlasov-Poisson equations.

A hybrid lattice Boltzmann - Navier-Stokes method for unsteady aerodynamic and aeroacoustic computations

A hybrid numerical method coupling the standard lattice Boltzmann method and a compressible finite-volume Navier-Stokes solver is proposed in the context of unsteady aerodynamic and aeroacoustic simulations. The trend being towards more realistic and detailed simulations in a reasonable amount of CPU time, lattice Boltzmann and Navier-Sokes methods can be combined to solve the same problem. The present hybrid method relies on a zonal decomposition of the computational domain thus allowing to exploit the numerical features of both methods to their optimal extent in specific flow regions.The key issue when combining the lattice Boltzmann method with a Navier-Stokes solver is to ensure a smooth transition of the flow variables at the two-way coupling interface. While existing approaches consider overlapping meshes, a direct grid coupling is here derived. The mapping from the macroscopic flow variables to the set of lattice Boltzmann distribution functions is performed analytically thanks to a Chapman-Enskog expansion and draws a direct link to advanced regularised collision operators. Unsteady computations are enabled by coupling the lattice Boltzmann stream and collide algorithm with explicit and implicit Navier-Stokes time-stepping schemes. The hybrid method is then assessed on four time-dependent test cases representative of aerodynamic and aeroacoustic problems. The proposed approach is proven to yield very accurate results while keeping the numerical advantages of both methods and reducing the overall computational cost of direct noise computations.

: A high-order discontinuous Galerkin solver for flow simulations and multi-physics applications

We present the latest developments of our High-Order Spectral Element Solver (▪), an open source high-order discontinuous Galerkin framework, capable of solving a variety of flow applications, including compressible flows (with or without shocks), incompressible flows, various RANS and LES turbulence models, particle dynamics, multiphase flows, and aeroacoustics. We provide an overview of the high-order spatial discretisation (including energy/entropy stable schemes) and anisotropic p-adaptation capabilities. The solver is parallelised using MPI and OpenMP showing good scalability for up to 1000 processors. Temporal discretisations include explicit, implicit, multigrid, and dual time-stepping schemes with efficient preconditioners. Additionally, we facilitate meshing and simulating complex geometries through a mesh-free immersed boundary technique. We detail the available documentation and the test cases included in the GitHub repository. Program summaryProgram Title:▪CPC Library link to program files:https://doi.org/10.17632/738py2shk4.1Developer's repository link:https://github.com/loganoz/horses3d.Licensing provisions: MITProgramming language: Fortran 2008External routines/libraries: METIS, MPI, HDF5, MKL, PETSc (all are optional).Nature of problem:▪ is a high-order discontinuous Galerkin framework, capable of solving a variety of flow applications, including compressible flows (with or without shocks), incompressible flows, various RANS and LES turbulence models, particle dynamics, multiphase flows, and aeroacoustics.Solution method: high-order discontinuous Galerkin Spectral Element (DGSEM) and explicit/implicit time-steppers.Additional comments including restrictions and unusual features:▪ is a multiphysics environment where the compressible Navier-Stokes equations, the incompressible Navier–Stokes equations, the Cahn–Hilliard equation and entropy–stable variants are solved. Arbitrary high–order, p–anisotropic discretisations are used, including static and dynamic p–adaptation methods (feature-based and truncation error-based). Explicit and implicit time-steppers for steady and time-marching solutions are available, including efficient multigrid and preconditioners. Numerical and analytical Jacobian computations with a colouring algorithm have been implemented. Multiphase flows are solved using a diffuse interface model: Navier–Stokes/Cahn–Hilliard. Turbulent models implemented include RANS: Spalart-Allmaras and LES: Smagorinsky, Wale, Vreman; including wall models. Immersed boundary methods can be used, to avoid creating body fitted meshes. Acoustic propagation can be computed using Ffowcs-Williams and Hawkings models. ▪ supports curvilinear, hexahedral, conforming meshes in GMSH, HDF5 and SpecMesh/HOHQMesh format. A hybrid CPU-based parallelisation strategy (shared and distributed memory) with OpenMP and MPI is followed.

High-order implicit RBF-based differential quadrature-finite volume method on unstructured grids: Application to inviscid and viscous compressible flows

This paper exploits the potential of a high-order implicit radial basis function-based differential quadrature-finite volume (IRBFDQ-FV) method for effective simulation of inviscid and viscous compressible flows using unstructured grids. The framework of IRBFDQ-FV method is based on finite volume discretization which can guarantee conservation of mass, momentum and energy. The flow field variables within each control cell are represented by a fourth-order approximation constructed by a Taylor series expansion to the cell center with spatial derivatives as the undetermined coefficients. All spatial derivatives are computed by the meshless and highly converged RBF-based differential quadrature (RBFDQ) technique. Regarding flux evaluation, the discrete gas-kinetic flux solver is applied to evaluate inviscid and viscous fluxes simultaneously for compressible viscous flows. Resultantly, the special mathematical treatment for the viscous flux, which is usually adopted in other high-order methods, can be avoided. A point extrema-based extended bounds (PEEB) shock-capturing limiter is introduced to eliminate spurious numerical oscillations near discontinuities. Due to the implicit nature of the proposed method, implicit time-stepping schemes compatible with the spatial accuracy are devised to efficiently solve the resulting discretized equations. Representative compressible flow problems are simulated to verify the high-order accuracy, excellent efficiency and robustness of the present method. Comparison with two other high-order finite volume methods further shows competitiveness of the present method in solving compressible flow problems on unstructured grids.

Numerical Simulation of Two-Phase Flow in Liquid Composite Moulding Using VOF-Based Implicit Time-Stepping Scheme

The filling stage in injection/infusion moulding processes plays a key role in composite manufacturing that can be influenced by the inlet and vent ports. This will affect the production of void-free parts and the desirable process time. Flow control is usually required in experiments to optimise such a stage; however, numerical simulations can be alternatively used to predict manufacturing-induced deficiencies and potentially remove them in the actual experiments. This study uses ANSYS Fluent software to model flow-front advancement during the impregnation of woven fabrics. A developed technique is applied by creating tracking points (e.g., on-line monitor) in the direction of the flow to report/collect data for flow-front positions as a function of time. The study adopts the FVM-VOF-based two-phase flow model together with an implicit time-stepping scheme, i.e., a dual-time formulation solution method with a preconditioned pseudo-time derivative. Initially, three time-step sizes, 5 s (small), 25 s, and 50 s (large), are evaluated to examine their impact on numerical saturation lines at various fabric porosities, 40%, 50%, and 60%, for a two-dimensional (2D) rectangular mould, and predictions are then compared with the well-known analytical Darcy. This is followed by a three-dimensional (3D) curved mould for a fillet L-shaped structure, wherein the degree-of-curvature of fibre preforms is incorporated using a User-Defined Function (UDF) to tailor the impregnation process. The developed approach shows its validation (1–5.7%) with theoretical calculations and experimental data for 2D and 3D cases, respectively. The results also stress that a shorter computational time can be achieved with a large time-step size while maintaining the same level of accuracy.

A modified thermodynamic sea ice model and its application

A modified thermodynamic sea ice model suitable for large-scale climate simulations is described. Originated from the Winton’s three-layer model framework, this new model includes several improvements in the vertical thermodynamics: (1) the number of ice layers increases from two to three; (2) the snow heat capacity is included; (3) a vertically varying salinity profile is implemented; and (4) a temperature- and salinity-dependent heat conductivity parameterization scheme is introduced. A non-iterative fully implicit time-stepping scheme similar to Winton’s model is used to calculate the temperature of ice and snow. Results from a series of one-dimensional experiments show that equilibrium ice thickness in the modified model is increased by 45 cm when compared with the original Winton’s model. All modifications mentioned above contribute to this change of ice thickness, among which the increase of ice layer has the most significant effect. Experiments using the Modular Ocean Model version 4 (MOM4) coupled with the modified model show an improved sea ice simulation which includes an increase in both the sea ice volume and thickness over the entire Arctic region, confirming the above founding. However, contrary model behavior exhibits when the snow heat capacity is considered that warrants further investigation.

Price options on investment project expansion under commodity price and volatility uncertainties using a novel finite difference method

In this paper we develop a PDE-based mathematical model for valuing real options on the expansion of an investment project whose underlying commodity price and its volatility follow their respective geometric Brownian motions. This mathematical model is of the form of a 2-dimensional Black-Scholes equation whose payoff condition is determined also by a PDE system. A novel 9-point finite difference scheme is proposed for the discretization of the spatial derivatives and the fully implicit time-stepping scheme is used for the time discretization of the PDE systems. We show that the coefficient matrix of the fully discretized system is an M-matrix and prove that the solution generated by this finite difference scheme converges to the exact one when the mesh sizes approach zero. To demonstrate the usefulness and effectiveness of the mathematical model and numerical method, we present a case study on a real option pricing problem in the iron-ore mining industry. Numerical experiments show that our model and methods are able to produce numerical results which are financially meaningful.

Metode Numerik Untuk Menentukan Harga Opsi Dengan Model Volatilitas Leland

Option price under transaction cost with leland volatility model is the solution of a non linear diferential equations. To solve this equation used numerical methods based on an upwind finite difference for spatial discretization as well as the use of explicit and implicit methods for discretizing time-stepping. upwind finite difference method with explicit time-stepping scheme proved to be unstable so as not konvegen. While the use of implicit time-stepping scheme is proved monotonous, consistent and stable so that converge to the viscosity solution.

Convergence acceleration for high-order shock-fitting methods in hypersonic flow applications with efficient implicit time-stepping schemes

The research on hypersonic boundary layer transition is a recent focal point in fluid dynamics. One of the effective numerical methods for transition study is the high-order shock-fitting method, which is especially suitable for dealing with shock-related interactions and removing post-shock oscillations as compared with shock-capturing methods. Currently, explicit time-stepping schemes are mainly used for shock-fitting methods. In comparison, implicit schemes have received less attention due to the extra complexity arising from dynamic shock and moving mesh, even though they can accelerate the convergence speeds for steady flows by up to an order of magnitude. In this paper, two implicit time-stepping methods, GMRES (generalized minimum residual) and LR (line-relaxation), are realized for the efficient calculation of steady flow fields using the 5th-order shock-fitting method. The gas models include both calorically perfect gas and 5-species thermal-chemical non-equilibrium gas required for high-speed and high-enthalpy flows. The results of test cases showed that the convergence rates of GMRES+LR and LR are 4 to 20 times faster than the explicit/semi-implicit Runge-Kutta schemes. Therefore, the usage of implicit GMRES and LR as time-stepping schemes raises the computational efficiency for obtaining steady boundary layer flows, which helps the important high-order shock-fitting methods play a greater role in the study of hypersonic boundary layer stability and transition.

Structure Preserving Discretization of Allen–Cahn Type Problems Modeling the Motion of Phase Boundaries

We study the systematic numerical approximation of a class of Allen–Cahn type problems modeling the motion of phase interfaces. The common feature of these models is an underlying gradient flow structure which gives rise to a decay of an associated energy functional along solution trajectories. We first study the discretization in space by a conforming Galerkin approximation of a variational principle which characterizes smooth solutions of the problem. Well-posedness of the resulting semi-discretization is established and the energy decay along discrete solution trajectories is proven. A problem adapted implicit time-stepping scheme is then proposed and we establish its well-posed and decay of the free energy for the fully discrete scheme. Some details about the numerical realization by finite elements are discussed, in particular the iterative solution of the nonlinear problems arising in every time-step. The theoretical results are illustrated by numerical tests which also provide further evidence for asymptotic expansions of the interface velocities derived by Alber et al. and support the observation that their hybrid Allen–Cahn model avoids the problem of mesh-locking to a large extent.

Time-stepping DPG formulations for the heat equation

For a wide range of PDEs, the discontinuous Petrov–Galerkin (DPG) methodology of Demkowicz and Gopalakrishnan provides discrete stability starting from a coarse mesh and minimization of the residual in a user-controlled norm, among other appealing features. Research on DPG for transient problems has mainly focused on spacetime discretizations, which has theoretical advantages, but practical costs for computations and software implementations. The sole examination of time-stepping DPG formulations was performed by Führer, Heuer, and Gupta, who applied Rothe’s method to an ultraweak formulation of the heat equation to develop an implicit time-stepping scheme; their work emphasized theoretical results, including error estimates in time and space.In the present work, we follow Führer, Heuer, and Gupta in examining the heat equation; our focus is on numerical experiments, examining the stability and accuracy of several formulations, including primal as well as ultraweak, and explicit as well as implicit and Crank–Nicolson time-stepping schemes. We are additionally interested in communication-avoiding algorithms, and we therefore include a highly experimental formulation that places all the trace terms on the right-hand side of the equation.

Chemotaxis on networks: Analysis and numerical approximation

We consider the Keller–Segel model of chemotaxis on one-dimensional networks. Using a variational characterization of solutions, positivity preservation, conservation of mass, and energy estimates, we establish global existence of weak solutions and uniform bounds. This extends related results of Osaki and Yagi to the network context. We then analyze the discretization of the system by finite elements and an implicit time-stepping scheme. Mass lumping and upwinding are used to guarantee the positivity of the solutions on the discrete level. This allows us to deduce uniform bounds for the numerical approximations and to establish order optimal convergence of the discrete approximations to the continuous solution without artificial smoothness requirements. In addition, we prove convergence rates under reasonable assumptions. Some numerical tests are presented to illustrate the theoretical results.

Second-order cone interior-point method for quasistatic and moderate dynamic cohesive fracture

Cohesive fracture is among the few techniques able to model complex fracture nucleation and propagation with a sharp (nonsmeared) representation of the crack. Implicit time-stepping schemes are often favored in mechanics due to their ability to take larger time steps in quasistatic and moderate dynamic problems. Furthermore, initially rigid cohesive models are typically preferred when the location of the crack is not known in advance, since initially elastic models artificially lower the material stiffness. It is challenging to include an initially rigid cohesive model in an implicit scheme because the initiation of fracture corresponds to a nondifferentiability of the underlying potential. In this work, an interior-point method is proposed for implicit time stepping of initially rigid cohesive fracture. It uses techniques developed for convex second-order cone programming for the nonconvex problem at hand. The underlying cohesive model is taken from Papoulia (2017) and is based on a nondifferentiable energy function. That previous work proposed an algorithm based on successive smooth approximations to the nondifferential objective for solving the resulting optimization problem. It is argued herein that cone programming can capture the nondifferentiability without smoothing, and the resulting cone formulation is amenable to interior-point algorithms. A further benefit of the formulation is that other conic inequality constraints are straightforward to incorporate. Computational results are provided showing that certain contact constraints can be easily handled and that the method is practical.

Simulation of shell buckling by implicit dynamics and numerically dissipative schemes

Shell buckling process is studied by implicit structural dynamics time-stepping schemes with numerical dissipation (in the high-frequency range), which fall into the class of Generalized energy-momentum method. In particular, modified Generalized-α, Bossak, HHT, and Newmark’s schemes are considered, as well as Energy-momentum conserving scheme, and Energy-decaying scheme with controllable numerical dissipation. We are interested in the assessment of the ability of these schemes to handle complex buckling and post-buckling processes of thin shells (even for cases when path-following methods fail). The schemes are specialized for geometrically exact shell formulation (with only displacement-like degrees of freedom). Computed numerical examples include classical shell buckling problems: snap-through, shell collapse, and buckling of perfect and imperfect cylinder under axial load. The examples illustrate that (high-frequency) numerical dissipation is absolutely necessary for an efficient implicit dynamic simulation of complex shell buckling and post-buckling processes. Since the Energy-momentum conserving scheme does not damp high-frequency oscillations that accompany shell buckling (and thus accumulates large error in the high-frequency range), it is not suitable for the simulation of shell instabilities (the scheme fails to converge for several computed examples). The dynamic results are compared with the static ones that were computed by the path-following method. It turns out that implicit dynamic analysis with (high-frequency) numerical dissipation is considerably more robust and efficient than static analysis for several computed examples, including snap trough problems and buckling of cylinder under axial load.

An implicit time-stepping scheme and an improved contact model for ice-structure interaction simulations

This paper presents a novel time-stepping scheme for the modelling of discrete ice-structure interaction. The scheme extends the non-smooth discrete element modelling (NDEM) technique to enable compliant continuous and discontinuous contacts. This increases the accuracy and expands the applicability range of the NDEM technique related to ice-structure interaction problems. We derive the parameters representing the compliant behaviour of contacts. The accuracy of the presented scheme for discontinuous contacts is compared to an existing, simpler scheme that limits the contact force based on a maximum force assumption. The comparison shows that the derived scheme results in more accurate contact forces, for the same time step size, as previously applied NDEM schemes in ice-structure interaction. An example simulation is compared against ice tank tests of a 4-legged, vertical-walled structure moving through a broken ice field.

Laser-induced heating of dynamic particulate depositions in additive manufacturing

Many applications in additive manufacturing involve the dynamic deposition of powders and the in-flight heating of such material by a laser. In order to characterize such systems, the Discrete Element Method (DEM) is employed. Specifically, this paper focuses on the initial stages of this process by developing a modular discrete-element type multiphysics simulation method for the particle dynamics and heating by a laser and detailed thermal behavior. The objective is to provide researchers with a framework to construct computational tools for this growing industry. In order to achieve this, from a simulation standpoint, the overall particle-mixture system is constructed by coupling submodels for each primary physical process (dynamics and heating) together. An iterative staggering scheme is developed whereby, within every time step, each individual particle is solved “freezing” the state of the remaining multiparticle system. The state of the particle is then updated and the algorithm moves to the next particle in the system and the process is repeated. The overall process sweeps through the entire system repeatedly until convergence is achieved in an appropriate norm. As the system evolves, an error estimate dictates the time-step size that is needed to induce convergence to below an appropriate error level. Thus, the process can be considered as an implicit time-stepping scheme, which is combined with an (internal) iterative staggering process. In order to control rates of convergence within a time-step, the algorithm adjusts the time-step size. If the iterative process does not converge within a desired number of iterations, below an error tolerance, the time-step is reduced. The degree of time-step reduction is determined by utilizing an estimate of the spectral radius of the coupled system. Since the construction of model and solution process is modular, one can easily replace physical submodels with other choices, making it easy to numerically experiment with a variety of models. Qualitative and quantitative analyses are provided, as well as three-dimensional numerical examples. This paper addresses the first of two overall “macrostages” of a complex manufacturing process with DEM. The first stage, addressed in this paper, concentrates on dry-powder materials, which are of a discrete particulate character, using DEM to characterize the deposition and in-flight particle heating. Other papers of the author (Zohdi, 2010, 2013, 2014, 2015 [81,82], 2008) have addressed stage two, namely phase transformations, curing (cooling) and stress analysis, using more appropriate continuum and hybrid DEM-continuum approaches.