Critical Time Step Size Research Articles

An automatic simulation pipeline for coupled simulations of acoustic damping materials

AbstractFoamed materials are widely used to reduce noise due to their comparably good acoustic damping behavior. However, out of a large variety of these materials a suitable candidate has to be identified for each application. This is a challenging process that is typically guided by experiments and experience. While numerical simulations could support these experiments and reduce the effort to a great extent, no suitable discretization approach has yet been established that can fully capture the complex geometry of the foam. A fully resolved model is desirable in order to yield reliable predictions that can then be used to establish homogenized models. We established a monolithic coupling approach based on the finite cell method (FCM) that realizes a vibroacoustic simulation in this sense. The fluid and the structure domain are discretized by Cartesian grids and the geometry defined based on computed tomography scans is accounted for during the quadrature of the weak form. Our simulation in the time domain makes use of explicit time marching schemes and is therefore limited by a critical time step size. This is known to be arbitrarily low for discretizations with the FCM containing cells with arbitrarily small support. As a remedy against this we use the classical ‐stabilization technique and investigate its potentials and limitations.

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.

An eigenvalue stabilization technique for immersed boundary finite element methods in explicit dynamics

The application of immersed boundary methods in static analyses is often impeded by poorly cut elements (small cut elements problem), leading to ill-conditioned linear systems of equations and stability problems. While these concerns may not be paramount in explicit dynamics, a substantial reduction in the critical time step size based on the smallest volume fraction χ of a cut element is observed. This reduction can be so drastic that it renders explicit time integration schemes impractical. To tackle this challenge, we propose the use of a dedicated eigenvalue stabilization (EVS) technique.The EVS-technique serves a dual purpose. Beyond merely improving the condition number of system matrices, it plays a pivotal role in extending the critical time increment, effectively broadening the stability region in explicit dynamics. As a result, our approach enables robust and efficient analyses of high-frequency transient problems using immersed boundary methods. A key advantage of the stabilization method lies in the fact that only element-level operations are required.This is accomplished by computing all eigenvalues of the element matrices and subsequently introducing a stabilization term that mitigates the adverse effects of cutting. Notably, the stabilization of the mass matrix Mc of cut elements – especially for high polynomial orders p of the shape functions – leads to a significant raise in the critical time step size Δtcr.To demonstrate the efficiency of our technique, we present two specifically selected dynamic benchmark examples related to wave propagation analysis, where an explicit time integration scheme must be employed to leverage the increase in the critical time step size.

Implicit-explicit time integration for the immersed wave equation

Efficient solution strategies for wave propagation problems with complex geometries are essential in many engineering fields. Immersed boundary methods simplify mesh generation by embedding the domain of interest into an extended domain that is easy to mesh, introducing the challenge of dealing with cells that intersect the domain boundary. We consider the finite cell method that extends the weak form from the physical into the extended domain, multiplied by a small value, to stabilize badly cut cells. Combined with explicit time integration methods, the presence of finite cells with very little support in the physical domain results in tiny critical time step sizes. While the finite cell stabilization limits how small the critical time step size can become, this limit is still restrictive for many applications. Explicit transient analyses commonly use the spectral element method due to its natural way of obtaining diagonal mass matrices through nodal lumping. The resulting method is very efficient since nodal lumping renders the solution of the equation systems trivial while not reducing the accuracy or the critical time step size. The combination of the spectral element and finite cell methods is called the spectral cell method. Unfortunately, a direct application of nodal lumping in the spectral cell method is impossible due to the special quadrature necessary to treat the discontinuous integrand inside the cut cells. The existing approaches to lump the mass matrices of cut cells significantly reduce the accuracy of the approximation.We analyze an implicit-explicit (IMEX) time integration method to exploit the advantages of the nodal lumping scheme for uncut cells on one side and the unconditional stability of implicit time integration schemes for cut cells on the other. In this hybrid, immersed Newmark IMEX approach, we use explicit second-order central differences to integrate the uncut degrees of freedom that lead to a diagonal block in the mass matrix and an implicit trapezoidal Newmark method to integrate the remaining degrees of freedom (those supported by at least one cut cell). The immersed Newmark IMEX approach preserves the high-order convergence rates and the geometric flexibility of the finite cell method while retaining the efficiency of the nodal lumping scheme for uncut cells without compromising the critical time step size. We analyze a simple system of spring-coupled masses to highlight some of the essential characteristics of Newmark IMEX time integration. We then solve the scalar wave equation with immersed boundaries on two- and three-dimensional examples with significant geometric complexity to show that our approach is more efficient than state-of-the-art time integration schemes when comparing accuracy and runtime. While we focus on the finite and spectral cell methods, we expect other immersed methods to benefit equally from Newmark IMEX time integration.

A Hu-Washizu variational approach to self-stabilized quadrilateral Virtual Elements: 2D linear elastodynamics

A recent mixed formulation of the Virtual Element Method in 2D elastostatics, based on the Hu-Washizu variational principle, is here extended to 2D elastodynamics. The independent modeling of the strain field, allowed by the mixed formulation, is exploited to derive first order quadrilateral Virtual Elements (VEs) not requiring a stabilization (namely, self-stabilized VEs), in contrast to the standard VEs, where an artificial stabilization is always required for first order quads. Lumped mass matrices are derived using a novel approach, based on an integration scheme that makes use of nodal values only, preserving the correct mass in the case of rigid-body modes. In the case of implicit time integration, it is shown how the combination of a self-stabilized stiffness matrix with a self-stabilized lumped mass matrix can produce excellent performances both in the compressible and quasi-incompressible regimes with almost negligible sensitivity to element distortion. Finally, in the case of explicit dynamics, the performances of the different types of derived VEs are analyzed in terms of their critical time-step size.

Critical time-step size analysis and mass scaling by ghost-penalty for immersogeometric explicit dynamics

In this article, we study the effect of small-cut elements on the critical time-step size in an immersogeometric explicit dynamics context. We analyze different formulations for second-order (membrane) and fourth-order (shell-type) equations, and derive scaling relations between the critical time-step size and the cut-element size for various types of cuts. In particular, we focus on different approaches for the weak imposition of Dirichlet conditions: by penalty enforcement and with Nitsche’s method. The conventional stability requirement for Nitsche’s method necessitates either a cut-size dependent penalty parameter, or an additional ghost-penalty stabilization term. Our findings show that both techniques suffer from cut-size dependent critical time-step sizes, but the addition of a ghost-penalty term to the mass matrix serves to mitigate this issue. We confirm that this form of ‘mass-scaling’ does not adversely affect error and convergence characteristics for a transient membrane example, and has the potential to increase the critical time-step size by orders of magnitude. Finally, for a prototypical simulation of a Kirchhoff–Love shell, our stabilized Nitsche formulation reduces the solution error by well over an order of magnitude compared to a penalty formulation at equal time-step size.

A transient global-local generalized FEM for parabolic and hyperbolic PDEs with multi-space/time scales

This paper presents a novel Generalized Finite Element Method with global-local enrichment (GFEMgl) to solve time-dependent parabolic and hyperbolic problems with subscale features in both space and time. Compared with currently available GFEMgl, the proposed method requires a solution of local problems to solve transient PDEs, possibly with a different time integrator from the global problem or other local problems. The fine-scale solution can be easily retrieved. The proposed method is tested for the solution of the heat, advection, and advection-diffusion equations whose accuracy, stability, scalability, and convergence are thoroughly analyzed. Compared to direct analysis, numerical results show that the proposed method has the following properties 1) the accuracy closely matches direct analysis result with a fine mesh; 2) critical time step size is loosened; 3) optimal convergence rate is achieved. Moreover, the proposed method allows local problems to be turned off if their solutions are not helpful for the current global time step, which can further improve its efficiency.

Improved XFEM (IXFEM): 3D dynamic crack propagation under impact loading

An improved XFEM (IXFEM) for three-dimensional non-planar crack propagation under impact loadings is developed. The proposed method couples the IXFEM (Tian and Wen, 2015; Wen and Tian, 2016) and the B-spline ruled surface method (BRSM) proposed previously by authors (Xiao, Wen and Tian (2021)). The issues of the standard XFEM for dynamic crack propagation, such as energy inconsistency, “null” critical time step size and optimal mass lumping at crack tip, are fundamentally avoided with an extra-DOF free PU approximation (Tian, 2013). The tip singularity enrichment is introduced to assure high accuracy through the extra-DOF free PU approximation. The BRSM, which is designed to solve the challenge about robust and efficient geometry representation of 3D non-planar surface crack growth, is applied to model the geometry of the 3D dynamic crack propagation in this paper. The interaction integral for 2D dynamic stress intensity factor (SIF) evaluation is extended to 3D. Efficiency of the proposed approach is demonstrated via numerical tests.

An asynchronous parallel explicit solver based on scaled boundary finite element method using octree meshes

Explicit time integration methods are an integral part of solving structural dynamics problems such as the propagation of elastic and acoustic waves, impact scenarios including crash tests and many more. One common limitation of this class of time integrators is, however, their conditional stability, meaning that highly distorted, small, or stiff finite elements typically govern the critical time step size. That is to say, even a small number of such elements will result in a significantly decrease of the feasible time step, which in turn drastically increases the computational costs of solving the semi-discrete equations of motion. This is especially true for non-uniform and unstructured meshes. Therefore, an asynchronous explicit solver in parallel is proposed in this article. The idea is to assign different time step sizes to different parts of the mesh, such that the overall computational effort is minimized. To improve the performance of the proposed solver, it is combined with a sophisticated octree meshing framework, where the balanced octrees are used, meaning that the ratio of the element sizes of adjacent elements cannot exceed a value of two. Thus, the critical time step of each element can be estimated straightforwardly. To avoid issues with hanging nodes, a special polyhedral element formulation, the scaled boundary finite element method, is employed. Since there are only a limited number of cell patterns in a balanced octree, the stiffness and mass matrices are pre-computed, significantly reducing the computational cost. Moreover, exploiting an element-by-element technique, the assembly of global stiffness matrices can be avoided. The main advantage of the described methodology is that the asynchronous explicit solver can be easily implemented in a high-performance computing environment. By means of several numerical examples, the accuracy and efficiency of the proposed method are demonstrated, as well as its versatility in handling complex engineering problems.

Variationally consistent mass scaling for explicit time-integration schemes of lower- and higher-order finite element methods

In this paper, we propose a variationally consistent technique for decreasing the maximum eigenfrequencies of structural dynamics related finite element formulations. Our approach is based on adding a symmetric positive-definite term to the consistent mass matrix that follows from the integral of the traction jump across element boundaries. The added term is weighted by a small factor, for which we derive a suitable, and simple, element-local parameter choice. We perform numerical experiments for the linear wave equation in one and two dimensions, for quadrilateral elements and triangular elements, and for up to fourth order polynomial basis functions. Despite the increase in critical time-step size, we do not observe adverse effects in terms of spatial accuracy and orders of convergence. To extend the method to non-linear problems, we introduce a linear approximation. Our three-dimensional experiments with tetrahedral and hexahedral elements show that a sizeable increase in critical time-step size can be achieved while only causing minor (even beneficial) influences on the dynamic response.

A variational approach based on perturbed eigenvalue analysis for improving spectral properties of isogeometric multipatch discretizations

A key advantage of isogeometric discretizations is their accurate and well-behaved eigenfrequencies and eigenmodes. For degree two and higher, however, optical branches of spurious outlier frequencies and modes may appear due to boundaries or reduced continuity at patch interfaces. In this paper, we introduce a variational approach based on perturbed eigenvalue analysis that eliminates outlier frequencies without negatively affecting the accuracy in the remainder of the spectrum and modes. We then propose a pragmatic iterative procedure that estimates the perturbation parameters in such a way that the outlier frequencies are effectively reduced. We demonstrate that our approach allows for a much larger critical time-step size in explicit dynamics calculations. In addition, we show that the critical time-step size obtained with the proposed approach does not depend on the polynomial degree of spline basis functions.

A high-order generalized Finite Element Method for multiscale structural dynamics and wave propagation

This paper introduces a high-order multiscale Generalized/eXtended Finite Element Method (GFEM) tailored for the solution of structural dynamics and wave propagation problems exhibiting fine-scale and/or localized solution features such as singularities and discontinuities. The proposed method uses an explicit time-marching scheme together with a block-diagonal lumped mass matrix applicable to arbitrary patch approximation spaces and adopts shape functions computed numerically on-the-fly through the solution of local initial/boundary value problems. Numerical results in one-, two-, and three-dimensions for problems exhibiting singularities in the spatial gradient show that the proposed GFEM is able to capture relevant features of the response using structural-scale meshes that are much coarser than those required by Direct Generalized Finite Element Analyses (DGFEAs) of comparable accuracy. In addition, a detailed study of the critical time step size of the method is presented and compared against the one for DGFEA discretizations that provide similar levels of accuracy. It shows that the proposed method has considerably looser time step size restrictions than available DGFEAs.

Transient thermal analysis of brake disc in c++

This work examines the development of a C++ program to analyze one dimensional thermal conduction for asymmetrical domain using Finite Element Method in C++ environment. Finite Element Library has been used to incorporate various pre-defined Finite Element classes. The various classes have been further grouped into several modules which serve as the fundamental units of any program that encapsulated Finite Element Method. The work strives to illuminate the abstraction behind the computation of the desired unknown solution to a complex real-world problem by capturing the physics of heat conduction and highlights the mathematics required to arrive at an understandable transient solution, which adheres to the conditions of stability, and thus can be relied upon to be implemented in various automotive engineering applications. It also identifies the underlying mechanism through which simple one-dimensional and multi-dimensional partial differential equation problem are generally solved. The elliptic, classic form of partial differential equation has been developed in the form of C++ code. The development begins with elliptic partial differential equations in multi-dimension dimension for steady state heat conduction and then is expanded to included time discretization, which is requirement of a time-dependent problem. Then output of the code is compared to the solution is validated iteratively so as to find the critical time step size for stability, which also validates whether the code actually works for different numerical techniques such as Forward Difference, Backward Difference and the Crank-Nicholson Method, all three distinguished by certain value of a common parameter. The C++ code is both flexible as well as expandable and can be thought of as a mathematical model to solve time-dependent heat conduction problem.

Simulation of discrete cracks driven by nearly incompressible fluid via 2D combined finite‐discrete element method

SummaryIn this work, we propose a novel hydraulic solver in order to simulate key mechanisms that control fluid‐driven cracks in the framework of the combined finite‐discrete element method (FDEM). The main innovative aspect of the present work is the independence of the fluid's critical time step size with the fracture opening. This advantage is extremely important because it means that very fine meshes can be used around areas of interest, such as boreholes, without penalizing the computational cost as fractures propagate (ie, open) and the fluid flows through them. This is a great advantage over other recently introduced approaches that exhibit a dependency of the time step in the form of Δtcrit ∝ (l/a)2 where l is the element size and a is the fracture aperture. This paper presents a series of benchmark cases for the proposed solver. The rationale adopted by the authors was to benchmark and validate the implementation of the hydraulic solver in an incremental fashion, starting from the simplest cases and building in complexity. The results shown in this work clearly demonstrate that the proposed approach is able to reproduce analytical results for fluid flow through a single crack. The results presented in this paper also demonstrate that the new approach is robust enough to deal with complex fracture patterns and complex geometries; the obtained fluid‐driven fracture patterns in the vicinity of a borehole certainly stand to the scrutiny of human visual perception.

Explicit higher-order accurate isogeometric collocation methods for structural dynamics

The objective of the present work is to develop efficient, higher-order space- and time-accurate, methods for structural dynamics. To this end, we present a family of explicit isogeometric collocation methods for structural dynamics that are obtained from predictor–multicorrector schemes. These methods are very similar in structure to explicit finite-difference time-domain methods, and in particular, they exhibit similar levels of computational cost, ease of implementation, and ease of parallelization. However, unlike finite difference methods, they are easily extended to non-trivial geometries of engineering interest. To examine the spectral properties of the explicit isogeometric collocation methods, we first provide a semi-discrete interpretation of the classical predictor–multicorrector method. This allows us to characterize the spatial and modal accuracy of the isogeometric collocation predictor–multicorrector method, irrespective of the considered time-integration scheme, as well as the critical time step size for a particular explicit time-integration scheme. For pure Dirichlet problems, we demonstrate that it is possible to obtain a second-order-in-space scheme with one corrector pass, a fourth-order-in-space scheme with two corrector passes, and a fifth-order-in-space scheme with three corrector passes. For pure Neumann and mixed Dirichlet–Neumann problems, we demonstrate that it is possible to obtain a second-order-in-space scheme with one corrector pass and a third-order-in-space scheme with two corrector passes, and we observe that fourth-order-in-space accuracy may be obtained pre-asymptotically with three corrector passes. We then present second-order-in-time, fourth-order-in-time, and fifth-order-in-time fully discrete predictor–multicorrector algorithms that result from the application of explicit Runge–Kutta methods to the semi-discrete isogeometric collocation predictor–multicorrector method. We confirm the accuracy of the family of explicit isogeometric collocation methods using a suite of numerical examples.

On stability and reflection‐transmission analysis of the bipenalty method in contact‐impact problems: A one‐dimensional, homogeneous case study

SummaryThe stability and reflection‐transmission properties of the bipenalty method are studied in application to explicit finite element analysis of one‐dimensional contact‐impact problems. It is known that the standard penalty method, where an additional stiffness term corresponding to contact boundary conditions is applied, attacks the stability limit of finite element model. Generally, the critical time step size rapidly decreases with increasing penalty stiffness. Recent comprehensive studies have shown that the so‐called bipenalty technique, using mass penalty together with standard stiffness penalty, preserves the critical time step size associated to contact‐free bodies. In this paper, the influence of the penalty ratio (ratio of stiffness and mass penalty parameters) on stability and reflection‐transmission properties in one‐dimensional contact‐impact problems using the same material and mesh size for both domains is studied. The paper closes with numerical examples, which demonstrate the stability and reflection‐transmission behavior of the bipenalty method in one‐dimensional contact‐impact and wave propagation problems of homogeneous materials.

A mixed explicit–implicit time integration approach for nonlinear analysis of base-isolated structures

The paper investigates the accuracy, the stability and the computational efficiency of a mixed explicit–implicit time integration approach proposed for predicting the nonlinear response of base-isolated structures subjected to earthquake excitation. Adopting the central difference method for evaluating the response of the nonlinear base isolation system and the Newmark’s constant average acceleration method for estimating the superstructure linear response, the proposed partitioned solution approach is used to analyze a 3D seismically isolated structure subjected to a bidirectional earthquake excitation. Both isolation systems adopting lead rubber and friction pendulum bearings are considered. Numerical results show that the computational time required by the proposed method, in spite of its conditional stability arising from the use of the central difference method in the explicit integration substep, is clearly reduced in comparison to the widely used implicit time integration method adopted in conjunction with the pseudo-force approach (i.e., pseudo-force method). As a matter of fact, the typical low stiffness of the isolation system leads to a critical time step larger than the one used to define the ground acceleration accurately and the proposed method preserves its computational efficiency even in the case of isolators with very high initial stiffness (i.e., friction pendulum bearings) for which the critical time step size could become smaller.

On the stability of projection methods for the incompressible Navier–Stokes equations based on high-order discontinuous Galerkin discretizations

The present paper deals with the numerical solution of the incompressible Navier–Stokes equations using high-order discontinuous Galerkin (DG) methods for discretization in space. For DG methods applied to the dual splitting projection method, instabilities have recently been reported that occur for small time step sizes. Since the critical time step size depends on the viscosity and the spatial resolution, these instabilities limit the robustness of the Navier–Stokes solver in case of complex engineering applications characterized by coarse spatial resolutions and small viscosities. By means of numerical investigation we give evidence that these instabilities are related to the discontinuous Galerkin formulation of the velocity divergence term and the pressure gradient term that couple velocity and pressure. Integration by parts of these terms with a suitable definition of boundary conditions is required in order to obtain a stable and robust method. Since the intermediate velocity field does not fulfill the boundary conditions prescribed for the velocity, a consistent boundary condition is derived from the convective step of the dual splitting scheme to ensure high-order accuracy with respect to the temporal discretization. This new formulation is stable in the limit of small time steps for both equal-order and mixed-order polynomial approximations. Although the dual splitting scheme itself includes inf–sup stabilizing contributions, we demonstrate that spurious pressure oscillations appear for equal-order polynomials and small time steps highlighting the necessity to consider inf–sup stability explicitly.

An enriched finite element model for wave propagation in fractured media

In this paper a new numerical method has been developed in the context of enriched finite element methods (FEMs) to analyze wave propagation in fractured media. The method combines the advantages of global enrichment with harmonic functions via the Generalized FEM (GFEM) with the efficacy of the Phantom Node Method (PNM), an eXtended FEM (XFEM) variant, to model cracks independently of the mesh. The GFEM enrichment suppresses the spurious oscillation that appear in regular FEM analysis of transient wave propagations due to numerical dispersion and Gibb's phenomenon. For use in explicit simulations, a mass lumping methodology has been introduced with a critical time step size that is both similar to that of the underlining FEM and independent of the location of the fracture. Through three examples, the developed PNM-GFEM is demonstrated to more accurately model wave propagation in fractured media than either the FEM or the PNM/XFEM.

Von Neumann stability analysis of the u–p reproducing kernel formulation for saturated porous media

This paper introduces the von Neumann method to investigate the temporal stability of the displacement---pressure $$(u{-}p)$$(u-p) reproducing kernel formulation for saturated porous media. Both dynamic and quasi-static formulations are considered and the critical time steps are derived. The effect of lumped and consistent matrices on temporal stability is analyzed under explicit temporal discretization. It is shown that lumped matrices have better temporal stability than consistent matrices. The study also shows that nodal support size greatly affects the critical time step size of the formulations. For consistent matrices, larger support size results in smaller critical time step size; however, opposite relation occurs if lumped scheme is used. The numerical study shows that stabilization parameter of the stabilized nodal integration methods reduces the critical time step size. Transient analyses are performed to verify the results from von Neumann analysis.