# Subdomain Precise Integration Method for Periodic Structures

Authors: F. Wu; Q. Gao; W. X. Zhong
Journal: Shock and Vibration
 (2014)

---

## Abstract

A subdomain precise integration method is developed for the dynamical responses of periodic structures comprising many 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.


---

## Body

## 1. 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 lumped-mass 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 elastic-piezoelectric 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–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–26], and the repeatability of periodic structures. This method reduces the computational efforts and reduces the amount of memory required during the computation.


## 2. 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 Figure1) produced using the finite element method (FEM) can be written as
(1)MX¨+CX˙+KX=f
with the initial conditions
(2)Xt0=X0,X˙t0=V0,
where M, C, and K are the N×N mass and damping and stiffness matrices of the periodic structure, respectively, and N is the number of degrees of freedom (DOFs) of the periodic structure. X, X˙, and X¨ are the N×1 displacement and velocity and acceleration vectors, respectively, and the dot over the variable denotes differentiation with respect to time t. f∈RN×1 is the load vector, and t0 denotes the initial time.Figure 1
A periodic structure.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:
(3)U˙=HU+Nf,
where
(4)U=XV,H=0I-M-1K-M-1C,N=0M-1.
In (4), V=X˙ denotes the velocity vector, I∈RN×N is the unit matrix, and 0∈RN×N is a zero matrix.For numerical integration, the time domain is divided into a series of time intervals of equal lengthη; that is,
(5)t0=0,t1=η,…,tk=kη,tk+1=k+1η,….
By letting
(6)Uk=Utk,
the solution to (3) at discrete times is given by
(7)Uk+1=TUk+∫0ηexp⁡Hη-ξNftk+ξdξ,
in which
(8)T=exp⁡Hη.The matrixT 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 ftk+ξ can first be approximated using the Lagrange interpolation polynomial [27] in the time interval tk,tk+1. Hence, if we select m+1 interpolation points, denoted by  tk+ξi, in the time interval tk,tk+1, the load vector can be approximated as
(9)ftk+ξ=∑i=0mLm,iξftk+ξi
for which
(10)Liξ=ξ0-ξξ1-ξ⋯ξi-1-ξξi+1-ξ⋯ξm-ξξ0-ξiξ1-ξi⋯ξi-1-ξiξi+1-ξi⋯ξm-ξi,
where Li is the ith Lagrange interpolation polynomial, i=0,1,…,m. In terms of (7) and (9), the state vector at tk+1 is obtained:
(11)Uk+1=TUk+Ψ0f0,k+Ψ1f1,k+⋯+Ψmfm,k
for which
(12)fi,k=ftk+ξi,Ψi=∫0ηLiξexp⁡⁡Hη-ξNdξ
and Ψi 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 matrixH and time step η can be given as follows. Let H1=Hη/2Q, in which Q is an integer; (8) then becomes
(13)T=exp⁡H12Q.
If Q is large enough, such that the infinite norm of H1 satisfies H1∞≪1, then exp⁡H1 can be approximated by the Taylor series of order g; that is,
(14)exp⁡H1≈I+H1+⋯+H1gg!.
Theoretically, if both Q and g are large enough, (14) yields a very accurate approximation of exp⁡H1. 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⁡H1 is divided into two parts; that is,
(15)exp⁡H1=I+R0,R0=H1+⋯+H1gg!.
Then, using (13), (14), and (15), T can be written as
(16)T=I+R02Q=I+Ri-1I+Ri-12Q-i=I+Ri2Q-i,
where
(17)Ri=2Ri-1+Ri-12,1≤i≤Q.
After computing RQ, T can be given by
(18)T=I+RQ.PIM is an accurate algorithm for computing the matrix exponential. However, the computational effort of PIM isON3, 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.


## 3. The Subdomain PIM

Suppose that a periodic structure comprisesL identical structural cells, and the mass, damping, and stiffness of the ith identical structural cell are Mi, Ci, and Ki, respectively. The number of the DOFs of the ith structural cell is Ni, and the displacement and load vectors are denoted by Xi and fi, respectively. The dynamical model of the ith structural cell is given by
(19)MiX¨i+CiX˙i+KiXi=fi.The DOFs of each structural cell are contained within the DOFs of the whole periodic structure. Hence, the displacement vectorXi of the ith structural cell can be found in the displacement vector X of the whole structure; that is,
(20)Xi=DiTX,
where DiT∈RNi×N is a matrix whose elements are 0 or 1. DiT  describes the relationship between Xi and X and satisfies
(21)DiTDi=Ii
in which Ii∈RNi×Ni is a unit matrix. Multiplying both sides of (19) by Di and substituting (20) into (19), the dynamical equation of the ith structural cell can be rewritten as
(22)DiMiDiTX¨+DiCiDiTX˙+DiKiDiTX=Difi.
In terms of (22), the entire dynamical equation for the whole periodic structure can be written as
(23)∑i=1LDiMiDiTX¨+∑i=1LDiCiDiTX˙+∑i=1LDiKiDiTX=∑i=1LDifi,
where
(24)∑i=1LDiMiDiT=M,∑i=1LDiCiDiT=C,(25)∑i=1LDiKiDiT=K,∑i=1LDifi=f.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 fi. In (19), the force vector fi actually contains two parts, denoted by pi and qi. pi is the force produced by outside factors, such as gravity, and qi is the force produced by the neighboring structural cells. For the whole periodic structure, qi is the so-called internal force, which means
(26)∑i=1LDipi=f,∑i=1LDiqi=0.
The force vectors ft, pit, and qit can also be approximated as
(27)qitk+ξ=∑n=02q-1L2q,nξqi,n,pitk+ξ=∑n=02q-1L2q,nξpi,nfitk+ξ=∑n=02q-1L2q,nξfi,n,fi,n=qi,n+pi,n
by selecting 2q Lagrange interpolating points [27] in the time interval [tk,tk+1]. In (27), pi,n, qi,n, and fi,n are the Ni×1 force vectors at the nth interpolating time point, and qi,n is an undetermined vector. Substituting (27) into (19) yields
(28)MiX¨i+CiX˙i+KiXi=∑n=02q-1L2q,nξfi,n
which can also be rewritten as
(29)U˙i=HiUi+Ni∑n=02q-1L2q,nξfi,n,(30)Ui=XiVi,Ni=0Mi-1,Hi=0Ii-Mi-1Ki-Mi-1Ci.Note thatpi is known and qi is undetermined. Selecting 2q interpolating points to approximate the force vector qi implies that there are 2q unknown vectors (i.e., qi,n,  n = 1~2q) that should be determined. To determine these vectors, take l derivatives of (29) with respect to time:
(31)Uil+1=HiUil+Ni∑n=02q-1L2q,nlξfi,n,
where
(32)Uil=dlUitdtl,l=0,1,…,q-1.
Then, applying PIM to (31) yields
(33)Ui,k+1l=TiUi,kl+∑n=02q-1ψi,l,nfi,n,
where ψi,l,n∈R2Ni×Ni is given by
(34)ψi,l,n=∫0ηexp⁡⁡Hiη-ξNiL2q,nlξdξ
which can be computed using the Gaussian quadrature method [28]. In terms of (33), we have
(35)Ui,k+10Ui,k+11⋮Ui,k+1q-1=TiUi,k0TiUi,k1⋮TiUi,kq-1+ψi,0,0ψi,0,1⋯ψi,0,2q-1ψi,1,0ψi,1,1⋯ψi,1,2q-1⋮⋮⋱⋮ψi,q-1,0ψi,q-1,1⋯ψi,q-1,2q-1fi,0fi,1⋮fi,2q-1.
For convenience, let (35) be rewritten as
(36)Wi,k+1=bi,k+YiQi
in which
(37)Wi,k+1=Ui,k+10Ui,k+11⋮Ui,k+1q-1,Yi=ψi,0,0ψi,0,1⋯ψi,0,2q-1ψi,1,0ψi,1,1⋯ψi,1,2q-1⋮⋮⋱⋮ψi,q-1,0ψi,q-1,1⋯ψi,q-1,2q-1,Qi=qi,0qi,1⋮qi,2q,bi,k=TiUi,k0TiUi,k1⋮TiUi,kq-1+ψi,0,0ψi,0,1⋯ψi,0,2q-1ψi,1,0ψi,1,1⋯ψi,1,2q-1⋮⋮⋱⋮ψi,q-1,0ψi,q-1,1⋯ψi,q-1,2q-1pi,0pi,1⋮pi,2q,
where Yi∈R2qNi×2qNi is a square matrix and pi,n∈RNi×1 is the force vector produced by the neighboring structural cells at the nth interpolating point. Therefore, in the vector pi,n, the elements corresponding to the DOFs inside the structural cell are all zeros. For the ith structural cell, the displacements and velocities of the DOFs inside the cell are denoted by Xi,i and Vi,i, respectively. The displacements and velocities of the DOFs that are on the boundary of the cell are denoted by Xb,i and Vb,i, respectively. The subscripts i and b indicate internal and boundary DOFs, respectively. Meanwhile, let the state vector Ui,k+1 be partitioned into
(38)Ui,k+1=Ui,i,k+1Ub,i,k+1,Ui,i,k+1=Xi,iVi,i,Ub,i,k+1=Xb,iVb,i
and rewrite (36) in the block matrix form
(39)Wi,i,k+1Wb,i,k+1=bi,i,kbb,i,k+Yii,iYib,iYbi,iYbb,i0Qb,i
in which Wi,i,k+1 corresponds to the DOFs inside the structural cell, and Wb,i,k+1 corresponds to the DOFs on the boundary; that is,
(40)Wi,i,k+1=Ui,i,k+10Ui,i,k+11⋮Ui,i,k+1q-1,Wb,i,k+1=Ub,i,k+10Ub,i,k+11⋮Ub,i,k+1q-1,Qb,i=qb,i,0qb,i,1⋮qb,i,2q.
In terms of (39), we have
(41)Ybb,i-1Wb,i,k+1=Ybb,i-1bb,i,k+Qb,i,(42)Wi,i,k+1=bi,i,k+Yib,iYbb,i-1Wb,i,k+1-bb,i,k.Equation (41) can be seen as the dynamical stiffness equation of the ith structural cell, and Qb,i 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 ith structural cell is reduced, and the ith 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 Xb and Vb, respectively, and we have
(43)Xb,i=Db,iTXb,Vb,i=Db,iTVb,
where Db,iT describes the relationship between Xb,i and Xb and Db,iTDb,i is a unit matrix. Combining (38) and (40) with (43) yields
(44)Ub,i,k+1l=Eb,iTUb,k+1l,Wb,i,k+1=Gb,iTWb,k+1,
where
(45)Eb,iT=Db,iT00Db,iT,Ub,k+1l=Xb,k+1lVb,k+1l,Ub,i,k+1l=Xb,i,k+1lVb,i,k+1l,Wb,k+1=Ub,k+10Ub,k+11⋮Ub,k+1q-1,Gb,iT=EbTEbT⋱EbT.
Substituting (44) into (41) and multiplying (41) with Gb,i yield
(46)Gb,iYbb,i-1Gb,iTWb,k+1=Gb,iYbb,i-1bb,i,k+Gb,iQb,i.
In terms of (46), we have
(47)∑i=1LGb,iYbb,i-1Gb,iTWb,k+1=∑i=1LGb,iYbb,i-1bb,i,k+∑i=1LGb,iQb,i.Qb,i 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
(48)∑i=1LGb,iQb,i=0.
Combining (47) with (48) yields
(49)∑i=1LGb,iYbb,i-1Gb,iTWb,k+1=∑i=1LGb,iYbb,i-1bb,i,k.In terms of (49), Wb,k+1 can be computed. Once Wb,k+1 is obtained, the vector Wb,i,k+1 can be computed using (44). Then, Wi,i,k+1 can also be computed by applying Wb,i,k+1 to (42). When Wb,i,k+1 and Wi,i,k+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 [tk,tk+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 number2q of the interpolation points for approximating qit 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 q on the performance of SPIM will be discussed using two numerical examples.


## 4. Numerical Examples

Example 1.
Consider the periodic structure consisting of masses and springs depicted in Figure2. The entire periodic structure is composed of 50 structural cells, shown in Figure 2(a), with its masses and stiffnesses being m1=1 (kg) and m2=2 (kg) and k=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.The spring-mass model.
(a)
(b)To study the effect of the number2q of interpolation points on the accuracy of SPIM, q=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
(50)ed=log10X-XRK2XRK2,ev=log10V-VRK2VRK2,
where XRK and VRK 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 iss, the numerical error for time step η is Cηs, in which C is a constant. Thus, the errors e1 and e2 for two different time steps η1 and η2 can be given by
(51)e1=Cη1s,e2=Cη2s
and the order of the algorithm can be shown by
(52)s=log⁡10e1-log⁡10e1log⁡10η1-log⁡η2.In this example, the errors of the displacement vectors are defined by (50); thus, the order of the proposed method is defined as
(53)s=ed,1-ed,2log⁡10η1-log⁡η2.
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.Figures3(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.Table 1
Comparison of CPU times for Example1.
SPIM(q=1)
SPIM(q=2)
ODE45
Newmark
Order
4
5
—
2
η  (s)
1
0
-
1
1
0
-
1
—
1
0
-
3
1
0
-
4
CPU times (s)
35
42
374
317
3111The displacements and velocities of all nodes att=1000 s.
(a)
Displacements of the nodes att=1000 s
(b)
Velocities of the nodes att=1000 sThe relative errors for Example1. SPIM with q=2 and time step 10-1s (solid line), SPIM with q=1 and time step 10-1s (•), Newmark method with time step 10-4 s (□), and Newmark method with time step 10-3 s (△).
(a)
Relative errors of the displacement vectors
(b)
Relative errors of the velocity vectorsFigure 5
The variations of the errors of the displacement vectors versus the time step for Example1.Figures4(a) and 4(b) demonstrate that the results obtained using the SPIM with q=2 are more accurate than those obtained using the SPIM with q=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 2q+1. However, in terms of Table 1, the SPIM with q=2 requires slightly more CPU time than the SPIM with q=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 q=1 and q=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 q=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.Example 2.
Consider a two-dimensional phononic crystal consisting of two different materials. The corresponding material parameters areE1=1.4×1010 Pa, E2=7.0×1010 Pa, ρ1=2000 kg/m3, ρ2=2700 kg/m3, ν1=0.3, and ν2=0.3. 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 ft=2×107 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.Figure 6
The structural cell.Figure 7
The phononic crystal.Figure 8
The finite element model of the structural cell.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 interval0,0.02 s. The time step for SPIM is η=10-5 s, while q=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).Figures9(a) and 9(b) present the displacement distributions at 0.02 s in the x and y directions computed using the proposed method, respectively, and Figures 9(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 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.Table 2
Comparison of CPU times for Example2.
SPIM(q=1)
SPIM(q=2)
SPIM(q=3)
Newmark
ODE45
Order
4
6
8
2
—
η (s)
1
0
-
5
1
0
-
5
1
0
-
5
1
0
-
6
1
0
-
7
—
CPU times (s)
68
169
321
291
4779
2338Figure 9
The displacement and velocity distributions att=0.02 s.The relative errors for Example2. SPIM with q=3 and time step 10-5 s (solid line), SPIM with q=2 and time step 10-5 s (•), SPIM with q=1 and time step 10-5 s (□), Newmark method with time step 10-7 s (○), and Newmark method with time step 10-6 s (△).
(a)
Relative errors of the displacement vectors
(b)
Relative errors of the velocity vectorsIn terms of Figures10(a) and 10(b) and Table 2, it can be observed that the SPIM with q=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 2q+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 q=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 q=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 q=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 q=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.Figure 11
The variations of the errors of the displacement vectors versus the time step for Example2.


## 5. Conclusions

A subdomain precise integration method (SPIM) has been developed to solve the dynamical responses of periodic structures comprising many identical structural cells. SPIM does not require the entire mass, damping, and stiffness matrices of the periodic structure. Using the original precise integration method for the structural cell and considering the repeatability of periodic structures, the method avoids repeatedly computing and storing many involved identical exponential matrices. Hence, the presented method not only inherits the precision of the original method but also improves the computational efficiency in terms of both CPU time and storage. Numerical examples demonstrate that the proposed method is more efficient than the R-K method or Newmark method and is more precise than the Newmark method. It can be concluded that the proposed method is particularly suitable for periodic structures containing many identical structural cells.



---

*Source: 657589-2014-11-09.xml*