An Explicit-Implicit Mixed Staggered Asynchronous Step Integration Algorithm in Structural Dynamics

Many engineering applications need to analyse the system dynamics on the macro and micro level, which results in a larger computational effort. An explicit-implicit asynchronous step algorithm is introduced to solve the structural dynamics in multi-scale both the space domain and time domain. The discrete FEA model is partitioned into explicit and implicit parts using the nodal partition method. Multiple boundary node method is adopted to handle the interface coupled problem. In coupled region, the implicit Newmark coupled with an explicit predictor corrector Newmark whose predictive wave propagates into the implicit mesh. During the explicit subcycling process, the variables of boundary nodes are solved directly by dynamics equilibrium equation. The dissipation energy is dynamically determined in accordance with the energy balance checking. A cantilever beam and a building two numerical examples are proposed to verify that the method can greatly reduce the computing time while maintaining a high accuracy.


Introduction
Explicit and implicit time integration algorithms are the two direct integration algorithms in the transient dynamic analysis. With the improvement of numerical simulation method, there is a growing demand to solve dynamic problems with different temporal and spatial scale. The explicit and implicit algorithms are often combined to solve the problem, in order to take advantage of different algorithms [Hughes (2012)]. There are two kinds of combination approach, one is explicit and implicit algorithms used sequentially to solve the problem in the time domain, the other is to use different family algorithm in different spatial domain. In the first case, a reasonable approach is to integrate explicit and implicit algorithms into a unified program. The algorithm shift from one family to another according to the appropriate transition condition. It is necessary to deal with the numerical oscillation of explicit and implicit switching [Jung and Yang (1998)]. This kind of method is used to predict springback in sheet metal forming by Narkeeran and Lovell [Narasimhan and Lovell (1999)]. Noels used this explicit-implicit combined method to handle crashworthiness analysis [Noels, Stainier and Ponthot (2004); Noels (2006)]. Bock presented a unified framework for the numerical solution of optimal control problems solved with both implicit and explicit switches [Bock, Kirches and Meyer (2018)]. For another type of combination, the element partition and nodal partition are applied. The structure is divided into a number of subdomains. Each subdomain has its own integration scheme and time step. This kind of research originates from Belytschko and Hughes [Belytschko and Mullen (1978); Hughes and Liu (1978)]. In the later decades, many representative algorithms have emerged. The algorithms can be divided into two categories: with or without overlapping region. In the context of non-overlapping subdomain coupling method, Schur-type approaches and dual Schur-type approaches are two typical methods. The Schur-type method enforces the displacement or acceleration continuity at the interface [Mahjoubi, Gravouil and Combescure (2009)]. Liu et al. [Liu and Belytschko (1982)] developed a new family of implicit-explicit algorithm which permits different time integration and different time steps to be used simultaneously by element partitioning. Smolinski et al. [Smolinski, Belytschko and Neal (1988); Wu and Smolinski (2000)]. Daniel proposed an implicit subcycling algorithm [Daniel (2003)] developed an explicit subcycling algorithm by nodal partitioning. For dual Schur-type approsches, Gravouil et al. [Gravouil and Combescure (2002)] proposed a multi-time-step dual Schur approach with a velocity continuity condition on the interface at the fine time step,which referenced as GC method. Prakash et al. [Prakash and Hjelmstad (2004)] developed the PH method which enforced continuity of velocity at the coarse time step but only valid for two subdomains. Karimi et al. [Karimi and Nakshatrala (2014)] compared the GC method and PH method in detail and a new multitime-step monolithic coupling method combining the advantages of both methods is developed. For overlapping subdomain coupling method, The Arlequin framework is a representative method [Ghanem, Torkhani, Mahjoubi et al. (2013)]. In recent years, explicit-implicit mixed method is gradually combined with adaptive time step or improve the performance and application scope of the algorithm [Soares (2018) ;Fernandes, Cardoso and Mansur (2017)]. Fekak extended the explicit-implicit method to the field of nonsmooth dynamics [Fekak, Brun, Gravouil et al. (2017)]. Dimarco developed a high order implicit-explicit linear multistep methods for kinetic equations [Dimarco and Pareschi (2017)]. In simulation of three-dimensional grapheme, a hybrid implicit-explicit finite-difference time-domain method is presented [Chen and Wang (2016)]. When mixed with different time steps, the assumption of linear acceleration or linear velocity is applied to tackle the partition boundary conditions. The interpolation process of Lagrange multipliers is used to deal the interface problem .This approach is often regarded as the source of error accumulation and the development of instability, which has led to a bigger error of the calculation result and a harsh condition [Biesiadecki and Skeel (1993)]. In this article, we use the dynamic multiple boundary method to handle the explicit and implicit boundary information transmission. The variable of the boundary nodes is solved directly by explicit transfer process. Numerical cases show that our modification results in a higher accuracy than classical interpolation method. The proposed method can increase computational efficiency by using nodal partition and asynchronous time step method. The rest of this paper is organized as follows: In Section 2, the Newmrak integrator and explicit Newmark based on predictor-corrector formula schemes are introduced to solve the structural dynamic equation. Then, In Section 3 The asynchronous explicit-implicit integration process is proposed to use different time scales and integral methods in different regions. In Section 4, The check of energy balance for the explicit-implicit mixed asynchronous step method is described. In Section 5, two numerical examples are conducted to verify the rationality and feasibility of the explicit-implicit mixed asynchronous method. Finally, conclusions are drawn in Section 6.

The Newmark discretized solution procedure
The finite element discretized of a structure which leads to the semi-discretized dynamic equilibrium equation can be written as Cook [Cook (2007) Where M is the mass matrix, a is the nodal acceleration vector, C is the damping matrix, f int and f ext are the vectors of internal and external force respectively. The given initial displacement and velocity conditions is u(t 0 )=u 0 and v(t 0 )=v 0 . For linear elasticity, the equilibrium equation is a second-order time differential equation and the internal force vector f int can be written as int ( ( )) Where K is the symmetric stiffness matrix, v and u are the nodal velocity and displacement vectors. The Eq. (1) can be solved by discretized in time among which the Newmark based time integrator is mostly used.

The Newmark integrator scheme
The Newmark method has two parameters γ and β . Displacement, velocity and acceleration vectors of generalized structural nodes from time t n to time t n+1 are based on the equilibrium of the dynamic equation [Cook (2007)]. Subscript n is used to denote the time step. At the time step t n+1 ( ) , a n+1 are the displacement ,velocity and acceleration vectors at time step n+1 respectively. u n , v n , a n are the displacement ,velocity and acceleration vectors at time step n. t ∆ is the time step. In the Newmark integral method, the displacement u n+1 in time t n+1 is solved from the dynamic equation The above equation can be rewritten as The matrix K eff in the above equation is referred as the effective stiffness matrix. The matrix K eff and external force vector 1 n+ f are defined as In all time steps the effective stiffness matrix K eff remains constant in linear dynamic analysis.

The explicit newmark based on predictor-corrector scheme
In the predictor-corrector time discretion of explicit Newmark scheme, the predictors of the scheme can be written as Hughes [Hughes (2012)] If the lumped mass matrix is adopted, the mass matrix M is diagonal and the method is explicit. The nodal acceleration vector a n+1 can be solved from the dynamic equation Then, the corrector displacement u n+1 and velocity v n+1 vectors of n+1 time step can be obtained from the equations With the provided initial displacement u 0 and velocity v 0 conditions, the initial acceleration a 0 can be solved as follows 3 The explicit-implicit coupling method with overlapping node

Notation and tag for explicit-implicit coupling
There are two partition methods, namely mesh partition and node partition. The overlapping node partition method is used. The elements fall into three groups: Explicit region, implicit region and interface elements. The schematic of node partition with overlapping node is shown in Fig. 1.
The following notation represents the value of a quantity at subdomain time level. X can indicates displacement, velocity or acceleration n is the number for system time step and j is for subcycling time step. The asynchronous time step in system time level is shown in Fig. 2 Figure 2: Suddomain and system time step

The asynchronous step explicit-implicit integration coupling method
The time steps for implicit and explicit region are ∆t and ∆t r . The step ratio satisfies ∆t=m∆t r . The multi-layer boundary node method is used to tackle the asynchronous step explicit-implicit coupling algorithm. Multi-layer overlapping node consists boundary node and external node. In order to explain the method clearly, we take the step ratio m=3 for example. Data matching for multi overlapping node is shown in Fig. 3.  I  E  EI  E  EI  n  n  n  n  I  I  I  I  I  IE  I  IE  I  n The damping matrix C is assumed in the form of Rayleigh damping. While different from the Belytschko method [Belytschko and Mullen (1978)], the predictor corrector form explicit Newmark scheme is used to solve the nodal acceleration of explicit region. During the explicit subcycling process, the explicit region need to be updated all nodes through time step m-1 The supermark E represents all explicit nodes. The nodal acceleration of explicit region can be solved as Where j=1, 2, 3…m. The displacement and velocity of the explicit non-overlapping nodes are corrected by the Eq. (10). The displacement of implicit region is calculated after explicit subdomain. The expression is as follows The subscript n+1 means system time ( With the increase of explicit subcycling time step, the external node which can be solved correctly in explicit region is decrease. The internal and boundary nodes for explicit region can be solved correctly during the subcycling time step. Before the next system step, external nodes need data transmission as I  The proposed method uses the equilibrium equation to solve the acceleration of the boundary node which completely preserves the precision of boundary node.

Energy balance check for numerical stability
In the compute process, any unstable result can lead to pseudo energy. By checking the energy balance, the calculation of stability can be checked in real time. For the linear without structural damping case, the system internal energy w int , work of the external force w ext and kinetic energy w kin are defined as Krenk [Krenk (2006) e is a very small error tolerance limit whose general order of magnitude is 10 -2 . The energy balance calculation can be implemented in each sub region if the overall dynamic system is very large.

Case of a cantilever beam in bending
A cantilever beam with rectangle section is presented. A bending load is subjected to the free end point A. The load time curve is shown in Fig. 5. The geometry and finite element mesh are shown in Fig. 6. The maximum force is 10000 N.  The thickness of the rectangular cross-section is 0.03 m. The material is isotropic, linear elastic with Young's modulus E=210×10 9 Pa, mass density ρ=7800 kgm -3 and Poisson ratio ν=0.3. The beam is discretized into isoparametric quadrilateral element. The mesh size of the element in the coarse grid region is about 3 times of the size in the fined area. The whole element number of the region is 1080. The step of the explicit region is defined by the CFL condition: ∆t cr =5.78×10 -5 s. In this case, we use the step ∆t r =5×10 -5 s for explicit region. The implicit region adopts a large time step ∆t. Under the condition of time step ratio ∆t=6∆t r , vertical displacement of point A is compared with the classical mixed method [Wu and Smolinski (2000)]. The vertical displacements of point A are shown in Fig. 7.  Figure 7: The vertical displacements of point A for different methods Then we compare the accuracy of the algorithm under three different time step ratios m=3, m=12 and m=24. Fig. 8 compares the vertical displacement of point A. Fig. 9 is the local enlarged of Fig. 8. It can be find that the displacement curve is basically coincident and the accuracy of the algorithm is not significantly decreased with the step ratio increase.
Where u the means the theoretical displacement, u prop means the displacement calculated by different method and i is time step. The displacement error for different methods are counted as shown in Tab. 1. Although with the increase of the step ratio m, the two precision indexes is increasing for both methods. The interpolation process of boundary data is not involved in the proposed method and it has higher computational accurarcy than traditional mixed explicit implicit method. The energy curve Fig. 10 shows that the dissipation energy almost zero for step ratio m=12. Namely, the development of pseudo energy does not appear in the computational process. The method is stability for a larger step ratio condition.

Figure 10:
The energy curve of cantilever beam with step ratio m=12 The computation time for different methods is counted as shown in Tab. 2. The two methods can reduce the calculation time. But the traditional mixed method needs to solve the global stiffness matrix at system time level. This has caused an increase in computing time to a certain extent. The proposed multiple partition boundary method can also reduce the computational time while maintain a higher accuracy.

Case of a building in exploding
This case considers the dynamic response of a building under the explosion shock wave. The shock wave pressure, the overpressure peak value, impulse and duration are used to describe the explosion shock wave. The explosion shock can be simplified as a right triangle pulse load on the front of the structure, retaining its peak and duration characteristics [Ngo, Mendis, Gupta et al. (2007)]. The blasting shock wave conditions are set as showed in Tab. 3. The external force is loaded in the form of equivalent nodal load. The damping ratio of the structure is 5%, which load in the form of Rayleigh damping. The parameters of concrete are taken as follows: Mass density ρ=2500 kgm -3 , the material Young modulus E=30 GPa and Poisson ratioν=0.2. The model is discretized into three dimensional hexahedral linear fully integrated elements. The mesh number of the model is 144384. There are 21657 elements in explicit region. The critical time step of explicit region is ∆t cr =1.43×10 -4 s. The time step ∆t r =1.15×10 -4 s is adopted in this case. Explicit-implicit asynchronous time step algorithm is used to solve the whole finite element simulation model and the explicit-implicit partition is shown in Fig. 11.

Figure 11: A building model under explosion shock wave
For the response of building structure under earthquake or other impacts, the inter layer displacement angle is an important parameter. Fig. 12 compares the distribution of the maximum inter layer displacement angle of the building structure under five loading conditions [Jayasooriya, Thambiratnam, Perera et al. (2011)]. We choose the step ratio m=3, which means the time step ∆t=3∆t r is used in the implicit region. The results of Fig. 12 show that the response of the building structure is influenced significantly by the mass of the explosive and the distance of the explosion point. When the 2 and 4 of the working conditions, the maximum inter story drifts angle of the structure is close to 4%. And the structural deformation under the load condition 3 is much larger than that of the other working conditions. The maximum lateral drift of the building is more than 10%, which will produce great damage to the building. Fig. 13 compares the inter layer displacement angle of three different step ratios m at the same load case 2. With the increase of the step ratio m, the three calculation results have the same change rule.
Fig. 14 is the displacement time curve at the top of the building in the direction of the impact of the shock wave. At the same load case 2, displacement time curve of the explicit method, the proposed method and the classical mixed method (reference method) are compared [Wu and Smolinski (2000)]. The maximum displacement reaches 124 millimeters. It is necessary to further analyze the displacement response of the building under various load conditions.  The computation time required for solving the same physical problem is greatly reduced by using the multi-scale discrete grid and the different step size in the time domain. The more nodes are used in large step size, the greater computational efficiency is obtained. The proposed method has the potential for practical application in large scale computation.

Conclusions
In many actual dynamic engineering applications, the system has different physical and mechanical properties. The accuracy requirement of the different parts of the system determines the scale of the finite element discrete mesh, which limits the time step. In this paper, an explicit-implicit multi-time-step method of finite element is proposed to solve this problem. Compared with other time integration analysis methods, this method has several characteristics: (1) According to the mechanical properties and computational accuracy of a different computing domain, in both space and time domain, the model is simulated with different scales. The more nodes are used in large step size, the greater computational efficiency of the theory can be obtained.
(2) The algorithm uses the changing boundary nodes to deal with the coupling between large and small steps. There is no additional assumption on the boundary node information transfer. The larger time step ratio can be used in theory.
(3) The energy balance method can be used to find the generation of pseudo-energy and unstable development in calculation. It is helpful to get more reasonable calculation results in real time to adjust the calculation step.
To a certain extent, the changing boundary nodes method increases the difficulty of preprocessing. If the algorithm can be reasonably parallel and the calculation domain is optimized to achieve load balance, the method must be able to obtain a good performance in solving the super large scale application. It will be the focus of further research.