Temporal and Spatial Coupling Methods for the Efficient Modelling of Dynamic Solids

Abstract

This paper presents efficient coupling methods that accurately reduce the computational cost for modelling solids dynamically with finite elements. A multi-time-step integration algorithm is developed to leverage varying time steps throughout a domain. Interfaces between subdomains are resolved explicitly with the continuity of acceleration and tractions. A spatial coupling method is combined with multiple time steps, allowing for meshes that do not necessarily conform at their interfaces. The method avoids solving additional degrees of freedom at these interfaces, with parameter-free coupling operators defined between meshes. A speedup >12× is achieved in comparison to reference single-time-step methods.

1. Introduction

The temporal and spatial discretisation of structural dynamic problems is directly related to the accuracy and computational cost of the explicit finite element method. Constrained by the Courant–Friedrichs–Lewy (CFL) condition, the critical time step, $\Delta t_C$, is proportional to the element size, and inversely proportional to the dilatational wave speed [1]. This leads to simulations being restricted to $\min \{\Delta t\}$ of the element with the smallest size or highest wave speed. Initially referred to as subcycling, pioneering works allowed for the integration of multiple time steps in a single domain [2,3]. The coupling of various kinematic fields was explored, along with the stability of such algorithms [4,5]. Asynchronous variational integration is a non-trivial alternative, discretising the functional instead of the equations of motion [6]. The main drawback being its high complexity implementation. Heterogeneous asynchronous time integration extended methods to varying, non-integer and large time step ratios [7,8]. However, maintaining the continuity of kinematics at the interface remains a challenge, especially for all three fields [9]. In recent works, energy-conserving methods have been developed; however, they depart far from the CFL condition to resolve the interface conditions [10]. Spatially, non-matching mesh algorithms facilitate more flexible geometric modelling [11]. Nitsche’s method has shown to weakly enforce conditions on the interfaces of non-matching meshes without additional unknowns, but commonly suffers from sensitivity to parameters [12,13,14,15]. Following similar principles of weak continuity, Lagrange multipliers are called upon in the mortar-based methods [16,17,18,19,20]. These types of methods have proven to be very robust, however they struggle to fulfil the inf-sup stability condition, and incur a large computational cost with mapping master and slave nodes [21]. In comparison, the use of localised Lagrange multipliers introduces a frame that independently enforces compatibility constraints on each boundary [22,23,24,25,26,27]. However, as is common with other methods that introduce an additional discretised interface, the construction of this interface is neither trivial nor computationally cheap [28,29,30,31]. This brief review justifies the need for more efficient couplings, both temporally and spatially. Coupling algorithms that allow for varying time-step ratios, whilst stepping close to the CFL condition, concertedly those that solve for non-matching meshes, without increasing the degrees of freedom, remain a hot topic of research.

2. Governing Equations of Dynamic Solids

The problem of a solid body subject to impact is described through the partitioning of a domain as illustrated in Figure 1. The deformation is governed by the momentum balance equation acting on the solid domain $\Omega$:

$$\begin{array}{c}\hfill \rho \ddot{\mathbf{u}}=\nabla ·\sigma +\rho \mathbf{b},\mathrm{in}\Omega \times [0,T]\end{array}$$

(1)

where $\rho$, $\ddot{\mathbf{u}}$, $\sigma$ and $\mathbf{b}$ denote the density, acceleration field, Cauchy stress tensor and body forces. Deformation is described at time $t\in [0,T]$ for a specified constant $T>0$. Updated Lagrangian formulations of a single body are found in the following [32,33,34]; here, we extend this to a partitioned formulation to solve multiple solid subdomains.

Defining multiple subdomains, we state $\mathcal{B}$ is a solid body in an open region $\Omega \subset {\mathbb{R}}^3$, with its boundary denoted $\partial \Omega$. $\Omega$ is partitioned into $\mathcal{S}$ non-overlapping subdomains:

$$\begin{array}{c}\hfill \Omega =\bigcup _{i=1}^S{\Omega }_i\mathrm{and}{\Omega }_i\cap {\Omega }_j=Ø\mathrm{for}i\ne j\end{array}$$

(2)

Starting from $\mathcal{S}=2$, for a two-subdomain partitioning, where ${\Omega }_L$ and ${\Omega }_s$ denote large and small subdomains. On the boundary of the two matching subdomains, $\Gamma$, we enforce:

$$\begin{array}{c}\hfill {\ddot{\mathbf{u}}}_L={\ddot{\mathbf{u}}}_s\mathrm{on}\Gamma \end{array}$$

(3)

$$\begin{array}{c}\hfill {\mathbf{t}}_L={\mathbf{t}}_s\mathrm{on}\Gamma \end{array}$$

(4)

to ensure the continuity of acceleration $\ddot{\mathbf{u}}$ and tractions $\mathbf{t}$, as well as enforcing Dirichlet and Neumann boundary conditions. The variational formulation of the dynamic equilibrium in Equation (1) can be described for both ${\Omega }_L$ and ${\Omega }_s$ as the following:

$$\begin{array}{c}\hfill \underset{{\Omega }_L}{\int }{\rho }_L{\ddot{\mathbf{u}}}_L·\delta \dot{\mathbf{u}}d\Omega =\underset{{\Gamma }_L}{\int }{\mathbf{t}}_L·\delta \dot{\mathbf{u}}d\Gamma -\underset{{\Omega }_L}{\int }{\sigma }_L:\delta {\mathbf{D}}_{L}d\Omega +\underset{{\Omega }_L}{\int }{\rho }_L{\mathbf{b}}_L·\delta \dot{\mathbf{u}}d\Omega \end{array}$$

(5)

$$\begin{array}{c}\hfill \underset{{\Omega }_s}{\int }{\rho }_s{\ddot{\mathbf{u}}}_s·\delta \dot{\mathbf{u}}d\Omega =\underset{{\Gamma }_s}{\int }{\mathbf{t}}_s·\delta \dot{\mathbf{u}}d\Gamma -\underset{{\Omega }_s}{\int }{\sigma }_s:\delta {\mathbf{D}}_{s}d\Omega +\underset{{\Omega }_s}{\int }{\rho }_s{\mathbf{b}}_s·\delta \dot{\mathbf{u}}d\Omega \end{array}$$

(6)

where we denote a variational velocity $\delta \dot{\mathbf{u}}\in {\mathcal{V}}_0$, in a space ${\mathcal{V}}_0$ where $\forall \delta \dot{\mathbf{u}}\in H^1\left({\Omega }_i\right)$, and $\mathbf{D}$ as the rate of deformation. A discrete approximation of the variational form can be reduced to the following ordinary differential equations:

$$\begin{array}{c}\hfill {\mathbf{M}}_L{\ddot{\mathbf{u}}}_L={\mathbf{f}}_L^{\mathrm{ext}}-{\mathbf{f}}_L^{\mathrm{int}};{\mathbf{M}}_s{\ddot{\mathbf{u}}}_s={\mathbf{f}}_s^{\mathrm{ext}}-{\mathbf{f}}_s^{\mathrm{int}}\end{array}$$

(7)

where for subdomains ${\Omega }_L$ and ${\Omega }_s$, we sum, over a number of finite elements, ${\mathcal{N}}_L$ and ${\mathcal{N}}_s$, in each subdomain:

$$\begin{array}{c}\hfill {\mathbf{M}}_L=\sum _{e=1}^{{\mathcal{N}}_L}{\int }_{{\Omega }_e}{\rho }_L{\mathbf{N}}^{\mathrm{T}}\mathbf{N}d{\Omega }_e;{\mathbf{M}}_s=\sum _{e=1}^{{\mathcal{N}}_s}{\int }_{{\Omega }_e}{\rho }_s{\mathbf{N}}^{\mathrm{T}}\mathbf{N}d{\Omega }_e\end{array}$$

(8)

$$\begin{array}{c}\hfill {\mathbf{f}}_L^{\mathrm{ext}}=\sum _{e=1}^{{\mathcal{N}}_L}{\int }_{{\Omega }_e}{\mathbf{N}}^{\mathrm{T}}{\mathbf{t}}_{L}d{\Omega }_e;{\mathbf{f}}_s^{\mathrm{ext}}=\sum _{e=1}^{{\mathcal{N}}_s}{\int }_{{\Omega }_e}{\mathbf{N}}^{\mathrm{T}}{\mathbf{t}}_{s}d{\Omega }_e\end{array}$$

(9)

$$\begin{array}{c}\hfill {\mathbf{f}}_L^{\mathrm{int}}=\sum _{e=1}^{{\mathcal{N}}_L}{\int }_{{\Omega }_e}{\mathbf{B}}^{\mathrm{T}}{\sigma }_{L}d{\Omega }_e{\mathbf{f}}_s^{\mathrm{int}}=\sum _{e=1}^{{\mathcal{N}}_s}{\int }_{{\Omega }_e}{\mathbf{B}}^{\mathrm{T}}{\sigma }_{s}d{\Omega }_e\end{array}$$

(10)

The Lagrangian shape functions are represented by $\mathbf{N}$, with their derivatives denoted as $\mathbf{B}$. As is common within the explicit finite element method, $\mathbf{M}$ is lumped for each subdomain, with ${\mathbf{f}}^{\mathrm{ext}}$ and ${\mathbf{f}}^{\mathrm{int}}$ computed with vectors too. We elect to use the leapfrog time integration scheme to step through time, staggering the solution of each kinematic quantity such that:

$$\begin{array}{c}\hfill {\ddot{\mathbf{u}}}_L^n={\mathbf{M}}_L^{-1}({\mathbf{f}}_L^{\mathrm{ext}}-{\mathbf{f}}_L^{\mathrm{int}});{\ddot{\mathbf{u}}}_s^n={\mathbf{M}}_s^{-1}({\mathbf{f}}_s^{\mathrm{ext}}-{\mathbf{f}}_s^{\mathrm{int}})\end{array}$$

(11)

$$\begin{array}{c}\hfill {\dot{\mathbf{u}}}_L^{n+1/2}={\dot{\mathbf{u}}}_L^{n-1/2}+{\ddot{\mathbf{u}}}_L^n·\Delta t_L;{\dot{\mathbf{u}}}_s^{n+1/2}={\dot{\mathbf{u}}}_s^{n-1/2}+{\ddot{\mathbf{u}}}_s^n·\Delta t_s\end{array}$$

(12)

$$\begin{array}{c}\hfill {\mathbf{u}}_L^{n+1}={\mathbf{u}}_L^n+{\dot{\mathbf{u}}}_L^{n+1/2}·\Delta t_L;{\mathbf{u}}_s^{n+1}={\mathbf{u}}_s^n+{\dot{\mathbf{u}}}_s^{n+1/2}·\Delta t_s\end{array}$$

(13)

noting that a diagonal mass matrix allows for a direct computation of acceleration. Velocity $\dot{\mathbf{u}}$ is computed on the half time step, with displacement $\mathbf{u}$ found for each full time step. Next, we summarise the temporal coupling, enabled by multi-time-step (MTS) integration.

3. Multi-Time-Step Integration

Multi-time stepping enables partitioned subdomains ${\Omega }_L$ and ${\Omega }_s$ to integrate with $\Delta t_L$ and $\Delta t_s$, respectively. However, to allow for this difference, special attention must be given to the solution of the interface $\Gamma$. Crucially, the conditional stability of explicit methods requires an element’s time step to obey the CFL condition for a linear undamped system as follows:

$$\Delta t_C={\displaystyle \frac{2}{{\omega }_C}}\le \underset{e}{min}\left({\displaystyle \frac{h_e}{c_e}}\right)$$

(14)

Here, we represent $\Delta t_C$ as the critical time step, ${\omega }_C$ as the maximum eigenfrequency, $h_e$ as the characteristic length of an element e and $c_e$ the dilatational (longitudinal) wave speed.

3.1. Salient Multi-Time-Stepping Features

The asynchronous integration is enabled with three important groups of computations. The first is the explicit computation of the acceleration on the interface $\Gamma$:

$$\begin{array}{c}\hfill {\mathbf{M}}_{\Gamma }=\left({\mathbf{C}}_s^{\mathrm{T}}{\mathbf{M}}_s{\mathbf{C}}_s\right)+\left({\mathbf{C}}_L^{\mathrm{T}}{\mathbf{M}}_L{\mathbf{C}}_L\right)\end{array}$$

(15)

$$\begin{array}{c}\hfill {\mathbf{f}}_{\Gamma }^{\mathrm{int}}={\mathbf{C}}_s^{\mathrm{T}}{\mathbf{f}}_s^{\mathrm{int}}+{\mathbf{C}}_L^{\mathrm{T}}{\mathbf{f}}_L^{\mathrm{int}};{\mathbf{f}}_{\Gamma }^{\mathrm{ext}}={\mathbf{C}}_s^{\mathrm{T}}{\mathbf{f}}_s^{\mathrm{ext}}+{\mathbf{C}}_L^{\mathrm{T}}{\mathbf{f}}_L^{\mathrm{ext}}\end{array}$$

(16)

We define indicator matrices (vectors in 1-D) for each subdomain $\mathbf{C}$ to identify the degrees of freedom on the interface of subdomains with dimensions ${\mathbf{C}}_i\in {\mathbb{R}}^{{\mathcal{N}}_{n_i}\times {\mathcal{N}}_{\Gamma }}$ for nodal number $\mathcal{N}$. These summations provide the ingredients for computing the interface acceleration:

$$\begin{array}{c}\hfill {\ddot{\mathbf{u}}}_{\Gamma }={\mathbf{M}}_{\Gamma }^{-1}({\mathbf{f}}_{\Gamma }^{\mathrm{ext}}-{\mathbf{f}}_{\Gamma }^{\mathrm{int}})\end{array}$$

(17)

where we compute ${\ddot{\mathbf{u}}}_{\Gamma }$ at each large time step $\Delta t_L$. The integration of the subdomains is controlled by the definition of the time-step ratios. Suppose the two subdomains begin at a similar point in time $t_L^N=t_s^n$, where N and n are the small and large steps, respectively. The time after the maximum stable integration step (governed by the CFL condition) on each subdomain is referred to as the trial time $t_T$, such that:

$$\begin{array}{c}\hfill t_{TL}^{N+1}=t_L+\Delta t_{CL};t_{Ts}^{n+k}=t_s+\Delta t_{Cs}\end{array}$$

(18)

where for every $\Delta t_L$, k small time steps elapse since the last point in time where subdomains are equal in time. Now we can define the current time-step ratio, $t_{ratio}^{n+k}$, and next time-step ratio, $t_{ratio}^{n+k+1}$ for the advancement of ${\Omega }_s$ with:

$$\begin{array}{c}\hfill t_{ratio}^{n+k}={\displaystyle \frac{t_s^{n+k}-t_L^N}{t_{TL}^{N+1}-t_L^N}};t_{ratio}^{n+k+1}={\displaystyle \frac{(t_s^{n+k}+\Delta t_s^{n+k})-t_L^N}{t_{TL}^{N+1}-t_L^N}}\end{array}$$

(19)

Starting from each common time step with the integration of ${\Omega }_s$, the number of small time steps $\Delta t_s$ is determined by the evaluation of time-step ratios $t_{ratio}^{n+k}$ and $t_{ratio}^{n+k+1}$. If the condition of $t_{ratio}^{n+k+1}\le 1$ or ($t_{ratio}^{n+k}\le 1$ and $t_{ratio}^{n+k+1}\le 1$) is satisfied, further steps on $\Delta t_s$ can proceed before integrating ${\Omega }_L$ over $\Delta t_L$. As a consequence of subdomains integrating over their own respective time step, we ensure that the proposed method still finds a common time between all subdomains. Following the small trial time $t_{Ts}$ exceeding the large trial time $t_{TL}$, the method computes two additional ratios ${\alpha }_L$ and ${\alpha }_s$. These denote the reduction factors required to maintain the subdomains in synchronisation where $0\le \alpha \le 1$. Hence:

$$\begin{array}{c}\hfill {\alpha }_L=1-{\displaystyle \frac{(t_{TL}^N-t_s^{n+k})}{(t_{TL}^{N+1}-t_L^N)}};{\alpha }_s=1-{\displaystyle \frac{(t_{Ts}^{n+k}-t_{TL}^{N+1})}{(t_{Ts}^{n+k+1}-t_s^{n+k})}}\end{array}$$

(20)

Through computing $\alpha$ on both subdomains, we can compare reduction factors, such that:

$$\begin{array}{c}\hfill t_L^{N+1}=t_s^{n+k}=\left\{\begin{array}{c}{t}_s^n+\sum _{k=0}^{k-1}\Delta t_{Cs}^{n+k}+{\alpha }_s·\Delta t_{Cs}^{n+k},{\alpha }_s>{\alpha }_L\hfill \\ t_s^n+{\alpha }_L·\Delta t_{CL}^N,{\alpha }_L\ge {\alpha }_s\hfill \end{array}\right.\end{array}$$

(21)

where we look to reduce the time step of the subdomain closest to the CFL condition. In the following section, we summarise the multi-time-stepping algorithm, with each of these key features, as well as its implementation.

3.2. Summary of Temporal Algorithm

We provide an overview of the method required to integrate ${\Omega }_L$ and ${\Omega }_s$ with $\Delta t_L$ and $\Delta t_s$. It shows a single large step N, exemplifying each of the features mentioned above. Note that $TOL=1\times {10}^{-6}$ is used for when $\Delta t_L\approx \Delta t_s$, and computation of ${\mathbf{M}}_{\Gamma }$ in Equation (15) is only required at $t=0$ for a mass-conserving problem. One full loop of the procedure times $t_L=t_s$ is followed, with the subdomain synchronised after each large step N. The proposed method is implemented in an open-source python code, found in the following repository: https://github.com/kinfungchan/multi-time-step-integration (accessed on 9 December 2024) [35]. It contains re-implementations of methods from the literature for the two-subdomain cases [9,10]. Whilst we depict the case of just two subdomains, the Algorithm 1 can be extended to multiple by processing subdomains as pairs.

Algorithm 1 Summary of Algorithm for Coupling in Time from N to $N+1$

1: procedure a two-subdomain multi-time integration step   2:       while $t_{ratio}^{n+k+1}\le 1$ or ($t_{ratio}^{n+k}\le 1$ and $t_{ratio}^{n+k+1}\le 1+\mathrm{TOL}$) do   3:          Integrate subdomain ${\Omega }_s$ with ${\ddot{\mathbf{u}}}_{\Gamma }$ and compute force vectors ${\mathbf{f}}_s^{\mathrm{int}},{\mathbf{f}}_s^{\mathrm{ext}}$ over $\Delta t_s$   4:          Compute trial times $t_{Ts}^{n+k}$, $t_{TL}^{N+1}$ and time step ratios $t_{ratio}^{n+k}$, $t_{ratio}^{n+k+1}$   5:       Compute time step reduction factors ${\alpha }_L,{\alpha }_s$   6:       if ${\alpha }_L\ge {\alpha }_s$ then   7:          $\Delta t_L={\alpha }_L·\Delta t_L$   8:       else   9:          $\Delta t_s={\alpha }_s·\Delta t_s$ 10:          Integrate subdomain ${\Omega }_s$ with ${\ddot{\mathbf{u}}}_{\Gamma }$ and recompute force vectors ${\mathbf{f}}_s^{\mathrm{int}},{\mathbf{f}}_s^{\mathrm{ext}}$ over $\Delta t_s$ 11:          Recompute trial times $t_{Ts}^{n+k}$, $t_{TL}^{N+1}$ and time step ratios $t_{ratio}^{n+k}$, $t_{ratio}^{n+k+1}$ 12:       Integrate subdomain ${\Omega }_L$ with ${\ddot{\mathbf{u}}}_{\Gamma }$ and compute force vectors ${\mathbf{f}}_L^{\mathrm{int}},{\mathbf{f}}_L^{\mathrm{ext}}$ over $\Delta t_L$ 13:       Compute interface acceleration ${\ddot{\mathbf{u}}}_{\Gamma }$ with Equations (15)–(17) for next time step 14:       Recompute trial times $t_{Ts}^{n+k}$, $t_{TL}^{N+1}$ and time step ratios $t_{ratio}^{n+k}$, $t_{ratio}^{n+k+1}$

3.3. Numerical Examples in Time

We present a numerical example in 1-D, with the elastic wave propagation through a heterogeneous bar. Suppose the domain $\Omega$ is split into two subdomains ${\Omega }_L$ and ${\Omega }_s$ of similar discretisation, with isotropic elastic properties of $E_L=207$ GPa, $E_s=1000$ GPa and ${\rho }_L={\rho }_s=7.83\times {10}^{-6}$ kgmm⁻³. These material properties result in a non-integer $m=2.19$, where the time-step ratio is solely driven by the dissimilar material properties. Figure 2 depicts the bar configuration. The velocity boundary condition is applied to ${\Omega }_L$ at $x=0$ with $\dot{u}\left(t\right)=0.01\sin \left(2\pi {\omega }_{L}t\right)$ ms⁻¹, where we define a half sine wave with a frequency of ${\omega }_L={\left(125((7.83\times {10}^{-6})/207)\right)}^{-1}$ rads⁻¹. The difference in material impedance results in a portion of the incident wave being transmitted and the remainder reflected in the opposite direction.

In Figure 3, a comparison between the coupled solution (${\Omega }_L$ with ${\Omega }_s$) and the monolithic (single-time-step) solution is presented at four separate time stamps. The multi-time-step solution solves ${\dot{\mathbf{u}}}_L$ and ${\dot{\mathbf{u}}}_s$ over $\Delta t_L$ and $\Delta t_s$, whereas ${\dot{\mathbf{u}}}_{mono}$ is limited by $\Delta t_s$. Consequently, for $m=2.19$, our method reduces the number of integration steps on ${\Omega }_L$ by half. From prescription of the full wave at $t1$, through to the transmission and reflection of the wave at $t4$, the MTS solution aligns very well with the single-time-step solution, despite halving the number of large time steps. This reduction in computational effort is even more prominent for highly heterogeneous configurations, as well as variance in spatial discretisation.

For the above simulation, we also compute the energy components, at each small time step, of each subdomain with the following:

$$\begin{array}{c}\hfill W_{\mathrm{ext}}^{n+1}=W_{\mathrm{ext}}^n+{\displaystyle \frac{\Delta t_i^{n+1/2}}{2}}{\left({\dot{\mathbf{u}}}_i^{n+1/2}\right)}^{\mathrm{T}}({\mathbf{f}}_{\mathrm{ext}}^n+{\mathbf{f}}_{\mathrm{ext}}^{n+1})=W_{\mathrm{ext}}^n+{\displaystyle \frac{1}{2}}\Delta {\mathbf{u}}_i^{\mathrm{T}}({\mathbf{f}}_{\mathrm{ext}}^n+{\mathbf{f}}_{\mathrm{ext}}^{n+1})\end{array}$$

(22)

$$\begin{array}{c}\hfill W_{\mathrm{int}}^{n+1}=W_{\mathrm{int}}^n+{\displaystyle \frac{\Delta t_i^{n+1/2}}{2}}{\left({\dot{\mathbf{u}}}_i^{n+1/2}\right)}^{\mathrm{T}}({\mathbf{f}}_{\mathrm{int}}^n+{\mathbf{f}}_{\mathrm{int}}^{n+1})=W_{\mathrm{int}}^n+{\displaystyle \frac{1}{2}}\Delta {\mathbf{u}}_i^{\mathrm{T}}({\mathbf{f}}_{\mathrm{int}}^n+{\mathbf{f}}_{\mathrm{int}}^{n+1})\end{array}$$

(23)

$$\begin{array}{c}\hfill W_{\mathrm{kin}}^n={\displaystyle \frac{1}{2}}{\left({\dot{\mathbf{u}}}_i^{n+1/2}\right)}^{\mathrm{T}}\mathbf{M}{\dot{\mathbf{u}}}_i^{n+1/2}\end{array}$$

(24)

where n can be interchanged with N when evaluating ${\Omega }_L$. The balance of energy can be evaluated in a similar way to the works of Neal and Belytschko [3], with the following:

$$\begin{array}{c}\hfill |W_{\mathrm{ext}}-W_{\mathrm{int}}+W_{\mathrm{kin}}|\le ||W_{bal}\left|\right|\end{array}$$

(25)

In Figure 4, we show each of the components of energy and its overall balance. For both monolithic and multi-time-step solutions, a smooth transition of energy is observed as the wave interacts with $\Gamma$. Remarkably, as the temporal coupling is enforced, the $W_{bal}$ for both ${\Omega }_L$ and ${\Omega }_s$ is of the order $1\times {10}^{-13}$ kNmm, significantly below each of the components at $1\times {10}^{-8}$ kNmm. This numerical example captures the propagation of a smooth pulse; however, severe loading cases, highly heterogeneous domains, and 3-D problems, with further details are provided in the following work [35].

4. Solving Non-Matching Meshes

The problem of non-matching meshes is commonly found when simulating the dynamical behaviour of solids. We present an algorithm, combined with multi-time stepping, that relaxes the constraint of these conforming nodes, allowing for a coarser representation of a subdomain to be utilised, hence reducing the computational overhead.

4.1. Combined Spatial and Temporal Coupling

The following section follows on from the governing equations defined in Section 2; however, we now allow for the non-overlapping interface $\Gamma ={\Gamma }_L\cap {\Gamma }_s$ to consist of two non-matching spatial discretisations. Their compatibility is maintained such that:

$$\begin{array}{c}\hfill {\ddot{\mathbf{u}}}_{\Gamma L}\left({\mathbf{C}}_L{\mathbf{x}}_L\right)={\ddot{\mathbf{u}}}_{\Gamma s}\left({\mathbf{C}}_s{\mathbf{x}}_L\right)={\ddot{\mathbf{u}}}_{\Gamma }\left({\mathbf{x}}_{\Gamma }\right)\end{array}$$

(26)

$$\begin{array}{c}\hfill {\mathbf{t}}_{\Gamma L}\left({\mathbf{C}}_L{\mathbf{x}}_L\right)={\mathbf{t}}_{\Gamma s}\left({\mathbf{C}}_s{\mathbf{x}}_s\right)={\mathbf{t}}_{\Gamma }\left({\mathbf{x}}_{\Gamma }\right)\end{array}$$

(27)

where positions use an incidence $\mathbf{C}$ matrix, such that ${\mathbf{C}}_L{\mathbf{x}}_L\in {\Gamma }_L$ and ${\mathbf{C}}_s{\mathbf{x}}_s\in {\Gamma }_s$ for each subdomain, and the externally assembled interface contains ${\mathbf{x}}_{\Gamma }\in \Gamma$ to describe the common boundary between the subdomains. The total virtual power can be summated for two subdomains to give $\delta P=\delta P_L+\delta P_s+\delta P_{\Gamma }$, where a general form is obtained:

$$\begin{array}{cc}\delta P\hfill & =\delta {\dot{\mathbf{u}}}_L^{\mathrm{T}}\{{\mathbf{f}}_L^{\mathrm{int}}-{\mathbf{f}}_L^{\mathrm{ext}}+{\mathbf{M}}_L{\ddot{\mathbf{u}}}_L+{\mathbf{N}}_L^{\mathrm{T}}{\mathbf{f}}_{\Gamma L}\}+\delta {\dot{\mathbf{u}}}_s^{\mathrm{T}}\{{\mathbf{f}}_s^{\mathrm{int}}-{\mathbf{f}}_s^{\mathrm{ext}}+{\mathbf{M}}_s{\ddot{\mathbf{u}}}_s+{\mathbf{N}}_s^{\mathrm{T}}{\mathbf{f}}_{\Gamma s}\}\hfill \\ \hfill & +\delta {\mathbf{f}}_{\Gamma L}^{\mathrm{T}}\{{\mathbf{N}}_L^{\mathrm{T}}{\dot{\mathbf{u}}}_{\Gamma L}-{\mathbf{L}}_L^{\mathrm{T}}{\dot{\mathbf{u}}}_{\Gamma }\}+\delta {\mathbf{f}}_{\Gamma s}^{\mathrm{T}}\{{\mathbf{N}}_s^{\mathrm{T}}{\dot{\mathbf{u}}}_{\Gamma s}-{\mathbf{L}}_s^{\mathrm{T}}{\dot{\mathbf{u}}}_{\Gamma }\}-\delta {\dot{\mathbf{u}}}_{\Gamma }^{\mathrm{T}}\{{\mathbf{N}}_L^{\mathrm{T}}{\mathbf{f}}_{\Gamma L}+{\mathbf{N}}_s^{\mathrm{T}}{\mathbf{f}}_{\Gamma s}\}\hfill \end{array}$$

(28)

where variations in velocity coupling force on the interface $\dot{\mathbf{u}}$ and ${\mathbf{f}}_{\Gamma i}$ are accounted for. ${\mathbf{N}}_i^{\mathrm{T}}$ and ${\mathbf{L}}_i^{\mathrm{T}}$ are interpolation (or prolongation) and incidence operators ($\Gamma$ to ${\Omega }_i$), respectively. To map to the interface, we describe this spatial coupling operator ${\mathbf{N}}_i$ in more detail:

$$\begin{array}{c}\hfill {\mathbf{N}}_i\in {\mathbb{R}}^{{\mathcal{N}}_{\Gamma i}\times {\mathcal{N}}_{\Gamma }}\end{array}$$

(29)

with dimensions determined by ${\mathcal{N}}_{\Gamma i}$ and ${\mathcal{N}}_{\Gamma }$, as the number of nodes on the interface of the subdomain ${\Omega }_i$ and the number of nodes on the interface $\Gamma$, respectively. Therefore, ${\mathbf{N}}_{\Gamma }$ interpolates using Lagrangian shape functions for the two subdomains ${\Omega }_L$ and ${\Omega }_s$:

$$\begin{array}{c}\hfill {\mathbf{N}}_{\Gamma }\left(\mathbf{x}\right)=\delta ({\mathbf{x}}_{\Gamma }-{\mathbf{x}}_{\Gamma i}),\mathrm{for}i=L,s\end{array}$$

(30)

where $\delta$ is viewed as a dirac Delta function for coincident nodes. It is convenient to define a restriction operator ${\mathbf{R}}_i\in {\mathbb{R}}^{{\mathcal{N}}_{\Gamma }\times {\mathcal{N}}_{\Gamma i}}$ as the transpose of ${\mathbf{N}}_i$, to map both forces and mass from subdomain interfaces ${\Gamma }_L$ and ${\Gamma }_s$ onto $\Gamma$. The summation on the interfaces now becomes:

$$\begin{array}{c}\hfill {\mathbf{M}}_{\Gamma }={\mathbf{R}}_L{\mathbf{C}}_L^{\mathrm{T}}{\mathbf{M}}_L+{\mathbf{R}}_s{\mathbf{C}}_s^{\mathrm{T}}{\mathbf{M}}_s={\mathbf{R}}_L{\mathbf{M}}_{\Gamma L}+{\mathbf{R}}_s{\mathbf{M}}_{\Gamma s}\end{array}$$

(31)

$$\begin{array}{c}\hfill {\mathbf{f}}_{\Gamma }^{\mathrm{int}}={\mathbf{R}}_L{\mathbf{C}}_L^{\mathrm{T}}{\mathbf{f}}_L^{\mathrm{int}}+{\mathbf{R}}_s{\mathbf{C}}_s^{\mathrm{T}}{\mathbf{f}}_s^{\mathrm{int}}={\mathbf{R}}_L{\mathbf{f}}_{\Gamma L}^{\mathrm{int}}+{\mathbf{R}}_s{\mathbf{f}}_{\Gamma s}^{\mathrm{int}}\end{array}$$

(32)

$$\begin{array}{c}\hfill {\mathbf{f}}_{\Gamma }^{\mathrm{ext}}={\mathbf{R}}_L{\mathbf{C}}_L^{\mathrm{T}}{\mathbf{f}}_L^{\mathrm{ext}}+{\mathbf{R}}_s{\mathbf{C}}_s^{\mathrm{T}}{\mathbf{f}}_s^{\mathrm{ext}}={\mathbf{R}}_L{\mathbf{f}}_{\Gamma L}^{\mathrm{ext}}+{\mathbf{R}}_s{\mathbf{f}}_{\Gamma s}^{\mathrm{ext}}\end{array}$$

(33)

where we compute mass ${\mathbf{M}}_{\Gamma }$, internal force ${\mathbf{f}}_{\Gamma }^{\mathrm{int}}$ and external force ${\mathbf{f}}_{\Gamma }^{\mathrm{ext}}$ to allow for the explicit computation of ${\ddot{\mathbf{u}}}_{\Gamma }$ in Equation (17) to be recalled. These operators are analogous to concepts in multigrid methods and localised Lagrange multipliers (LLMs) [36,37]. Subsequently, we map ${\ddot{\mathbf{u}}}_{\Gamma }$ from $\Gamma$, back to the subdomains’ interfaces:

$$\begin{array}{c}\hfill {\ddot{\mathbf{u}}}_{\Gamma L}={\mathbf{N}}_L{\ddot{\mathbf{u}}}_{\Gamma };{\ddot{\mathbf{u}}}_{\Gamma s}={\mathbf{N}}_s{\ddot{\mathbf{u}}}_{\Gamma }\end{array}$$

(34)

We illustrate a non-matching mesh in Figure 5 and compute its operators through exemplifying a linear isoparametric mapping in 2-D, where $\Gamma$ is discretised with line elements to depict the simplicity of this coupling.

For the interpolation matrix, we elucidate that from $\Gamma$ to ${\Gamma }_L$ the mapping is simply one-to-one and ${\mathbf{N}}_L$ will always take the form of an identity matrix ${\Omega }_L$, whereas ${\mathbf{N}}_s$ requires computation of the shape functions:

$$\begin{array}{c}\hfill {\mathbf{N}}_L=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right];{\mathbf{N}}_s=\left[\begin{array}{ccc}1& 0& 0\\ 1/3& 2/3& 0\\ 0& 2/3& 1/3\\ 0& 0& 1\end{array}\right];{\mathbf{R}}_L={\mathbf{N}}_L^{\mathrm{T}};{\mathbf{R}}_s={\mathbf{N}}_s^{\mathrm{T}}\end{array}$$

(35)

where, for this example case, it can also be shown that the restriction is the transpose of the interpolation. The interface $\Gamma$ is assumed to share the same geometrical description on both subdomain interfaces ${\Gamma }_L$ and ${\Gamma }_s$, without overlap or separation. The spatial and temporal methods are combined to give the following Algorithm 2:

Algorithm 2 Summary of Non-Matching Mesh Algorithm with Multi-Time Stepping

1: procedureintegrate a two-domain non-matching mesh with MTS   2:       while $t_{ratio}^{n+k+1}\le 1$ or ($t_{ratio}^{n+k}\le 1$ and $t_{ratio}^{n+k+1}\le 1+TOL$) do   3:          Compute ${\ddot{\mathbf{u}}}_{{\Gamma }_s}$ with operator ${\mathbf{N}}_s$ in Equation (34) on ${\Gamma }_s$   4:          Integrate small domain ${\Omega }_s$ and compute vectors ${\mathbf{f}}_s^{\mathrm{int}},{\mathbf{f}}_s^{\mathrm{ext}}$   5:          Compute trial times $t_{Ts}^{n+k}$, $t_{TL}^{N+1}$ and time step ratios $t_{ratio}^{n+k}$, $t_{ratio}^{n+k+1}$   6:       Compute time step reduction factors ${\alpha }_L,{\alpha }_s$   7:       if ${\alpha }_L\ge {\alpha }_s$ then   8:          $\Delta t_L={\alpha }_L·\Delta t_L$   9:       else 10:          $\Delta t_s={\alpha }_s·\Delta t_s$ 11:          Compute ${\ddot{\mathbf{u}}}_{{\Gamma }_s}$ with operator ${\mathbf{N}}_s$ in Equation (34) on ${\Gamma }_s$ 12:          Integrate small domain ${\Omega }_s$ and compute vectors ${\mathbf{f}}_s^{\mathrm{int}},{\mathbf{f}}_s^{\mathrm{ext}}$ 13:          Recompute trial times $t_{Ts}^{n+k}$, $t_{TL}^{N+1}$ and time step ratios $t_{ratio}^{n+k}$, $t_{ratio}^{n+k+1}$ 14:       Compute ${\ddot{\mathbf{u}}}_{{\Gamma }_L}$ with operator ${\mathbf{N}}_L$ in Equation (34) on ${\Gamma }_L$ 15:       Integrate large domain ${\Omega }_L$ and compute vectors ${\mathbf{f}}_L^{\mathrm{int}},{\mathbf{f}}_L^{\mathrm{ext}}$ 16:       Summate kinetics with ${\mathbf{R}}_L,{\mathbf{R}}_s,{\mathbf{C}}_L,{\mathbf{C}}_s$ to find ${\mathbf{M}}_{\Gamma },{\mathbf{f}}_{\Gamma }^{\mathrm{int}}$ and ${\mathbf{f}}_{\Gamma }^{\mathrm{ext}}$ with Equation (31)–(33) 17:       Compute trial times $t_{Ts}^{n+k}$, $t_{TL}^{N+1}$ and time step ratios $t_{ratio}^{n+k}$, $t_{ratio}^{n+k+1}$

4.2. Numerical Examples in Space and Time

The following numerical study looks to represent the stress gradients prior to fracture in a compact-tension (CT) specimen test, utilising a similar geometry to the literature [38,39]. Figure 6a portrays the geometry modelled in the following example. As the specimen is loaded, stress concentrates about the specimen’s crack tip. We apply a ramped velocity $\dot{\mathbf{u}}\left(t\right)$ boundary condition on nodes that create a semicircle for upper and lower pins, as shown in Figure 6b, with a maximum magnitude of 0.2 ms⁻¹. Whilst uniformly distributed velocities are applied to each of the nodes in the pins, to replicate the contact pressure on the pins, methods such as those applied by Triclot et al. should be considered [40]. Material properties are similar to alumina with $E=370$ GPa, $\rho =3.9\times {10}^{-6}$ kgmm⁻³ and Poisson’s ratio $\nu =0.22$. We model the CT specimen with three simulations: one using a fine mesh throughout the entire domain, one coupling a coarse ${\Omega }_L$ and fine ${\Omega }_s$ mesh with a single $\Delta t$, and another with ${\Omega }_L$ and ${\Omega }_s$ integrating with multiple time steps $\Delta t_L$ and $\Delta t_s$, respectively. Structured meshes are used in all cases where an element size of $0.58$ mm for fine and 1 mm for coarse. All simulations use $Co=0.5$, running for a maximum of $t_{final}=0.02$ ms.

In Figure 7 and Figure 8, we plot the stress contours ${\sigma }_{\mathit{yy}}$ at $t=0.01408$ ms and $t=0.01966$ ms, respectively, where the stress concentration can be visualised at the crack tip of the specimen. Figure 7a and Figure 8a capture the reference mesh, whilst Figure 7b and Figure 8b capture the spatially coupled mesh on a single time step $\Delta t_s$. Figure 8 accurately predicts maximum stresses of 0.0396 GPa vs. 0.0410 GPa for reference against non-matching mesh.

Further quantification on the performance of the non-matching mesh algorithm can be made with the evaluation of stress intensity factors, near the crack tip. Utilising the stress extrapolation method, stress values in the vicinity of the crack tip can be used to simultaneously solve for in-plane $K_I$ and shear $K_{II}$ factors [41,42]. In rectangular coordinates, these are given as follows:

$$\begin{array}{c}\hfill {\sigma }_{yy}={\displaystyle \frac{K_I}{\sqrt{2\pi r}}}cos{\displaystyle \frac{\theta }{2}}\left(1-sin{\displaystyle \frac{\theta }{2}}sin{\displaystyle \frac{3\theta }{2}}\right)-{\displaystyle \frac{K_{II}}{\sqrt{2\pi r}}}sin{\displaystyle \frac{\theta }{2}}\left(2+cos{\displaystyle \frac{\theta }{2}}cos{\displaystyle \frac{3\theta }{2}}\right)\end{array}$$

(36)

$$\begin{array}{c}\hfill {\sigma }_{xx}={\displaystyle \frac{K_I}{\sqrt{2\pi r}}}cos{\displaystyle \frac{\theta }{2}}\left(1+sin{\displaystyle \frac{\theta }{2}}sin{\displaystyle \frac{3\theta }{2}}\right)+{\displaystyle \frac{K_{II}}{\sqrt{2\pi r}}}\left(cos{\displaystyle \frac{\theta }{2}}sin{\displaystyle \frac{\theta }{2}}cos{\displaystyle \frac{3\theta }{2}}\right)\end{array}$$

(37)

where r represents the radial distance from the crack, and $\theta$ the angle relative to the crack plane. Whilst J-integral and virtual crack closure techniques (VCCTs) have been developed for more accurately evaluating K, here the consistency in the method for both reference and coupled solution proves a sufficient comparison [43,44]. The stress ${\sigma }_{xx}$ and ${\sigma }_{yy}$ in the nearest element to the crack correspond to Figure 9a where strong alignment between spatially coupled and reference solution can be observed. In Figure 9b, the analysis computes $K_I$ and $K_{II}$ with $r^{\mathrm{ref}.}=0.1294$ mm, ${\theta }^{\mathrm{ref}.}=1.{040}^{\circ }$, $r^{\mathrm{coup}.}=0.1320$ mm and ${\theta }^{\mathrm{coup}.}=1.{015}^{\circ }$ with the slight variation in the spatial discretisation of ${\Omega }_s$ and summarising in

Table 1. Comparison of stress intensity factors $K_I$ and $K_{II}$ for the compact-tension specimen test with reference mesh and spatially coupled mesh at final time $t=0.02$ ms.

Simulation	${\mathit{K}}_{\mathit{I}}$ [MPa·${\mathbf{m}}^{\frac{\mathbf{1}}{\mathbf{2}}}$]	${\mathit{K}}_{\mathit{II}}$ [MPa·${\mathbf{m}}^{\frac{\mathbf{1}}{\mathbf{2}}}$]
Reference (monolithic)	0.02714	−0.002857
Spatial coupling	0.02810	−0.003036

. $K_I$ and $K_{II}$ compute a difference of 3.5% and 6.3% in value when comparing uncoupled and coupled solutions. This is deemed reasonable granted the approximation of the stress extrapolation method. For future comparisons, the contribution of multiple Gauss integration points could be considered with the aforementioned J-integral or VCCT methods.

Now with multi-time stepping, ${\Omega }_L$ and ${\Omega }_s$ step with $\Delta t_L=3.082\times {10}^{-5}$ ms and $\Delta t_s=5.315\times {10}^{-6}$ ms, producing a time-step ratio $m=5.80$. Through combining spatial and temporal coupling, we observe similar (<$1\times {10}^{-6}$ GPa) results in ${\sigma }_{\mathit{yy}}$ distribution, as seen in Figure 10a. The difference in the coupled simulations is captured via Figure 10b, where the root mean square error (RMSE) is plotted. The stress ${\sigma }_{\mathit{yy}}$ is linearly interpolated for the non-matching mesh to allow for comparison of stress on the coarse mesh’s Gauss points. The couplings capture nearly the entire specimen within a <5% error. Considerable error is located around the circumference of the pins, where the chequerboard pattern indicates potential hourglassing in ${\Omega }_L$. To mitigate such issues, higher order or fully integrated elements could be used. However, in the presented results, we note that no hourglass control methods have been applied for fair comparison of the reference and coupled solutions. The other portion of error can be seen in the far right of the specimen in Figure 10; however, like the pins, these Gauss points reside far from the area of interest at the crack tip. To avoid such errors, a smaller element size would be required in ${\Omega }_L$; however, this raises the trade-off between accuracy and computational efficiency. As per the temporal coupling, an energy balance in each of the two subdomains was achieved $W_{bal,i}<\times {10}^{-12}$ kNmm, orders of magnitude below individual components of energy. The computational efficiency is summarised with speedup achieved utilising both couplings, as presented in

Table 2. Computational runtimes and speedup vs. reference (single $\Delta t_s)$, of the dynamically loaded CT specimen. Reference mesh size is 0.58 mm throughout the domain, whereas coupled simulations utilise 0.58 and 1.0 mm for ${\Omega }_s$ and ${\Omega }_L$, respectively. $\Delta t_L$ and $\Delta t_s$ used for temporally coupled run.

Simulation	Runtime [s]	Speedup
Reference (monolithic)	7428	-
Spatially Coupled	2267	3.27×
Spatially and Temporally Coupled	572	12.98×

. Whilst a modest speedup is achieved with non-matching meshes, an even larger efficiency is found with coupling in both space and time. Considering that coupling with non-matching meshes and combined coupling in space and time result in similar accuracy, the addition of multi-time stepping to these simulations seems obvious.

5. Conclusions

We present coupling methods for the dynamic modelling of solids with explicit finite elements, temporally and spatially. When modelling composites, constituents have varying dilatational wave speeds, hence different time steps. Integrating over the smallest time step can prove highly inefficient, hence the need for multi-time stepping. Our method allows partitions of a domain to solve with their own respective time step, regardless of time-step ratio, hence reducing computational overhead. The method avoids solving a system of equations on the interface, unlike many of the methods that employ Lagrange multipliers. The stability of the method is assessed through evaluation of the subdomains’ energy balance. Very-high-frequency content or variations in motion below the time step of elements on the interface could result in spurious oscillations being generated. This potential limitation promotes the development of adaptive multi-time-stepping schemes that maintain the stability of the interface, withstanding high-frequency stress waves.

The addition of the coupling in space solves the issue of non-matching meshes so that small element sizes are only required in regions of interest. Coupling operators are easily implemented, without increasing the degrees of freedom on the interface. The method avoids the storage of large sparse matrices, reducing computational memory requirements. Ongoing work addresses the spatial coupling with quadrilateral non-matching interfaces in 3-D domains [45]. Numerical examples capture an increase in efficiency with stress wave propagation in a heterogeneous bar, and the modelling of a compact-tension specimen. Both couplings in time and space reduce computational runtimes when compared to their monolithic simulation, especially when combined. Limitations that concern the combination of spatial and temporal coupling include the computation of operators $\mathbf{N}$ and $\mathbf{R}$ solely at $t=0$. Whilst this suffices for the benchmarks shown, for larger deformations this assumption is likely to require further development. Geometric representations of a non-matching $\Gamma$ that are non-planar have still yet to be explored. This proves an important topic as these couplings are applied to real-world multi-scale problems.

Future work looks at the coupling between macro- and meso-scale meshes [46,47,48], with adaptivity a clear necessity for these multi-scale couplings [49]. Whilst linear elasticity is a fair assumption for the rates of deformation demonstrated, other constitutive models should be investigated to evaluate the performance of the couplings, with further reductions in time steps and element distortion. In parallel, experimental fields that require efficiently capturing wave propagation with explicit finite elements are widespread [50,51]. Other coupling opportunities are plentiful when considering dynamic applications; the modelling of contact [14,16], composite fractures [38,39], fluid–structure interactions [13,27], and other impact engineering scenarios are just a few worth mentioning.