High-order implicit time integration scheme based on Pad\'e expansions

A single-step high-order implicit time integration scheme for the solution of transient and wave propagation problems is presented. It is constructed from the Pad\'e expansions of the matrix exponential solution of a system of first-order ordinary differential equations formulated in the state-space. A computationally efficient scheme is developed exploiting the techniques of polynomial factorization and partial fractions of rational functions, and by decoupling the solution for the displacement and velocity vectors. An important feature of the novel algorithm is that no direct inversion of the mass matrix is required. From the diagonal Pad\'e expansion of order $M$ a time-stepping scheme of order $2M$ is developed. Here, each elevation of the accuracy by two orders results in an additional system of real or complex sparse equations to be solved. These systems are comparable in complexity to the standard Newmark method, i.e., the effective system matrix is a linear combination of the static stiffness, damping, and mass matrices. It is shown that the second-order scheme is equivalent to Newmark's constant average acceleration method, often also referred to as trapezoidal rule. The proposed time integrator has been implemented in MATLAB using the built-in direct linear equation solvers. In this article, numerical examples featuring nearly one million degrees of freedom are presented. High-accuracy and efficiency in comparison with common second-order time integration schemes are observed. The MATLAB-implementation is available from the authors upon request or from the GitHub repository (to be added).Comment: 43 pages, 19 figure

Automatic 3D modeling by combining SBFEM and transfinite element shape functions

The scaled boundary finite element method (SBFEM) has recently been employed as an efficient means to model three-dimensional structures, in particular when the geometry is provided as a voxel-based image. To this end, an octree decomposition of the computational domain is deployed and each cubic cell is treated as an SBFEM subdomain. The surfaces of each subdomain are discretized in the finite element sense. We improve on this idea by combining the semi-analytical concept of the SBFEM with certain transition elements on the subdomains' surfaces. Thus, we avoid the triangulation of surfaces employed in previous works and consequently reduce the number of surface elements and degrees of freedom. In addition, these discretizations allow coupling elements of arbitrary order such that local p-refinement can be achieved straightforwardly

High-order implicit time integration scheme with controllable numerical dissipation based on mixed-order Pad\'e expansions

A single-step high-order implicit time integration scheme with controllable numerical dissipation at high frequencies is presented for the transient analysis of structural dynamic problems. The amount of numerical dissipation is controlled by a user-specified value of the spectral radius $\rho_\infty$ in the high frequency limit. Using this user-specified parameter as a weight factor, a Pad\'e expansion of the matrix exponential solution of the equation of motion is constructed by mixing the diagonal and sub-diagonal expansions. An efficient timestepping scheme is designed where systems of equations, similar in complexity to the standard Newmark method, are solved recursively. It is shown that the proposed high-order scheme achieves high-frequency dissipation, while minimizing low-frequency dissipation and period errors. The effectiveness of the provided dissipation control and the efficiency of the scheme are demonstrated by numerical examples. A simple guideline for the choice of the controlling parameter and time step size is provided. The source codes written in MATLAB and FORTRAN are available for download at: https://github.com/ChongminSong/HighOrderTimeIntegration.Comment: 37 pages, 36 figures, 89 equation

Three-dimensional image-based numerical homogenisation using octree meshes

The determination of effective material properties of composites based on a three-dimensional representative volume element (RVE) is considered in this paper. The material variation in the RVE is defined based on the colour intensity in each voxel of an image which can be obtained from imaging techniques such as X-ray computed tomography (XCT) scans. The RVE is converted into a numerical model using hierarchical meshing based on octree decompositions. Each octree cell in the mesh is modelled as a scaled boundary polyhedral element, which only requires a surface discretisation on the polyhedron's boundary. The problem of hanging (incompatible) nodes – typically encountered when using the finite element method in conjunction with octree meshes – is circumvented by employing special transition elements. Two different types of boundary conditions (BCs) are used to obtain the homogenised material properties of various samples. The numerical results confirm that periodic BCs provide a better agreement with previously published results. The reason is attributed to the fact that the model based on the periodic BCs is not over-constrained as is the case for uniform displacement BCs

A time-domain approach for the simulation of three-dimensional seismic wave propagation using the scaled boundary finite element method

A direct time-domain approach to simulate seismic wave propagation in three-dimensional unbounded media is proposed based on the Scaled Boundary Finite Element Method (SBFEM). A domain of interest is commonly partitioned into a far field and a near field. The far field is modelled by the semi-analytical SBFEM satisfying rigorously the radiation conditions at infinity. Separate scaled boundary finite elements are employed to reach a balance between computational efficiency and accuracy. The near field is discretized into arbitrarily-shaped scaled boundary finite elements without the occurrence of hanging nodes. This advantage of the SBFEM in mesh generation is leveraged by incorporating the automatic octree-based meshing technique. By exploiting the geometrical similarity of both bounded and unbounded SBFE subdomains the computational cost is reduced. Inspired by the Domain Reduction Method (DRM), seismic waves are introduced to the system via a single layer of elements in the near field. This formulation of seismic input is mathematically convenient as it avoids the direct participation of the formulation of the far field. The proposed approach is attractive in a reliable simulation of the far field, flexible mesh generation of the near field and simple formulation of the seismic excitations. These merits are demonstrated through numerical simulations of seismic wave propagation in a free field and different examples featuring complex geometries in the near fields

Towards higher-order accurate mass lumping in explicit isogeometric analysis for structural dynamics

We present a mass lumping approach based on an isogeometric Petrov-Galerkin method that preserves higher-order spatial accuracy in explicit dynamics calculations irrespective of the polynomial degree of the spline approximation. To discretize the test function space, our method uses an approximate dual basis, whose functions are smooth, have local support and satisfy approximate bi-orthogonality with respect to a trial space of B-splines. The resulting mass matrix is ``close'' to the identity matrix. Specifically, a lumped version of this mass matrix preserves all relevant polynomials when utilized in a Galerkin projection. Consequently, the mass matrix can be lumped (via row-sum lumping) without compromising spatial accuracy in explicit dynamics calculations. We address the imposition of Dirichlet boundary conditions and the preservation of approximate bi-orthogonality under geometric mappings. In addition, we establish a link between the exact dual and approximate dual basis functions via an iterative algorithm that improves the approximate dual basis towards exact bi-orthogonality. We demonstrate the performance of our higher-order accurate mass lumping approach via convergence studies and spectral analyses of discretized beam, plate and shell models

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 $\chi$ 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 $\mathbf{M}_\mathrm{c}$ of cut elements -- especially for high polynomial orders $p$ of the shape functions -- leads to a significant raise in the critical time step size $\Delta t_\mathrm{cr}$. To demonstrate the efficacy 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.Comment: 45 pages, 25 figure

Revisiting Mindlin's theory with regard to a gradient extended phase‐field model for fracture

The application of generalized continuum mechanics is rapidly increasing in different fields of science and engineering. In the literature, there are several theories extending the classical first-order continuum mechanics formulation to include size-effects [1]. One approach is the strain gradient theory with the intrinsic features of regularizing singular stress fields occurring, e.g., near crack tips. It is crucial to realize that using this theory, the strain energy density is still localized around the crack tip, but does not exhibit any signs of a singularity. Therefore, these models seem to be appropriate choices for studying cracks in mechanical problems. Over the past several years, the phase-field method has gathered considerable popularity in the computational mechanics community, in particular in the field of fracture mechanics [2]. Recently, the authors have shown that integrating the strain gradient theory into the phase-field fracture framework is likely to improve the quality of the final results due to the inherent non-singular nature of this theory [3]. In the present work, we will focus on a general formulation of the first strain gradient theory. To this end, the homogenization approach introduced in Ref. [4] is employed. It is based on a series of systematic finite element simulations using different loading cases to determine the equivalent material coefficients on the macro-scale (i.e., for a strain gradient elastic material) by taking the underlying micro-structure into account

A strain gradient enhanced model for the phase‐field approach to fracture

Phase-field modelling has been shown to be a powerful tool for simulating fracture processes and predicting the crack path under complex loading conditions. Note that the total energy of fracture in the classical phase-field formulations includes the strain energy density from the linear elasticity theory resulting in singular stresses at the crack tip. Recently, we have demonstrated that integrating the strain gradient elasticity into the conventional phase-field fracture formulations may improve the final results by alleviating the effects of a singular stress field around the crack tip [1]. The current contribution focuses on a more general formulation of strain gradient elasticity