Acta Polytechnica CTU Proceedings https://doi.org/10.14311/APP.2023.42.0006 Acta Polytechnica CTU Proceedings 42:6–11, 2023 © 2023 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague ASYNCHRONOUS TIME INTEGRATION WHILE ACHIEVING ZERO INTERFACE ENERGY Radim Dvořáka,b,∗, Radek Kolmanb, Jan Faltaa, Michaela Neühauserováa a Czech Technical University in Prague, Faculty of Transportation Sciences, Department of Mechanics and Materials, Na Florenci 25, 110 00 Praha 1, Czech Republic b Czech Academy of Sciences, Institute of Thermomechanics, v. v. i., Dolejškova 1402/5, 182 00 Praha 8, Czech Republic ∗ corresponding author: dvorara9@fd.cvut.cz Abstract. 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. Keywords: Finite element method, direct time integration, asynchronous integration, domain decom- position, localized lagrange multipliers. 1. 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–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. 2. Methods Both strong and weak formulations of the domain- decomposed problem are introduced. 2.1. Strong formulation We recall classical problem of linear elasto-dynamics [10]: div σ + b = ρü on Ω × Υ , σ = C : ε , ε = 1 2 [ (grad u)T + grad u ] , u = u0(x, t = 0) on Ω , u̇(x, t = 0) = u̇0 on Ω , u = û on ΓD , n · σ = t̂ on ΓN . (1) where σ(x, t) and ε(x, t) are the stress and strain ten- sor, 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 4th or- der tensor of elasticity is given by C(x). Further, u0 and u̇0 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 com- posed of ns sub-domains Ωs, s = 1 . . . ns (see Figure 1). At the interface ΓI between subdomains we enforce the condition of the continuity of the displacement field us(x) = uI(x) on ΓI ∀s , (2) where uI(x) stands for the displacement of the inter- face. Note, that further, after the discretization of the weak problem, we use the second time derivative of 6 https://doi.org/10.14311/APP.2023.42.0006 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en vol. 42/2023 Asynchronous time integration with zero interface energy Ωs Ωs+1 Ωs-1 GI i-1 GI i uIi us ls ls+1 us+1 Figure 1. Arbitrary part of the system Ω divided into subdomanins Ωs connected with i-th interface ΓI. (2) so the acceleration continuity is actually dictated, see (6). 2.2. Weak formulation According to the Hamilton’s principle [11] we search for the stationary solution of δH ( u, uI, ℓ ) = δHfree (u) + δWI ( u, uI, ℓ ) , (3) over time interval [0, T ], where δHfree 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 δWI can be expressed as δWI ( u, uI, ℓ ) = δ ni∑ i=1 ni s∑ s=1 ∫ ΓI [ ℓs(us − uI i) ] dΓ , (4) where us and uI 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 con- sists of ni interfaces and of ns subdomains. The i-th interface is connected to ni s of subdomains. The s-th subdomain is connected with ns i interfaces. 2.3. Spatial discretization Using classical spatial discretization by means of finite elements, one obtains the system of governing equations in the form Mü + Ku = f − Bλ, (5) BT ü − L üI = 0, (6) LT λ = 0, (7) i.e. Equation of motion, Kinematic interface con- straints and Interface equilibrium condition, respec- tively, where M, K, are the mass and stiffness matri- ces, f is the load vector, B, L are the matrices of the connectivity [12], λ is the vector of localized Lagrange multipliers and uI it the interface displacement. Re- sulting from the existence of individual subdomains, matrices and vectors have block-diagonal and block form respectively, where each block represents one subdomain. 2.4. Isolation of the interface problem We will demonstrate that the λ and üI unknowns related to the i-th interface can be expressed explic- itly from the knowledge of the interface solution of connected ni s subdomains, i.e. from the knowledge of u on ΓI. Now, assume the problem of two subdomains con- nected via one interface. Using equations (5)–(7), one can extract the problem of this interface asBT 1 M−1 1 B1 0 L1 0 BT 2 M−1 2 B2 L2 L1T L2T 0  L1λ1 L2λ2 üI  = BT 1 ˜̈u1 BT 2 ˜̈u2 0  , (8) where ˜̈u is the acceleration predictor defined as ˜̈u = M−1(f − Ku) = ü + M−1Bλ . (9) Finally, unknowns of the system (8) are expressed as üI = MI−1 f I , λs = Mλs −1 fλs , (10) where s = 1 . . . 2 and MI = 2∑ s=1 [ LsT ( BT s M−1 s Bs )−1 Ls ] , f I = 2∑ s=1 [ LsT ( Bi s TM−1 s Bi s )−1 Bi s T ˜̈us ] , Mλs = ( BT s M−1 s Bs )−1 , fλs = BT s ˜̈us − LsüI . (11) are the interface mass matrix, interface load, boundary mass matrix and boundary load, respectively. 2.5. Temporal discretization Since the proposed asynchronous time scheme can be combined with arbitrary method of discretiza- tion [9], we present only the most simple case of central difference method [10] un+1 = un + ∆t u̇n + ∆t2 2 ün , u̇n+1 = u̇n + ∆t 2 ( ün + ün+1) . (12) To obtain the acceleration ün+1 one has to solve the system of equations (5)–(7) in time tn+1. Stable [8] interface temporal discretization is given by (u̇I)n+x1 = (u̇I)n + ∆t 2 ( (üI)n + (üI)n+1) , (uI)n+1 = (uI)n + ∆t (u̇I)n+1 . (13) 2.6. 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. 7 R. Dvořák, R. Kolman, J. Falta, M. Neühauserová Acta Polytechnica CTU Proceedings ∫ t+DtC dt = t W1 W2 ∫ dt t WR WR1 2t+DtC on GI GI GI Figure 2. The equality of the i-th interface solution assuming whole subdomains Ωs vs. interface regions ΩRs only W1 W2WR WR1 2GI t1 tI t2 t n1 t n2 t n+1 2 t n+1 1 ΔtI Δt2 Figure 3. Initial state of the algorithm 2.6.1. Initialization of interface regions By interface region ΩRi 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 ∆tCRi 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 + ∆tC from reduced system ΩR1 − ΓI − ΩR2 instead of assuming full problem Ω1 − ΓI − Ω2 (see Figure 2). By the substitution u = RT ur to (5)–(7) we ob- tain formally same system of equations, however, for reduced problem on interface regions. In the system, new objects appear Kr = RKRT , Mr = RMRT , f r = Rf , Br = RB, (14) which are the interface region matrices and vectors. 2.6.2. Algorithm We have the problem of two subdomains Ω1 and Ω2 with the interface ΓI (see Figure 3). The blue- highlighted line stands for the already known solution (note, that the last known solution of the interface exists at time level tn 1 ). Our challenge is to obtain the solution of Ω2 at tn+1 2 . Actual interface time step has the value of ∆tI = tn+1 2 − tn 1 . We propose the following algorithm: ρ 2 = 1 kg m-3, E 2 = 2.25 Paρ 1 = 1 kg m-3, E 1 = 1 Pa 0 L = 1 m E1 = 1 Pa E2 = 2.25 Pa 1 2 3 4 5 Δh = 0.01 m (t) = 1 Paσ Figure 4. Bar setup with linearly varying Young modulus (1.) Solve the interface in time tn+1 2 (a) (ur)n+1 = (ur)n + ∆tI (u̇r)n + ∆tI 2 2 (ür)n, (b) (˜̈ur)n+1 = Mr−1 (( fr)n+1 − Kr (ur )n+1 ) , (c) (üI)n+1 = MI−1 (f I)n+1, (d) (u̇I)n+1 = (u̇I)n + ∆tI 2 (( üI ) n + ( üI ) n+1 ) , (e) (uI)n+1 = (uI)n + ∆tI (u̇I)n+1. (2.) Complete the solution of ΩR1 in time tn+1 2 (a) (λ1)n+1 = Mλ1 1 −1 ( BrT 1 (˜̈ur 1)n+1 − L1(üI)n+1 ) , (b) (ür 1)n+1 = (˜̈ur 1)n+1 − Mr 1 −1Br 1(λ1)n+1, (c) (u̇)n+1 = (u̇)n + ∆tI 2 ( (ü)n + (ü)n+1 ) , (d) Correct {ur, u̇r}n+1 1 to satisfy (6) and its time deriva- tive respectively. (3.) Solve free subdomain Ω2 in time tn+1 2 (a) un+1 2 = un 2 + ∆t2 u̇n 2 + ∆(t2)2 2 ün 2 , (b) ˜̈un+1 2 = M−1(fn+1 2 − K2 un+1 2 ), (c) rewrite ˜̈un+1 2 → ün+1 2 with respect to BT 2 ün+1 2 − L2 (üI)n+1 = 0, (d) u̇n+1 2 = u̇n 2 + ∆t2 2 ( ün 2 + ün+1 2 ) , (e) Correct {un+1 2 , u̇n+1 2 } to satisfy (6) and its time deriva- tive respectively. (f) Reset {ur, u̇r, ür}n+1 2 = R2{u, u̇, ü}n+1 2 . Depending on real further time discretization, the algorithm is repeated or alternates between both sub- domains. The computation ends once each of the subdomains has reached the final time tend. 3. Numerical and experimental validation Wave propagation is investigated through the lin- early graded material (1D formulation). Furthermore, the impact of 3 thin rods is simulated (2D axisymmet- ric 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 deter- mined 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. 3.1. 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. The stress at time 8 vol. 42/2023 Asynchronous time integration with zero interface energy Bar Impact velocity Length Diameter Material Density Young mod. Poisson no. [ms−1] [mm] [mm] ρ [kgm−3] E [GPa] ratio ν [-] 1. 15.2 1750 20 PMMA 1180 5 0.37 2. 0 796 20 Steel 7850 210 0.3 3. 0 1600 20 Al 2806 72 0.334 Table 1. Initial velocities and parameters of the bars used in the experimental setup 0 0.2 0.4 0.6 0.8 1 x [m] 0 0.5 1 1.5 σ [P a] conventional asynchronous analytical Figure 5. Propagated stress pulse at time t = 0.75 s pressure max. 1 MPa barrel solenoid valve 1st bar PMMA laser through-beam sensor damper linear systemguidance (rail and carriage) air flow speckled patterns for DIC high speed cameras strain-gauge 2nd bar steel 3rd bar Aluminium Figure 6. The scheme of used Open Hopkinson pressure bar setup 0 0.2 0.4 0.6 0.8 1 1.2 t [s] ×10-3 0 0.5 1 1.5 2 2.5 3 u [m ] ×10-3 x = 1.750 m (PMMA-Steel - i.e. 1st interface) exp. PMMA exp. Steel asynchronous (a). PMMA-Steel interface 0 0.2 0.4 0.6 0.8 1 1.2 t [s] ×10-3 0 0.5 1 1.5 2 2.5 3 u [m ] ×10-3 x = 2.546 m (Steel-Al - i.e. 2nd interface) exp. Steel exp. Al asynchronous (b). Steel-Al interface Figure 7. The axial displacement u1 at bar interfaces t = 0.75 s is plotted in Figure 5. Spurious oscillations around the analytical solution are caused by the use of a noncritical computational time step ∆ts < 0.5 ∆tC s . 3.2. Experiment (2D axisymmetry) The setup consists of 3 aligned bars of various mate- rials (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 (perpendicu- lar to the interface) are permanently coupled. More- over, only the experimental time window when all bars are in the contact is considered (i.e. ≈ 1 ms). 3.2.1. 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 com- pared with experimental results up to this time. 9 R. Dvořák, R. Kolman, J. Falta, M. Neühauserová Acta Polytechnica CTU Proceedings 0 0.2 0.4 0.6 0.8 1 1.2 t [s] ×10-3 -5 -4 -3 -2 -1 0 1 2 ε [- ] ×10-4 x = 2.746 m (Al bar: 0.2 m from Steel-Al interface) experiment asynchronous Figure 8. Strain ε11 response close to the location of the 2nd interface 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 Fig- ure 8). We can see sufficient agreement between both time evolutions during time period of the interest (0–1.2 m s). 3.2.2. 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. 4. 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 solv- ing various problems of physics, e.g., fluid structure interaction, where the method used for decomposi- tion [13, 14] has already been used successfully. The weakness of the algorithm is the distortion of the en- ergy balance of the subdomains, which is, however, sufficiently small. Acknowledgements The work of Radim Dvořák and Radek Kolman was supported by grant projects with No. 22-00863K of the Czech Science Foundation (CSF) within institutional sup- port RVO:61388998. The work of Radim Dvořák was also supported by the Grant Agency of the Czech Technical University in Prague, grant No. SGS22/196/OHK2/3T/16. Radek Kolman and Radim Dvořák also realized the work under the Czech-Estonian mobility project EstAV-21-02. References [1] T. Belytschko, R. Mullen. Mesh partitions of explicit-implicit time integration. Formulations and computational algorithms in finite element analysis : 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 t [s] ×10-3 -1 0 1 2 3 4 % of E n er gy lo ss [- ] x = global asynchronous 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 t [s] ×10-3 -1 -0.5 0 0.5 1 E n er gy [J ] ×10-32 x = global asynchronous Figure 9. Total energy loss (compared to initial kinetic energy of the PMMA bar) evolution (top) and interface energy evolution (bottom). U.S.-Germany symposium. Massachusetts Institute of Technology, Cambridge, 1977. [2] P. Smolinski. Subcycling integration with non-integer time steps for structural dynamics problems. Computers & Structures 59(2):273–281, 1996. https://doi.org/10.1016/0045-7949(95)00256-1 [3] P. Smolinski, S. Sleith, T. Belytschko. Stability of an explicit multi-time step integration algorithm for linear structural dynamics equations. Computational Mechanics 18:236–244, 1996. https://doi.org/10.1007/BF00369941 [4] A. Combescure, A. Gravouil. A numerical scheme to couple subdomains with different time-steps for predominantly linear transient analysis. Computer Methods in Applied Mechanics and Engineering 191:1129–1157, 2002. https://doi.org/10.1016/S0045-7825(01)00190-6 [5] W. Subber, K. Matouš. Asynchronous space–time algorithm based on a domain decomposition method for structural dynamics problems on non-matching meshes. Computational Mechanics 57(2):211–235, 2015. https://doi.org/10.1007/s00466-015-1228-0 [6] F.-E. Fekak, M. Brun, A. Gravouil, B. Depale. A new heterogeneous asynchronous explicit–implicit time integrator for nonsmooth dynamics. Computational Mechanics 60(1):1–21, 2017. https://doi.org/10.1007/s00466-017-1397-0 [7] A. Prakash, K. D. Hjelmstad. A FETI-based multi-time-step coupling method for Newmark schemes in structural dynamics. International Journal for Numerical Methods in Engineering 61(13):2183–2204, 2004. https://doi.org/10.1002/nme.1136 [8] S. S. Cho, R. Kolman, J. A. González, K. C. Park. Explicit multistep time integration for discontinuous elastic stress wave propagation in heterogeneous solids. International Journal for Numerical Methods in Engineering 118(5):276–302, 2019. https://onlineli brary.wiley.com/doi/pdf/10.1002/nme.6027 https://doi.org/10.1002/nme.6027 10 https://doi.org/10.1016/0045-7949(95)00256-1 https://doi.org/10.1007/BF00369941 https://doi.org/10.1016/S0045-7825(01)00190-6 https://doi.org/10.1007/s00466-015-1228-0 https://doi.org/10.1007/s00466-017-1397-0 https://doi.org/10.1002/nme.1136 https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6027 https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6027 https://doi.org/10.1002/nme.6027 vol. 42/2023 Asynchronous time integration with zero interface energy [9] S. S. Cho, K. C. Park, H. Huh. A method for multidimensional wave propagation analysis via component-wise partition of longitudinal and shear waves. International Journal for Numerical Methods in Engineering 95(3):212–237, 2013. https://doi.org/10.1002/nme.4495 [10] T. J. R. Hughes. The finite element method: linear static and dynamic finite element analysis. Dover Publications, Mineola, NY, 2000. [11] J. A. González, J. Kopačka, R. Kolman, K.-C. Park. Partitioned formulation of contact-impact problems with stabilized contact constraints and reciprocal mass matrices. International Journal for Numerical Methods in Engineering 122(17):4609–4636, 2021. https://doi.org/10.1002/nme.6739 [12] K. C. Park, C. A. Felippa. A variational principle for the formulation of partitioned structural systems. International Journal for Numerical Methods in Engineering 47(1-3):395–418, 2000. https://doi.org/10.1002/(sici)1097-0207(200001 10/30)47:1/3<395::aid-nme777>3.0.co;2-9 [13] K. C. Park, C. Felippa, U. Gumaste. A localized version of the method of Lagrange multipliers and its applications. Computational Mechanics 24:476–490, 2000. https://doi.org/10.1007/s004660050007 [14] K. C. Park, C. A. Felippa, G. Rebel. A simple algorithm for localized construction of non-matching structural interfaces. International Journal for Numerical Methods in Engineering 53(9):2117–2142, 2002. https://doi.org/10.1002/nme.374 11 https://doi.org/10.1002/nme.4495 https://doi.org/10.1002/nme.6739 https://doi.org/10.1002/(sici)1097-0207(200 00110/30)47:1/3<395::aid-nme777>3.0.co;2-9 https://doi.org/10.1002/(sici)1097-0207(200 00110/30)47:1/3<395::aid-nme777>3.0.co;2-9 https://doi.org/10.1007/s004660050007 https://doi.org/10.1002/nme.374 Acta Polytechnica CTU Proceedings 42:6–11, 2023 1 Introduction 2 Methods 2.1 Strong formulation 2.2 Weak formulation 2.3 Spatial discretization 2.4 Isolation of the interface problem 2.5 Temporal discretization 2.6 The asynchronous time integration 2.6.1 Initialization of interface regions 2.6.2 Algorithm 3 Numerical and experimental validation 3.1 Graded material (1D) 3.2 Experiment (2D axisymmetry) 3.2.1 Displacement of the interfaces and strain 11 response 3.2.2 Energy balance and dissipation 4 Conclusions Acknowledgements References