Time Step Size Research Articles

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.

A GPU based accelerated solver for simulation of heat transfer during metal casting process

The metal casting process, which is one of the key drivers of the manufacturing industry, involves several physical phenomena occurring simultaneously like fluid flow, phase change, and heat transfer which affect the casting yield and quality. Casting process modeling involves numerical modeling of these phenomena on a computer. In recent decades, this has become an inevitable tool for foundry engineers to make defect-free castings. To expedite computational time graphics processing units (GPUs) are being increasingly used in the numerical modeling of heat transfer and fluid flow. Initially, in this work a CPU based implicit solver code is developed for solving the 3D unsteady energy equation including phase change numerically using finite volume method which predicts the thermal profile during solidification in the metal casting process in a completely filled mold. To address the computational bottleneck, which is identified as the linear algebraic solver based on the bi-conjugate gradient stabilized method, a GPU-based code is developed using Compute Unified Device Architecture toolkit and was implemented on the GPU. The CPU and GPU based codes are then validated against a commercial casting simulation code FLOW-3D CAST® for a simple casting part and against in-house experimental results for gravity die casting of a simple geometry. Parallel performance is analyzed for grid sizes ranging from 10 × 10 × 10 to 210 × 210 × 210 and for three time-step sizes. The performance of the GPU code based on occupancy and throughput is also investigated. The GPU code exhibits a maximum speedup of 308× compared to the CPU code for a grid size of 210 × 210 × 210 and a time-step size of 2 s.

Nonsmooth Reduced Interface Models and Their Use in Co-Simulation of Mechanical Systems

Abstract In a co-simulation setup, the entire system is decomposed into a collection of individual subsystems that are interfaced together, with each subsystem being modeled and integrated separately according to its own requirements. To maintain the interconnectivity and consolidation of the primary system, these subsystems must communicate with each other through the interface and transfer certain information at the end points of a defined time interval termed macro time step. Inside the macro time step, the evolution of the interface variables has to be approximated as information about them will only be available again at the end of the step. In real-time simulations, the size of the macro time step and the accuracy of the approximated interface variables are critical factors; if the interface variables are approximated accurately, the size of the macro time step can be kept large enough to provide interactive rates without loss of accuracy and stability. This work focuses on systems where unilateral contact interactions are important and proposes reduced interface model concepts for such nonsmooth systems. The use of the proposed reduced interface model (RIM) is demonstrated in co-simulation to provide model-based approximation of the interface variables. The advantages of the proposed method are demonstrated through two representative case studies.

The Finite Element Method in Time for Multibody Dynamics

Abstract The generalized-α scheme has become the approach of choice for the time integration of the equations of motion of multibody systems. Despite its simplicity, this scheme presents drawbacks: the time-step size cannot be changed easily, making it difficult to implement time adaptivity, and the solution of periodic problems cannot be found easily. This paper explores an alternative approach based on the finite element method in time. The basic principles underpinning the approach are presented and both time-continuous and time-discontinuous approaches are investigated. Two types of Galerkin schemes will be presented here: the time-continuous and the time-discontinuous schemes. In the former, the displacement field is continuous across interelement boundaries, whereas discontinuities or “jumps” are allowed across interelement boundaries for the latter. Simple problems are treated to identify the best schemes. Families of schemes of various accuracies are presented. The first family, based on time-continuous elements, features schemes that do not present numerical dissipation. Asymptotic annihilation is achieved by the time-discontinuous elements that form the second family. The problem of kinematic constraints is treated within the framework of the finite element method in time. Special emphasis is devoted to the satisfaction of the kinematic constraints and their time derivative within a time element.

Dynamics of the Fractional-Order Lorenz System Based on Adomian Decomposition Method and Its DSP Implementation

Dynamics and digital circuit implementation of the fractional-order Lorenz system are investigated by employing Adomian decomposition method (ADM). Dynamics of the fractional-order Lorenz system with derivative order and parameter varying is analyzed by means of Lyapunov exponents (LEs), bifurcation diagram, chaos diagram and phase diagram. Results show that the fractional-order Lorenz system has rich dynamical behavior and it is a potential model for application. It is also found that the minimum order is affected by numerical algorithm and time step size. Finally, the fractional-order system is implemented on DSP digital circuit. Phase diagrams generated by the DSP are consistent with that generated by simulation.

A high-order multi-time-step scheme for bond-based peridynamics

A high-order multi-time-step scheme (MTS) for the bond-based peridynamic (PD) model, an extension of classical continuous mechanics widely used for analyzing discontinuous problems like cracks, is proposed. The MTS scheme discretizes the spatial domain with a meshfree method and advances in time with a high-order Runge–Kutta method. To effectively handle discontinuities (cracks) that appear in a local subdomain in the solution, the scheme employs the Taylor expansion and Lagrange interpolation polynomials with a finer time step size, that is, coarse and fine time step sizes for smooth and discontinuous subdomains, respectively, to achieve accurate and efficient simulations. By eliminating unnecessary fine-scale resolution imposed on the entire domain, the MTS scheme outperforms the standard STS scheme for PD by significantly reducing computational costs, particularly for problems with discontinuous solutions, as demonstrated by comprehensive theoretical analysis and numerical experiments.

Determination of Maneuvering Force Coefficients for a Destroyer Model with OpenFOAM

_ Ship maneuvering performance can be predicted using various methods, including physical model tests, rapid simulations using coefficient-based forces, and high-fidelity computational fluid dynamics. This article considers the application of computational fluid dynamics for the prediction of hull force coefficients to be used as inputs to rapid maneuvering simulations. The open-source software OpenFOAM was used to simulate forces on the destroyer model DTMB 5415 for steady drift, oscillatory pure sway motion, and oscillatory yaw motion. Recommendations are provided regarding best practices in a number of areas, including domain size, mesh refinement, time step size, and turbulent modeling. Predicted forces for steady drift angles up to 12 degrees are typically within 10%of experimental values. For oscillatory sway and yaw motions, predicted forces are typically within 25%of experimental values. Keywords maneuvering; OpenFOAM; best modeling practices

Comprehensive investigation of the explicit, positivity preserving methods for the heat equation

In this paper-series, we investigate the performance of 12 explicit non-conventional algorithms. All of them have the convex combination property, thus they are unconditionally stable and preserve the positivity of the solution when they applied to the heat equation. In this part of the series, we construct several 2D systems to find how the errors depend on the time step size. Sweeps for other key parameters will be presented in the next part of the series.

Comprehensive investigation of the explicit, positivity preserving methods for the heat equation

In this paper-series, we investigate the performance of 12 explicit non-conventional algorithms in several 2D systems. All of them have the convex combination property, thus they are unconditionally stable and preserve the positivity of the solution when they are applied to the heat equation. In the first part of the series, we examined how the errors depend on the time step size and running times. Now we present additional numerical test results, where sweeps for parameters such as the stiffness and the wavelength of the initial function will be performed.

A highly parallel algorithm for simulating the elastodynamics of a patient-specific human heart with four chambers using a heterogeneous hyperelastic model

In this paper a highly parallel method is developed for simulating the elastodynamics of a four-chamber human heart with patient-specific geometry. The heterogeneous hyperelastic model is discretized by a finite element method in space and a fully implicit adaptive method in time, and the resulting nonlinear algebraic systems are solved by a scalable domain decomposition algorithm. The deformations of the cardiac muscles are quite complex due to the realistic geometry, the heterogeneous hyperelasticity of the cardiac tissue, and the myocardial fibers with active stresses. Moreover, the deformations in different chambers and at different phases of the cardiac cycle are very different. To simulate all the muscle movements including the atrial diastole, the atrial systole, the isovolumic contraction, the ventricular ejection, the isovolumic relaxation, and the ventricular filling, the temporal-spatial mesh needs to be sufficiently fine, but not too fine so that the overall computing time is manageable, we introduce a baseline mesh in space and a two-level time stepping strategy including a uniform baseline time step size to obtain the desired time accuracy and an adaptive time stepping method within a baseline time step to guarantee the convergence of the nonlinear solver. Through numerical experiments, we investigate the performance of the proposed method with respect to the material coefficients, the fiber orientations, as well as the mesh sizes and the time step sizes. For an unstructured tetrahedral mesh with more than 200 million degree of freedoms, the method scales well for up to 16,384 processor cores for all steps of an entire cardiac cycle.

One-Way Fluid–Structure Interaction Simulation of the Added Damping of T(0,1) Modes of Two Partially Overlapped Identical Plates in Still Water

Abstract Structures with a partially overlapped status in water can be seen in some engineering applications, and the fluid-structure coupling vibration behavior of two partially overlapped identical plates has been studied previously through the experimental method. In this study, the added damping of the T(0,1) modes of two submerged partially overlapped plates is numerically investigated through the one-way fluid-structure interaction (FSI) method. The relevant numerical settings, like the mesh size, calculated periods, time-step size, and vibration amplitude, were tested first. Then, the numerical results were compared with experimental results, and good agreements were found. Finally, numerical results were analyzed. The vibration status of two plates in the joint abutting area or the overlapped area has an important influence on the added mass variation. When the added mass is higher, the phase difference between modal force and vibration displacement is also greater, which is the main reason for the higher added damping. The relationship between the phase difference and the frequency in water can be approximately fitted to a straight line, which can probably be used to predict the added damping variations caused by fluid boundary changes of submerged structures.

Unconditional stability and error estimates of the modified characteristics FEMs for the Micropolar Navier–Stokes Equations

In this paper, the unconditional stability and optimal error estimate of the velocity, pressure and angular velocity for the modified characteristics FEMs of the unsteady Micropolar Naiver–Stokes Equations (MNSE) are presented. In this method, the nonlinear equation is linearized by the characteristic finite element method for dealing with the time derivative term and the convection term. Basing on the characteristic time-discrete system, the error function is split into a temporal error and a spatial error. With a rigorous analysis to the characteristic time-discrete system, we prove that the error between the numerical solution and the solution of the time-discrete system is τ-independent, where τ denotes the time stepsize. The stability results and optimal error estimates in L2 norm and H1 norm will be given. Finally, some numerical results will be provided to confirm our theoretical analysis.

Improved dynamic analysis for plane elasticity problems with a boundary meshfree method and distributed sources

In this paper, we present an improved dynamic analysis method to address plane elasticity problems. Distributed sources are employed in a boundary-free mesh-based approach known as Boundary Distributed Source (BDS). This method, which builds upon a modified Method of Fundamental Solutions (MFS), aims to eliminate singularities without the reliance on fictitious boundaries. Distributed sources are strategically positioned directly on the real boundary, centered over disks. The method's efficiency and simplicity are highlighted by its ability to operate without the need for precise data concerning adjacent points. The dual reciprocity method (DRM) is utilized to approximate the non-homogeneous inertial term. Furthermore, the present formulation investigates several factors that influence the convergence of the proposed method, including the number of collocation points, internal points, the radius of the circular disk, and the time step size. Based on the numerical results, it is evident that the method exhibits remarkable accuracy and efficiency when applied to plane elasticity problems including forced, harmonic, and free analyses. As compared to conventional numerical methods for dynamic analysis, the method outperforms in both accuracy and convergence rate.

Monolithic and local time-stepping decoupled algorithms for transport problems in fractured porous media

Abstract The objective of this paper is to develop efficient numerical algorithms for the linear advection-diffusion equation in fractured porous media. A reduced fracture model is considered where the fractures are treated as interfaces between subdomains and the interactions between the fractures and the surrounding porous medium are taken into account. The model is discretized by a backward Euler upwind-mixed hybrid finite element method in which the flux variable represents both the advective and diffusive fluxes. The existence, uniqueness, as well as optimal error estimates in both space and time for the fully discrete coupled problem are established. Moreover, to facilitate different time steps in the fracture-interface and the subdomains, global-in-time, nonoverlapping domain decomposition is utilized to derive two implicit iterative solvers for the discrete problem. The first method is based on the time-dependent Steklov–Poincaré operator, while the second one employs the optimized Schwarz waveform relaxation (OSWR) approach with Ventcel-Robin transmission conditions. A discrete space-time interface system is formulated for each method and is solved iteratively with possibly variable time step sizes. The convergence of the OSWR-based method with conforming time grids is also proved. Finally, numerical results in two dimensions are presented to verify the optimal order of convergence of the monolithic solver and to illustrate the performance of the two decoupled schemes with local time-stepping on problems of high Péclet numbers.

Development and benchmarking of space-time kinetics solver in REDAC

A space–time kinetics solver based on finite volume approach has been developed for REactor Dynamics Analysis Code (REDAC). It solves 3D multi-group neutron diffusion equation, coupled with the heat transport equations. Code offers flexibility in terms of reactor geometry, boundary conditions and discretization. Efficient numerical methods have been employed for solving the governing equations. The performance of the code has been demonstrated by solving diverse benchmark problems. It is shown that REDAC is able to predict important static and dynamic core parameters with good accuracy. The framework of REDAC easily lends itself for sensitivity studies with respect to spatial discretization and temporal integration schemes. From these studies, it is observed that spatial resolution could have significant influence on accuracy and the semi-implicit method performs better than the implicit method. Improved quasi-static method requires less computational time compared to the direct method, however, its accuracy depends on the time step size.

Enhanced physics‐informed neural networks for efficient modelling of hydrodynamics in river networks

AbstractThis article enhances the physics‐informed neural networks (PINNs) method to effectively model the hydrodynamics of real‐world river networks with irregular cross‐sections. First, we pre‐process hydraulic parameters to optimize training speed without compromising accuracy, achieving a 91.67% acceleration compared with traditional methods. To address the vanishing gradient problem, layer normalization is also incorporated into the architecture. We also introduce novel physical constraints—water level range and junction node equations—to ensure effective training convergence and enrich the model with additional physical insights. Two practical case studies using HEC‐RAS benchmarks demonstrate that our improved PINN method can predict river network hydrodynamics with less data and is less sensitive to time step size, allowing for longer computational time steps. Incorporating physical knowledge, our enhanced PINN methodology emerges as an efficient and promising avenue for modelling the complexities of hydrodynamic processes in natural river networks.

Efficient discrete element modeling of heat generation and transfer in granular flows: Validation and application

This work presents an efficient Discrete Element Method (DEM) framework for the simulation of the thermal behavior of granular media. The main focus is on long-lasting granular flows, involving heat transfer and generation due to mechanical energy dissipation. The proposed approach uses efficient strategies to reduce the computational cost of the analyses and, therefore, to enable its application to problems of practical relevance. For instance, the contact area is adjusted to compensate for the artificial material softening that is typically considered in DEM to increase the time step size. After extended validation, the methodology is applied to the simulation of different setups of an experimental rotating drum. The numerical simulations presented good agreement with the experimental results and allowed a detailed analysis of the mechanisms and patterns of heat generation, which could not be extrapolated from the experimental campaign.

Physics-Informed Neural Networks for Solving High-Index Differential-Algebraic Equation Systems Based on Radau Methods

As is well known, differential algebraic equations (DAEs), which are able to describe dynamic changes and underlying constraints, have been widely applied in engineering fields such as fluid dynamics, multi-body dynamics, mechanical systems, and control theory. In practical physical modeling within these domains, the systems often generate high-index DAEs. Classical implicit numerical methods typically result in varying order reduction of numerical accuracy when solving high-index systems. Recently, the physics-informed neural networks (PINNs) have gained attention for solving DAE systems. However, it faces challenges like the inability to directly solve high-index systems, lower predictive accuracy, and weaker generalization capabilities. In this paper, we propose a PINN computational framework, combined Radau IIA numerical method with an improved fully connected neural network structure, to directly solve high-index DAEs. Furthermore, we employ a domain decomposition strategy to enhance solution accuracy. We conduct numerical experiments with two classical high-index systems as illustrative examples, investigating how different orders and time-step sizes of the Radau IIA method affect the accuracy of neural network solutions. For different time-step sizes, the experimental results indicate that utilizing a 5th-order Radau IIA method in the PINN achieves a high level of system accuracy and stability. Specifically, the absolute errors for all differential variables remain as low as 10−6, and the absolute errors for algebraic variables are maintained at 10−5. Therefore, our method exhibits excellent computational accuracy and strong generalization capabilities, providing a feasible approach for the high-precision solution of larger-scale DAEs with higher indices or challenging high-dimensional partial differential algebraic equation systems.

A new linearized second-order energy-stable finite element scheme for the nonlinear Benjamin-Bona-Mahony-Burgers equation

In this paper, a novel linearized second-order energy-stable fully discrete scheme for the nonlinear Benjamin-Bona-Mahony-Burgers (BBMB) equation is introduced, along with a superconvergence analysis of the conforming finite element method (FEM). Through skillful decomposition of the nonlinear term u∇⋅u, a new linearized second-order fully discrete scheme is developed. This scheme requires fewer iterations and exhibits higher computational efficiency compared to the nonlinear approach. Furthermore, it maintains energy stability, enabling direct numerical solution boundedness in the H1-norm, which represents an improvement over the L∞-norm boundedness reported in previous literature. Secondly, relying on this boundedness and the high-precision results of the bilinear element, we derive the unconditional error estimates for superclose and global superconvergence without any restrictions on the time step size Δt and spatial mesh size h. Finally, numerical examples are presented to validate the theoretical analysis.

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.