Energy-conserving interface dynamics with asynchronous direct time integration employing arbitrary time steps

https://doi.org/10.1016/j.cma.2023.116110Get rights and content

Highlights

  • FEM model is decomposed to subdomains which are individually integrated in time.

  • Each domain is integrated with individually set time step.

  • No restriction for ratio of time steps of connected domain pair are imposed.

  • Continuity of acceleration, velocity and displacement at the interface is satisfied.

  • Individual domains can be time integrated by different schemes.

Abstract

The derivation and implementation of an asynchronous direct time integration scheme for domain-decomposed finite element models is presented. To maximize clarity in the description of the proposed asynchronous integration, the scheme is restricted to the linear-elastic stress wave propagation case. The proposed method allows the integration of individual subdomains with independent time steps. There is no requirement for an integer time steps ratio of the interacting domains while maintaining zero interface energy. The subdomains are connected by the condition of the continuity of the acceleration field at the interface. In addition, the a posteriori technique is applied to satisfy the continuity of the displacement and velocity fields. Another important contribution of this paper lies in the description of the implementation — we offer the reader a general description of the implementation of the case of any number of subdomains with any number of constraints between them, while the basics of the algorithm are explained on a single domain pair. The functionality of the asynchronous integrator is verified by solving selected problems and comparing with analytical solutions and experimental measurements obtained using a Split Hopkinson pressure bar setup.

Introduction

Transient phenomena described by the transport, wave, heat conduction, or unsteady Navier–Stokes equations are generally challenging to solve numerically, since in addition to the spatial domain, the time domain comes into play.

Restricted to direct time integration, one cannot parallelize a solution in the time domain since each new state directly depends on the previous one. This is a fundamental difference from the solution in the spatial domain, where parallelization is a matter of common practice — as evidenced by domain decomposition methods implemented in commercial software such as LS-Dyna or Abaqus [1].

In today’s era of clusters and cloud solvers with theoretically infinite memory, the spatial decomposition parallelization of the solution process becomes a real game changer for huge complex problems with billions of degrees of freedom. Such a complicated problem may consist of strongly heterogeneous subdomains, and it is necessary to adjust the direct time integration process. Different subdomains may require different direct time integration methods (implicit/explicit) and/or even different computational time steps.

The object of the investigation in this paper is to integrate subdomains directly in time as independently as possible with arbitrary time steps and with minimal communication requirements between them. As a crucial requirement, our aim is to preserve the exact continuity of the displacement, velocity, and acceleration fields at the interfaces.

In the framework of this paper, for simplicity in demonstrating novelty, we focus only on the problem of elastodynamics. For spatial discretization, we use the finite element method with domain decomposition using Localized Lagrange multipliers (LLM) [2], [3], [4].

The specificity of LLM over the MM approach lies in the definition of an entity that copies the interface between subdomains — the frame object or simply the frame (see Fig. 1). The frame has its own discretization with degrees of freedom (displacement, velocity and acceleration). One can also explain the difference by comparing to the mentioned mortar method (MM). In MM (see Fig. 1(b)), one usually defines the slave nodes (red circles) on subdomain 1 and the master surface (red line) on subdomain 2. Then, the kinematic continuity is enforced at the intersection (blue crosses) of the master surface with slave nodes, whose kinematics is projected to the master surface. The nodes of the Lagrange multipliers field are also located at these intersections and the interface equilibrium is automatically satisfied due to the equality of both fields λ1=λ2 and hence the existence of a single Lagrange multipliers field λ=λ1=λ2. In the LLM method (see Fig. 1(a)), we denote both subdomains nodes as slave ones and project them at the additionally defined master surface with its own discretization (i.e. the frame). The kinematic continuity is enforced at this master surface. Lagrange multipliers λ1 and λ2 are located at the slave nodes of subdomains 1 and 2 and the interface equilibrium is achieved by the local summations of their projection to the master surface (frame).

In the elastodynamic formulation, the frame has its own mass matrix (averaged projection of the mass distribution on the boundaries of surrounding subdomains) so that one can compute the acceleration of the frame directly from its equation of motion and calculate a posteriori its velocity and displacement, e.g. by means of central differences. The result is that LLM advantageously isolates each subdomain, leading to very straightforward implementation. Another advantage over widely used MM is the solution in cross-points (points/lines of contact of more than two subdomains), which is generally problematic, but LLM solves this problem by definition (see Fig. 9).

The main issues and challenges come when one intends to integrate each subdomain with a different time step — referred to in the literature as multirate [5], multi-timestep [6], heterogeneous [7], [8] or generally, by ourselves preferred, asynchronous [7], [9], [10] integration.

Recalling hyperbolic problems, we have to mention the authors of very first thoughts about his topic, Belytschko and Mullen [11], [12] in 1976 and 1978 and Smolinsky [13] in 1996, which however (as have been historically shown [14], [15]) are of questionable stability and accuracy and are not usable for nonmatching meshes and lack of robustness. The common feature of the methods mentioned so far is that they do not use domain decomposition. The different regions integrated by different steps partially overlap. For finer time-scale regions, the constant velocity or acceleration of the interface obtained from the computation of the coarse time-scale region is enforced, which is purely an estimate. For that reason, let us leave aside these methods of the past century.

The first robust methods came in the 21st century. These methods involve the technique of domain decomposition via Lagrange multipliers. The kinematic constraints at the interfaces must be chosen to satisfy the continuity of the displacement, velocity, or acceleration field.

The very first robust, easy-to-implement and simple algorithm for common engineering needs was introduced by Combescure and Gravouil [9] in 2002. The versatility of the code is evidenced by its existing implementation in the FE code Europlexus [16] in 2003, where the time steps of each subdomain can be arbitrary and time varying, which is necessary when solving nonlinear problems. The trade-off for robustness is a relatively noticeable dissipation of energy at the interface. Prakash and Hjelmstad [6] in 2004 and also Fekak et al. [7] in 2017 improved the original algorithm to be energy-conserving, but only cases of integer ratio of time steps were presented.

In 2014 Gravouil et al. [17] proposed the BGC macro and micro method (linear variation of the Lagrange multipliers over macro/micro time step) based on the seeking the solution with the canceled interface “pseudoenergy”. While the BGC-macro method is energy conserving and challenging to solve since a large system of equations has to be solved globally, the BGC-micro scheme of direct time integration is performed explicitly, but dissipates the energy at the interface.

Another energy conserving method was presented by Subber and Matouš [10] in 2016 and by Seibold and Rixen [18] in 2019. Similarly to the BGC-macro method [17], these methods suffer from the necessity of forward knowledge of time levels, because of the global formulation within the macro time step and hence also from the necessity of the need for the existence of such common time levels, where the solution of all subdomains is sought at once, which is again a strong drawback especially in nonlinear analysis, where the time step of each subdomain may change. The consequence of these restrictions is also the arise of a large system of equations per each macro time step.

All so far listed methods of the 21st century are based on linear interpolation of the Lagrange multipliers in time, which has been proven as a consistent approach for special cases [17] with several restrictions. The second important feature of listed methods is the prescription of the condition of continuous velocity on the interface, which however do not ensure the continuity of the displacement and the in time increasing gap occurs at the interface — this phenomenon is called drifting.

The different approach is taken by Cho et al. [8] — instead of velocity continuity at the interface, acceleration continuity is enforced. The developed time integration method solves the interface problem by purely computing the saddle point problem in an explicit way without any time interpolation of the variables. The result is the conservation of energy of the system. Moreover, Cho et al. proposed the technique of avoiding drifting, which results in preserving the continuity of all three fields – displacement, velocity, acceleration – at the interface, instead of one field chosen. The key is the use of LLM method and hence the existence of the frame, which is time integrated separately from the subdomains and dictates the velocity and the displacement of the subdomain’s interface. On the other hand, a strong drawback is the existence only of the formulation for integer time step ratio and its challenging implementation.

Overall, methods based on Combescure approach [9] do not satisfy the continuity of the displacement field (drifting), while the algorithm proposed by Cho [8] is presented only for the integer ratio of time steps of the pair of subdomains. Detailed categorization is given by the following list of restrictions and associated table (see Table 1):

  • the necessity of common time levels for all subdomains during time integration (see Fig. 2(a)),

  • presented only an integer ratio of time steps of neighboring subdomains Ωs and Ωs+1 interconnected via interface Γ, i.e., Δts/Δts+1=m, where m is an integer (assuming Δts>Δts+1) and Δts is computational time step of the sth subdomain (see Fig. 2(b)),

  • global formulation within the macro time step Δtmacro and hence the necessity of knowledge about the future time levels of all subdomains during this macro time step, i.e., in the given situation (see Fig. 2(a)), all unknowns of time levels t1t5 are included in single large system of equations,

  • drifting, i.e., violation of the displacement and/or velocity continuity at the interfaces,

  • the algorithm was presented only for a single pair of subdomains.

From the above we can understand that asynchronous integration methods are many and still a hot topic. Although the state of the art presented is sufficiently representative for our purposes, a more detailed and clear description of the history and various methods of asynchronous integration can be found, for example, in [17] or [10].

The motivation for canceling restrictions ① and ② is obvious. The length of the time step compared to its critical value has significant influence on results, especially in non-linear cases. Therefore, any limitations about its length may negatively distort the solution (spurious oscillation of the stress, instability in contact analysis, etc.). A global formulation, that is limitation ③, requires knowledge of several further time steps. This is also very limiting in nonlinear analysis, where the length of the time step changes during the solution. We do not even mention the complexity and size of the system of equations in the general case of a huge number of subdomains and interfaces. The negative consequence of ④ is self-explanatory. Limitation ⑤ is a common reason for non-implementation by a potential user of a given method. In particular cases of a complex system of subdomains, such an implementation may not even be possible. The asynchronous integration process is so complicated that it is very difficult to present a clear and complete algorithm. Nevertheless, the author of the method should attempt to do so if he wants to make it easier to spread it to engineering applications.

For the purpose of motivation to eliminate restrictions and limitations, for example ② and ④, we chose two robust methods (Combescure et al. [16] suffering from ④ and Cho et al. [8] from ②) to show the exemplary difference in the solution precision compared to our proposed scheme: The material-homogeneous bar (see Fig. 3) is discretized with different element sizes Δh1 and Δh2 at two intervals 1 and 2. A discontinuous 4-step load prescribed by set of Heaviside step functions H(tt0) acts at the left end of the free bar.

In the solution of Combescure [16] we can clearly see drifting (see Fig. 4 green dash-dotted line), whereas Cho [8] (red solid) and the proposed method (blue dashed) exactly satisfy the displacement continuity.

However, Cho et al. [8] method is capable of solving only the case of the integer ratio of time steps. This means that here one cannot solve both regions using the critical time step, and one of the steps must be shortened. The summary of time steps used for the benchmark follows (see Table 2).

The result is that the region 2 cannot be time-integrated with its critical time step (unlike the region 1 — see Fig. 5(a)), which causes spurious stress oscillations within this region (see Fig. 5(b) red solid line). The detail of the wavefront of the last step is given below (see Fig. 6). On the other hand, Combescure [16] (green dash-dotted) method leads to the best – almost analytical – solution, i.e., rectangular pulses.

In conclusion, the main novelty of the proposed method is the combination of the following properties:

  • Arbitrary computational time steps of each subdomain (elimination of restrictions ① and ②).

  • An explicit formulation of the system is proposed, i.e., subdomains are stepping from the level n to n+1 directly, a large number of small problems is processed (elimination of restriction ③).

  • The exact continuity of displacement, velocity, and acceleration at the interface is maintained (elimination of restriction ④).

  • The algorithm is introduced for any number of arbitrarily interconnected subdomains (elimination of restriction ⑤).

  • Zero value of the interface energy.

  • Not interfering with the classical implementation [19] of Newmark methods or the method of Central Differences (CDM), straightforward implementation of the proposed asynchronous integrator in existing schemes of commonly used codes.

The rest of the paper is organized as follows. The strong formulation of the linear elastodynamics problem is presented in Section 2 followed by the weak formulation in Section 3. The direct solution of the resulting equations is described in Section 4. The central part of the paper is constituted by Sections 5 Preparation for asynchronous solution, 6 Asynchronous time integrator — strategy of the solution and algorithms, where the implementation of the asynchronous time integration is thoroughly described. In Section 7, several numerical examples are presented to demonstrate the accuracy and performance of the proposed method, and experimental validation is provided in Section 8.

Section snippets

Strong formulation

The problem of elastodynamics in a linear elastic continuum is governed by the balance of linear momentum, including Hooke’s law complemented by initial and boundary conditions [19]: σ+b=ρüinΩ×Υ,σ=:ɛ,ɛ=12(u)T+u,u(x,t=0)=u0inΩ,u̇(x,t=0)=u̇0inΩ,u=uˆonΩD×Υ,nσ=tˆonΩN×Υ,where σ(x,t) is the stress tensor, ɛ(x,t) is the infinitesimal strain tensor, ρ(x) is the mass density, b(x) is the vector of the intensity of the volume forces, u(x,t) is the displacement field, (x) is the elastic 4th

Weak formulation discretized via finite elements and localized Lagrange multipliers

Based on Hamilton’s principle for partitioned constrained elastodynamics, the variational form of the Hamiltonian can be written using the following three-field expression [20]: δHu,uf,=0TδTu̇Uu+Wu+Wfu,uf,dt=0,where T=ΩδuρüdV,U=Ωδɛ:σdV,W=ΩδubdV+ΩNδutˆdS,stands for the kinetic and the potential energy and for the work done by external loads. The last term of (8) represents the frame energy Wfu,uf,=i=1nisSiΦis(usuif)dΓ,where the vector fields us(x,t) and s(x,t) represent the

Direct solution of the system

The core of the problem lies in determining the actual acceleration of the system at known displacements and velocities. For a given formulation, one can express explicitly the acceleration of the frame üf, the acceleration of the subdomains ü and also the Lagrange multipliers λ. For this purpose, we can reformulate the equation of motion and define the acceleration predictor ü̃ as follows Mü+Bλ=fKu,ü+M1Bλ=ü̃=M1(fKu),BTü+BTM1Bλ=BTü̃,Substituting from to last equation of (21) we can

Preparation for asynchronous solution

For each subdomain Ωs, it is necessary to introduce nsi subareas ΩrsiΩs, one per each connected subframe Φi, which represents the “immediate surroundings” of this subframe (see Fig. 10). Let us call this subarea interface region.

We can obtain interface regions by a reduction of degrees of freedom of the subdomain’s complement ΩΩrsi. For this purpose, we define the boolean operator R (transformation matrix) such that Ru=R100Rnsu1uns=ur1urns=ur,where Rs and urs consist of Rs=Rs1Rsnis,urs

Asynchronous time integrator — strategy of the solution and algorithms

In this section, the solution strategy and a suitable integration scheme are proposed.

Numerical benchmarks

Two problems of transient stress waves propagating through a bi-material interface in 1D and 2D are calculated. The numerical results of the central difference method are compared in both cases for:

  • conventional computation, i.e., non-decomposed model integrated in time with a single global time step,

  • asynchronous computation, i.e., model subdomains integrated by different time steps (proposed method).

The results obtained by both numerical approaches are then compared with analytical solutions of

Experimental validation

The experimental results from an Open Hopkinson pressure bar (OHPB), which is a modification of the classical split Hopkinson pressure bar (SHPB) method [24], have been used for validation. Generally, the OHPB setup consists of an incident bar, which is used as a striker pushed by an air gun, the specimen (whose material characteristics are measured), and a transmission bar. For a detailed description of the general setup and the specific technology used in the laboratory, see [25], [26], [27].

Conclusions

The presented asynchronous time integration scheme allows the integration of any number of arbitrarily connected subdomains with the zero interface energy dissipation, while exact continuity of the displacement, velocity, and acceleration fields is maintained. The error of the method is revealed in the negligible distortion of the kinetic and potential energy of the subdomains and only during the propagation of the discontinuity of the rectangular stress wave through the interface. Moreover,

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The work of R. D., R. K. M. M. and V. A. were supported by the grant projects with No. 22-00863K of the Czech Science Foundation (CSF) within institutional support RVO:61388998, R.K. a R.R. realized also the work under the Czech-Estonian mobility project EstAV-21-02. The work of R. D. was also supported by the Grant Agency of the Czech Technical University in Prague , grant No. SGS22/196/OHK2/3T/16.

References (30)

  • B. Seny, J. Lambrechts, V. Legat, J.-F. Remacle, Efficient Parallel Multirate Time Stepping for Accelerating Explicit...
  • PrakashA. et al.

    A FETI-based multi-time-step coupling method for newmark schemes in structural dynamics

    Internat. J. Numer. Methods Engrg.

    (2004)
  • FekakF.-E. et al.

    A new heterogeneous asynchronous explicit–implicit time integrator for nonsmooth dynamics

    Comput. Mech.

    (2017)
  • ChoS.S. et al.

    Explicit multistep time integration for discontinuous elastic stress wave propagation in heterogeneous solids

    Internat. J. Numer. Methods Engrg.

    (2019)
  • SubberW. et al.

    Asynchronous space–time algorithm based on a domain decomposition method for structural dynamics problems on non-matching meshes

    Comput. Mech.

    (2016)
  • Cited by (0)

    View full text