ASYNCHRONOUS TIME INTEGRATION WHILE ACHIEVING ZERO INTERFACE ENERGY

. This contribution deals with an asynchronous direct time integration of the finite-element model. The proposed method is applied to the phenomenon of wave propagation through an elastic linear continuum. The numerical model is partitioned into individual subdomains using the domain decomposition method by means of localized Lagrange multipliers. For each subdomain, different time discretizations are used. No restrictions for relation between subdomain’s time steps are imposed. The coupling of the subdomains is forced by an acceleration continuity condition. Additionally, we use the a posteriori technique to also provide the displacement and velocity continuity at the interfaces, and hence we obtain exact continuity of all three kinematic fields. The proposed method is experimentally validated using the modified SHPB (split Hopkinson pressure bar) setup.


Introduction
The problem of asynchronous integration has been addressed by many authors throughout the second half of the last century [1,2].These methods were not robust enough and lacked the required accuracy or were unstable [3].We consider the robust and accurate methods to be those of the present day, such as the following [4][5][6][7][8][9].However, even these new methods have drawbacks.The most common ones include: overly large systems of equations, the necessity of integer ratio of time steps of adjacent subdomains, the necessity of having common time levels where the solution of the whole model must be determined, etc.The method presented below eliminates these shortcomings and, in addition, offers the possibility of guaranteeing the continuity of the displacement, velocity and acceleration field at the interfaces.

Methods
Both strong and weak formulations of the domaindecomposed problem are introduced.

Strong formulation
We recall classical problem of linear elasto-dynamics [10]: ( where σ(x, t) and ε(x, t) are the stress and strain tensor, respectively.The scalar function ρ(x) represents the mass density, vector functions b(x) and u(x, t) have the meaning of the intensity of the volume forces and the displacement field, respectively.The 4 th order tensor of elasticity is given by C(x).Further, u 0 and u0 are the initial displacement an initial velocity fields, respectively.Displacement at the boundary Γ D is given by û and t is known traction at the boundary Γ N .The symbol n represents the outward surface normal unit vector and x ∈ Ω is the position of the material point, and t is the time.Finally, Ω refers to the spatial domain of interest and Υ represents the the time domain of interest.The dot superscript has the meaning of the time derivative.Furthermore we assume that the region Ω is composed of n s sub-domains Ω s , s = 1 . . .n s (see Figure 1).At the interface Γ I between subdomains we enforce the condition of the continuity of the displacement field where u I (x) stands for the displacement of the interface.Note, that further, after the discretization of the weak problem, we use the second time derivative of vol.42/2023 Asynchronous time integration with zero interface energy (2) so the acceleration continuity is actually dictated, see (6).

Weak formulation
According to the Hamilton's principle [11] we search for the stationary solution of δH u, u I , ℓ = δH free (u) + δW I u, u I , ℓ , (3) over time interval [0, T ], where δH free consists of the virtual kinetic energy and potential of the system and of the virtual work done by external loads.The virtual work of the interface δW I can be expressed as where u s and u I s are the displacements of the s-th subdomain and of the i-th interface, respectively (see Figure 1).The ℓ s stands for the Lagrange multiplier field of the s -th subdomain.The whole system consists of n i interfaces and of n s subdomains.The i-th interface is connected to n i s of subdomains.The s-th subdomain is connected with n s i interfaces.

Spatial discretization
Using classical spatial discretization by means of finite elements, one obtains the system of governing equations in the form i.e.Equation of motion, Kinematic interface constraints and Interface equilibrium condition, respectively, where M, K, are the mass and stiffness matrices, f is the load vector, B, L are the matrices of the connectivity [12], λ is the vector of localized Lagrange multipliers and u I it the interface displacement.Resulting from the existence of individual subdomains, matrices and vectors have block-diagonal and block form respectively, where each block represents one subdomain.

Isolation of the interface problem
We will demonstrate that the λ and üI unknowns related to the i-th interface can be expressed explicitly from the knowledge of the interface solution of connected n i s subdomains, i.e. from the knowledge of u on Γ I .Now, assume the problem of two subdomains connected via one interface.Using equations ( 5)-( 7), one can extract the problem of this interface as where ũ is the acceleration predictor defined as Finally, unknowns of the system (8) are expressed as where s = 1 . . . 2 and are the interface mass matrix, interface load, boundary mass matrix and boundary load, respectively.

Temporal discretization
Since the proposed asynchronous time scheme can be combined with arbitrary method of discretization [9], we present only the most simple case of central difference method [10] To obtain the acceleration ün+1 one has to solve the system of equations ( 5)-( 7) in time t n+1 .Stable [8] interface temporal discretization is given by

The asynchronous time integration
Conventionally, whole domain is computed with one globally determined time step.In the case of asynchronous integration, each subdomain has its own time step.

∫ t+Dt
The equality of the i-th interface solution assuming whole subdomains Ωs vs. interface regions ΩR s only

Initialization of interface regions
By interface region Ω R i s we mean the area of the subdomain s affected by the propagating wave from the interface i during at least two critical steps of this area ∆t CR i s .Assume subdomains Ω 1 and Ω 2 connected via interface Γ I .Further assume that we know the solution of given interface at time level t.Then, because of the finite speed of wave, we can compute the solution of the interface Γ I at time level t + ∆t C from reduced system Ω R1 − Γ I − Ω R2 instead of assuming full problem Ω 1 − Γ I − Ω 2 (see Figure 2).By the substitution u = R T u r to ( 5)-( 7) we obtain formally same system of equations, however, for reduced problem on interface regions.In the system, new objects appear which are the interface region matrices and vectors.

Algorithm
We have the problem of two subdomains Ω 1 and Ω 2 with the interface Γ I (see Figure 3).The bluehighlighted line stands for the already known solution (note, that the last known solution of the interface exists at time level t n 1 ).Our challenge is to obtain the solution of Ω 2 at t n+1 2 .Actual interface time step has the value of ∆t I = t n+1 2 − t n 1 .We propose the following algorithm: (2.) Complete the solution of Ω R1 in time to satisfy ( 6) and its time derivative respectively. (3.)  Depending on real further time discretization, the algorithm is repeated or alternates between both subdomains.The computation ends once each of the subdomains has reached the final time t end .

Numerical and experimental validation
Wave propagation is investigated through the linearly graded material (1D formulation).Furthermore, the impact of 3 thin rods is simulated (2D axisymmetric formulation) and the results are compared with experimental outputs of OHPB setup.
Further in the plots, the conventional term stands for nondecomposed model with single globally determined time step.The term asynchronous stands for the proposed method.Computational time step is further always set as the half of the critical time step for each subdomain individually.

Graded material (1D)
The bar with linearly changing Young modulus is loaded by rectangular pulse of 1 Pa (see Figure 4).The bar is divided into 5 equaly sized subdomains.Elements size is set to ∆h = 0.01 m.

Experiment (2D axisymmetry)
The setup consists of 3 aligned bars of various materials (see Table 1 and Figure 6).The first bar made from polymethylmethacrylate (PMMA) has an initial velocity that simulates the impact.The second bar is made of steel and the third bar is made of aluminum alloy EN-AW-7075-T6.
In case of the asynchronous integrator each bar refers to one subdomain, that is, 3 subdomains and 2 interfaces in total are defined.The length of the edge of the 2D axisymmetric square element is set to the value of ∼ 1.6 mm.The computational time step is set to the half of the critical one.
To maximize the credibility of contact behavior, only the longitudinal degrees of freedom (perpendicular to the interface) are permanently coupled.Moreover, only the experimental time window when all bars are in the contact is considered (i.e.≈ 1 ms).

Displacement of the interfaces and
strain ε 11 response The displacement of both interfaces (PMMA-Steel and Steel-Al) is plotted (see Figure 7).One can clearly see that bars are in contact for the time period of 1.2 m s, so the numerical output can be directly compared with experimental results up to this time.The history of axial strain ε 11 on the bar's free surface measured at distance x = 2.746 m, close to the location of the 2nd interface, is plotted (see Figure 8).We can see sufficient agreement between both time evolutions during time period of the interest (0-1.2m s).

Energy balance and dissipation
Finally, the energy balance is investigated.The percentage loss of total energy given by the initial velocity of the PMMA bar (see Figure 9) and interface energy (see Figure 9) is shown.Although the term (4) for interface energy has an exactly zero value, the system is dissipating energy during the very beginning of the simulation.This means, that the energy balance distortion is present in the potential and kinetic energy of the bars.

Conclusions
We found the asynchronous scheme useful especially in cases of existence of multiple material region within the model.The scheme presented is suitable for solving various problems of physics, e.g., fluid structure interaction, where the method used for decomposition [13,14] has already been used successfully.The weakness of the algorithm is the distortion of the energy balance of the subdomains, which is, however, sufficiently small.

Figure 1 .
Figure 1.Arbitrary part of the system Ω divided into subdomanins Ωs connected with i-th interface Γ I .

Figure 4 .
Figure 4. Bar setup with linearly varying Young modulus

Figure 7 .
Figure 7.The axial displacement u1 at bar interfaces

4 xFigure 8 .
Figure 8. Strain ε11 response close to the location of the 2nd interface

Figure 9 .
Figure 9.Total energy loss (compared to initial kinetic energy of the PMMA bar) evolution (top) and interface energy evolution (bottom). 2 Figure 3. Initial state of the algorithm

Table 1 .
The stress at time Initial velocities and parameters of the bars used in the experimental setup