Subdomain Precise Integration Method for Periodic Structures

A subdomain precise integrationmethod is developed for the dynamical responses of periodic structures comprisingmany identical structural cells. The proposed method is based on the precise integration method, the subdomain scheme, and the repeatability of the periodic structures. In the proposed method, each structural cell is seen as a super element that is solved using the precise integration method, considering the repeatability of the structural cells. The computational efforts and the memory size of the proposed method are reduced, while high computational accuracy is achieved. Therefore, the proposed method is particularly suitable to solve the dynamical responses of periodic structures. Two numerical examples are presented to demonstrate the accuracy and efficiency of the proposed method through comparison with the Newmark and Runge-Kutta methods.


Introduction
A periodic structure consists of identical structural cells that are connected together end-to-end to form the entire structure.Precisely because of the repeatability of the periodic structures, these structures exhibit numerous interesting and useful physical properties and are widely applied in many types of engineering, such as railway engineering [1], the pantograph-catenary system [2], and photonic [3,4] and phononic [5] crystals, among others.Currently, because of the importance of periodic structures, many correlational studies have been reported.In [1], a method based on the symplectic mathematical scheme and Schur decomposition was proposed for the random responses of a vehicle moving on an infinitely long periodic track.In [4], Dobson applied finite element discretization coupled with a preconditioned subspace iteration algorithm to periodic dielectric photonic crystals.In 1993, Kushwaha et al. [5] presented the first full band-structure calculations for periodic, elastic composites.Their work introduced a new field of research on periodic phononic crystals.Zhong and Williams investigated the wave propagation problems of repetitive structures [6,7] and the localization phenomenon for the high-frequency vibration modes of imperfectly repetitive structures [8] using an analogy between computational structural mechanics theory and optimal control theory.Wang et al. introduced a lumpedmass method to study the propagation of elastic waves in one- [9] and two-dimensional [10] periodic systems.Wang et al. [11] investigated the free and forced vibration of certain periodic structures using the properties of the structural modes of periodic structures.Mead [12] developed a general theory for the forced vibration of multicoupled, one-dimensional periodic structures.The theory starts from the dynamic stiffness matrix of a single multicoupled periodic element and derives the matrix equations for the magnitudes of the characteristic free waves excited by harmonic forces and/or displacements acting at a single periodic junction.Ding et al. [13] analyzed the wave propagation in a periodic elasticpiezoelectric axial-bending coupled beam by employing the Lyapunov exponent method.In [14], an efficient algorithm was developed for computing the dynamic responses of one-dimensional periodic structures and periodic structures with defects.In [15], the exact solutions for the dynamic response of a periodic spring and mass structure were given.The solutions cover arbitrary initial conditions and both polynomial and harmonic external forces.While many reports concerning periodic structures have been published, developing a method to accurately and efficiently calculate the dynamical responses of periodic structures remains a noteworthy issue.
A primary goal of this paper is to investigate the time integration scheme for the linear dynamical equations of periodic structures produced using the finite element method.To solve the dynamical equation, common methods such as the Runge-Kutta (R-K) method [16][17][18] and the Newmark method [19] can be used.In 1994, Zhong and Williams [20] developed the precise integration method (PIM) based on the accurate computation of the exponential of a matrix, which is widely employed in many types of problems because of its high accuracy.The high accuracy of the method is derived from precisely computing the exponential of the matrix, which, however, requires significant computational efforts.To improve the computational efficiency of PIM, many works have been reported.In [21], an adaptive algorithm of precise integration was proposed by automatically generating the involved parameters in the original PIM in terms of the required accuracy.In [22], Gao et al. proposed the modified fast precision integration method (FPIM), which utilizes the sparse nature of the system matrices and the physical features of the structural dynamics problems to improve the computational efficiency of the exponential of the matrix without loss of precision.Later, the FPIM was extended to the hyperbolic heat conduction problems in [23].In this paper, a subdomain precise integration method (SPIM) is presented for the dynamical responses of periodic structures with many identical structural cells.SPIM is based on PIM, the subdomain scheme [24][25][26], and the repeatability of periodic structures.This method reduces the computational efforts and reduces the amount of memory required during the computation.

Precise Integration Method for Dynamical Systems
The subdomain precise integration method is based on the original PIM; thus, a brief introduction for the original PIM is given in this section.Suppose that the dynamical equation for a periodic structure (see Figure 1) produced using the finite element method (FEM) can be written as with the initial conditions where M, C, and K are the  ×  mass and damping and stiffness matrices of the periodic structure, respectively, and  is the number of degrees of freedom (DOFs) of the periodic structure.X, Ẋ, and Ẍ are the ×1 displacement and velocity and acceleration vectors, respectively, and the dot over the variable denotes differentiation with respect to time .f ∈  ×1 is the load vector, and  0 denotes the initial time.
Many methods can be employed to solve (1), such as the Newmark method, the R-K method, and PIM.If PIM is used to solve (1), it should first be rewritten in the state space: where In ( 4), V = Ẋ denotes the velocity vector, I ∈  × is the unit matrix, and 0 ∈  × is a zero matrix.For numerical integration, the time domain is divided into a series of time intervals of equal length ; that is,  0 = 0,  1 = , . . .,   = ,  +1 = ( + 1) , . . . .
By letting the solution to (3) at discrete times is given by in which The matrix T defined by ( 8) is called the exponential of the matrix H.To evaluate the second term on the right-hand side of (7), the load vector f(  + ) can first be approximated using the Lagrange interpolation polynomial [27] in the time interval [  ,  +1 ].Hence, if we select  + 1 interpolation points, denoted by   +   , in the time interval [  ,  +1 ], the load vector can be approximated as for which where   is the th Lagrange interpolation polynomial,  = 0, 1, . . ., .In terms of ( 7) and (9), the state vector at  +1 is obtained: for which Shock and Vibration 3 and Ψ  can be computed using numerical integration methods, such as the Gaussian quadrature method [28], which has been studied and proved to be an efficient method for determining the integral term in PIM.In terms of ( 8) and ( 12), it can be easily observed that computing the exponentials of matrices efficiently and accurately is the key issue when computing the responses of dynamical systems, which is most likely to be overcome using PIM.Based on the addition theorem, PIM uses the idea of only computing the incremental part of T to improve the computational precision.
The implementation of PIM for computing the exponential of the matrix H and time step  can be given as follows.
Let H 1 = H/2  , in which  is an integer; (8) then becomes If  is large enough, such that the infinite norm of H 1 satisfies ‖H 1 ‖ ∞ ≪ 1, then exp(H 1 ) can be approximated by the Taylor series of order ; that is, Theoretically, if both  and  are large enough, ( 14) yields a very accurate approximation of exp(H 1 ).However, the round-off error may be significant in numerical computation because of the sum of the unit matrix I and the small rest terms.Therefore, in PIM, the exponential matrix exp(H 1 ) is divided into two parts; that is, Then, using ( 13), (14), and (15), T can be written as where After computing R  , T can be given by PIM is an accurate algorithm for computing the matrix exponential.However, the computational effort of PIM is ( 3 ), which is very large for the FEM model of periodic structures with large DOFs.For a periodic structure comprising many identical structural cells, the mass, damping, and stiffness matrices of all the structural cells are the same.With an increase in the number of structural cells, the mass, damping, and stiffness matrices of the whole structure will be too large, which leads to large computation and memory requirements during the computation of the matrix exponential.This issue restricts the application of PIM.In the next section, the repeatability of the periodic structure is combined with the subdomain technique to improve the efficiency of PIM.

The Subdomain PIM
Suppose that a periodic structure comprises  identical structural cells, and the mass, damping, and stiffness of the th identical structural cell are M  , C  , and K  , respectively.The number of the DOFs of the th structural cell is   , and the displacement and load vectors are denoted by X  and f  , respectively.The dynamical model of the th structural cell is given by The DOFs of each structural cell are contained within the DOFs of the whole periodic structure.Hence, the displacement vector X  of the th structural cell can be found in the displacement vector X of the whole structure; that is, where D T  ∈    × is a matrix whose elements are 0 or 1.

D T
describes the relationship between X  and X and satisfies in which I  ∈    ×  is a unit matrix.Multiplying both sides of ( 19) by D  and substituting ( 20) into (19), the dynamical equation of the th structural cell can be rewritten as In terms of (22), the entire dynamical equation for the whole periodic structure can be written as where Actually, (23) is the same as (1).For a periodic structure comprising a large number of identical structural cells, the entire dynamical equation is too large to be solved using PIM directly due to the computational efforts and required storage.However, the scale of the dynamical equation of each structural cell is small compared to the scale of the entire dynamical equation.In addition, the dynamical equations of structural cells are the same, except for the structural cells near the boundary of the periodic structure.Therefore, using PIM for the dynamical equation ( 19) of the structural cell could be a better approach for two reasons.The first reason is that the size of the cell dynamical equation is much smaller than that of the entire dynamical equation.Hence, the cost for the computations of the exponential matrices corresponding to the structural cell is insignificant compared with the cost for that corresponding to the whole periodic structure.The second reason is that the repeatability of the periodic structure can be utilized.Because the structural cells are the same, the corresponding exponential matrices are also the same.Once the exponential matrix for one of the structural cells is computed, the rest are also obtained.Hence, the required storage would be reduced substantially.Meanwhile, a large time step could be selected to ensure accuracy because the exponential matrix is computed using PIM.The exponential matrices of the structural cell involved could be computed using FPIM [22,23], which was modified by Gao et al. based on the original PIM to further improve the computational efficiency.
Before using PIM for the cell dynamical equation ( 19), some issues should first be discussed.The first issue concerns the force vector f  .In (19), the force vector f  actually contains two parts, denoted by p  and q  .p  is the force produced by outside factors, such as gravity, and q  is the force produced by the neighboring structural cells.For the whole periodic structure, q  is the so-called internal force, which means The force vectors f(), p  (), and q  () can also be approximated as by selecting 2 Lagrange interpolating points [27] in the time interval [  ,  +1 ].In (27), p , , q , , and f , are the   × 1 force vectors at the th interpolating time point, and q , is an undetermined vector.Substituting ( 27) into (19) yields which can also be rewritten as ) , Note that p  is known and q  is undetermined.Selecting 2 interpolating points to approximate the force vector q  implies that there are 2 unknown vectors (i.e., q , ,  = 1∼2) that should be determined.To determine these vectors, take  derivatives of (29) with respect to time: where Then, applying PIM to (31) yields where  ,, ∈  2  ×  is given by which can be computed using the Gaussian quadrature method [28].In terms of (33), we have ,+1 U (1)   ,+1 . . .
) , where Y  ∈  2  ×2  is a square matrix and p , ∈    ×1 is the force vector produced by the neighboring structural cells at the th interpolating point.Therefore, in the vector p , , the elements corresponding to the DOFs inside the structural cell are all zeros.For the th structural cell, the displacements and velocities of the DOFs inside the cell are denoted by X , and V , , respectively.The displacements and velocities of the DOFs that are on the boundary of the cell are denoted by X , and V , , respectively.The subscripts  and  indicate internal and boundary DOFs, respectively.Meanwhile, let the state vector U ,+1 be partitioned into and rewrite (36) in the block matrix form in which W ,,+1 corresponds to the DOFs inside the structural cell, and W ,,+1 corresponds to the DOFs on the boundary; that is, ,,+1 U (1)   ,,+1 . . .
) ) , ) . (40) In terms of (39), we have Equation ( 41) can be seen as the dynamical stiffness equation of the th structural cell, and Q , can be seen as the boundary force produced by the neighboring cells.In terms of (41), only the response information about the boundary DOFs is involved; that is, the dynamical stiffness equation (41) of the th structural cell is reduced, and the th structural cell is treated as a super element.Let the displacements and velocities of the DOFs on the boundaries of all structural cells be denoted by X  and V  , respectively, and we have where D T , describes the relationship between X , and X  and D T , D , is a unit matrix.Combining (38) and ( 40) with (43) yields where ) , ) , U (1)   ,+1 . . .
In terms of (46), we have Q , denotes the boundary force produced by the neighboring cells; that is, it is the so-called internal force for the entire periodic structure; hence, we have Combining ( 47) with ( 48) yields In terms of (49), W ,+1 can be computed.Once W ,+1 is obtained, the vector W ,,+1 can be computed using (44).Then, W ,,+1 can also be computed by applying W ,,+1 to (42).When W ,,+1 and W ,,+1 are determined, the displacements and velocities of all the DOFs in each structural cell are obtained, which means that the computation over the time interval [  ,  +1 ] is complete.In the same manner, the displacements and velocities of all the DOFs at each time node can be computed once the initial conditions are given.
For a periodic structure consisting of many identical structural cells, the mass, damping, and stiffness matrices of the structural cell are the same except for the structural cells constrained by the boundary.Hence, for the identical structural cells, the involved exponent matrices only need to be computed once.It can be observed that the proposed SPIM is particularly suitable for the dynamical problems of periodic structures.The number 2 of the interpolation points for approximating q  () could affect the performance of SPIM.If the interpolation points in SPIM are too numerous, a high accuracy can be obtained; however, the computational efficiency will be unsatisfactory.In the next section, the effect of  on the performance of SPIM will be discussed using two numerical examples.

Numerical Examples
Example 1.Consider the periodic structure consisting of masses and springs depicted in Figure 2. The entire periodic structure is composed of 50 structural cells, shown in Figure 2(a), with its masses and stiffnesses being  1 = 1 (kg) and  2 = 2 (kg) and  = 1 (N/m), respectively.For each structural cell, there are 50 DOFs.The entire periodic structure contains 2500 DOFs (see Figure 2(b)) and is fixed at both ends.The initial conditions are that the initial displacement of the 1226th DOF is unity, the initial displacements of the rest of the DOFs are zeros, and the initial velocity is zero for all DOFs.The external force is not considered in this example.
To study the effect of the number 2 of interpolation points on the accuracy of SPIM,  = 1 and 2 are used.To investigate the performance of SPIM, the Newmark method (i.e., the average acceleration implicit scheme) and the variable step R-K method (ODE45 in MATLAB) are employed.The interval of the integration is set to be [0, 1000] s.The time step for SPIM is 0.1 s, two time steps of 10 −3 s and 10 −4 s are used for the Newmark method, and 10 −13 is used for both the absolute and relative tolerances in the R-K method (ODE45).The numerical results computed using ODE45 are considered the reference solution to test the accuracy of the SPIM results.The relative errors of the SPIM and Newmark methods are defined as where X RK and V RK are the displacement and velocity vectors computed using the R-K method, respectively.X and V are the displacement and velocity vectors, respectively, computed using the SPIM or Newmark methods.The subscript 2 represents the 2-norm of a vector.
In this example, the order of the proposed method is also tested by using different time steps.If the order of an algorithm is , the numerical error for time step  is   , in which  is a constant.Thus, the errors  1 and  2 for two different time steps  1 and  2 can be given by and the order of the algorithm can be shown by In this example, the errors of the displacement vectors are defined by (50); thus, the order of the proposed method is defined as  Four time steps  = 0.1 s, 0.2 s, 0.3 s, and 0.4 s are selected to determine the order of the proposed method.Figures 3(a) and 3(b) give the displacements and velocities of all the DOFs computed using the proposed method at 1000 s, respectively.The relative errors of the displacements and velocities computed using the SPIM and Newmark methods are shown in Figures 4(a) and 4(b).Figure 5 shows the variations of the errors of the displacement vectors versus the time step.The order of accuracy and CPU time for all methods are given in Table 1.
Figures 4(a) and 4(b) demonstrate that the results obtained using the SPIM with  = 2 are more accurate than those obtained using the SPIM with  = 1 when the time steps used for the two cases are both 0.1 s.From Figure 5 and Table 1, it can be concluded that the order of accuracy of the SPIM is approximately 2( + 1).However, in terms of Table 1, the SPIM with  = 2 requires slightly more CPU time than the SPIM with  = 1.Compared with the Newmark method, the SPIM results are more accurate, although the time step used for SPIM is 100 times that used for the Newmark method.It can be observed from Table 1 that the CPU time of the R-K method is 8.9 and 10.7 times that of the SPIM with  = 1 and  = 2, respectively.In addition, compared with the Newmark method, SPIM remains more efficient because the CPU times of the Newmark method with two different time steps are 7.5 and 74 times that of the SPIM with  = 2, respectively.It can be concluded from Figure 4 and Table 1 that SPIM using a large time step can achieve high accuracy and efficiency.The models of the structural cell and the entire structure are shown in Figures 6 and 7, respectively.The entire phononic crystal is composed of 15 × 15 structural cells.The bottom of the phononic crystal is fixed.On top of the phononic crystal acts the uniform load () = 2 × 10 7 N/m, which is independent of time.After using the four-node linear rectangular elements for the phononic crystal, the number of the total DOFs is 11400.The FEM model of the structural cell is given in Figure 8. SPIM, the Newmark method (i.e., the average acceleration implicit scheme) and the variable step R-K method are used for this example in the time interval [0, 0.02] s.The time step for SPIM is  = 10 −5 s, while  = 1, 2, and 3 are used to study the effect of the number of interpolation points on   the accuracy of SPIM.Four time steps  = 10 −5 s, 2 × 10 −5 s, 3 × 10 −5 s, and 4 × 10 −5 s are selected to determine the order of the SPIM.For the Newmark method, two time steps of 10 −6 s and 10 −7 s are used, and 10 −13 is used for both the absolute and relative tolerances in the R-K method (ODE45).The numerical results computed using ODE45 are used as the reference solution, and the relative errors of the SPIM and Newmark methods are defined by (50).present the velocity distributions.The relative errors of the displacements and velocities computed using the SPIM and Newmark methods are given by Figures 10(a) and 10(b), respectively.Figure 10 shows the variations of the errors of the displacement vectors versus the time step.The order of accuracy and CPU time for all methods are given in Table 2.
In terms of Figures 10(a) and 10(b) and Table 2, it can be observed that the SPIM with  = 3 has the best precision and  requires most amount of the CPU time.Hence, for the proposed SPIM, the precision is improved, but the computational efficiency is reduced with an increase in the interpolation points.From Figure 11 and Table 2, it can be concluded that the order of accuracy of the SPIM is approximately 2( + 1).However, compared with the Newmark and R-K methods, the computational efficiency of SPIM is improved, according to Table 2.The CPU time of the SPIM with  = 3 is 321 s, which is much less than that of the R-K method or the Newmark method with the time step of 10 −7 s.Although the CPU time of the SPIM with  = 3 is slightly greater than that of the Newmark method with the time step of 10 −6 s, the SPIM results are much more accurate than the Newmark method, as shown by Figures 10(a) and 10(b).It can also be observed that when  = 1, the precision of SPIM is already close to that obtained with the Newmark method with the time step of 10 −7 s.However, note that when  = 1, the CPU time of the proposed method is approximately 0.015 times that obtained with the Newmark method with the time step of 10 −7 s.

Conclusions
A subdomain precise integration method (SPIM) has been developed to solve the dynamical responses of periodic

Figure 3 :Figure 4 :
Figure 3: The displacements and velocities of all nodes at  = 1000 s.

Figure 5 :
Figure 5: The variations of the errors of the displacement vectors versus the time step for Example 1.

Figure 8 :
Figure 8: The finite element model of the structural cell.

Figures 9 (
Figures 9(a) and 9(b) present the displacement distributions at 0.02 s in the  and  directions computed using the proposed method, respectively, and Figures9(c) and 9(d) present the velocity distributions.The relative errors of the displacements and velocities computed using the SPIM and Newmark methods are given by Figures10(a) and 10(b), respectively.Figure10shows the variations of the errors of the displacement vectors versus the time step.The order of accuracy and CPU time for all methods are given in Table2.In terms of Figures10(a) and 10(b) and Table2, it can be observed that the SPIM with  = 3 has the best precision and Displacement in the x direction (b)Displacement in the y direction (c)Velocities in the x direction (d)Velocities in the y direction

Table 1 :
Comparison of CPU times for Example 1.

Table 2 :
Comparison of CPU times for Example 2.