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.

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 Figure 1) produced using the finite element method (FEM) can be written as with the initial conditions where , , and 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. , , and are the displacement and velocity and acceleration vectors, respectively, and the dot over the variable denotes differentiation with respect to time . is the load vector, and 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: where In (4), denotes the velocity vector, is the unit matrix, and is a zero matrix.

For numerical integration, the time domain is divided into a series of time intervals of equal length ; that is, By letting the solution to (3) at discrete times is given by in which

The matrix defined by (8) is called the exponential of the matrix . To evaluate the second term on the right-hand side of (7), the load vector can first be approximated using the Lagrange interpolation polynomial [27] in the time interval . Hence, if we select interpolation points, denoted by  , in the time interval , the load vector can be approximated as for which where is the th Lagrange interpolation polynomial, . In terms of (7) and (9), the state vector at is obtained: for which 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 to improve the computational precision.

The implementation of PIM for computing the exponential of the matrix and time step can be given as follows. Let , in which is an integer; (8) then becomes If is large enough, such that the infinite norm of satisfies , then 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 . However, the round-off error may be significant in numerical computation because of the sum of the unit matrix and the small rest terms. Therefore, in PIM, the exponential matrix is divided into two parts; that is, Then, using (13), (14), and (15), can be written as where After computing , can be given by

PIM is an accurate algorithm for computing the matrix exponential. However, the computational effort of PIM is , 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 comprises identical structural cells, and the mass, damping, and stiffness of the th identical structural cell are , , and , respectively. The number of the DOFs of the th structural cell is , and the displacement and load vectors are denoted by and , 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 of the th structural cell can be found in the displacement vector of the whole structure; that is, where is a matrix whose elements are 0 or 1.   describes the relationship between and and satisfies in which is a unit matrix. Multiplying both sides of (19) by 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 . In (19), the force vector actually contains two parts, denoted by and . is the force produced by outside factors, such as gravity, and is the force produced by the neighboring structural cells. For the whole periodic structure, is the so-called internal force, which means The force vectors , , and can also be approximated as by selecting Lagrange interpolating points [27] in the time interval . In (27), , , and are the force vectors at the th interpolating time point, and is an undetermined vector. Substituting (27) into (19) yields which can also be rewritten as

Note that is known and is undetermined. Selecting interpolating points to approximate the force vector implies that there are unknown vectors (i.e., ,   = 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 is given by which can be computed using the Gaussian quadrature method [28]. In terms of (33), we have For convenience, let (35) be rewritten as in which where is a square matrix and is the force vector produced by the neighboring structural cells at the th interpolating point. Therefore, in the vector , 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 and , respectively. The displacements and velocities of the DOFs that are on the boundary of the cell are denoted by and , respectively. The subscripts and indicate internal and boundary DOFs, respectively. Meanwhile, let the state vector be partitioned into and rewrite (36) in the block matrix form in which corresponds to the DOFs inside the structural cell, and corresponds to the DOFs on the boundary; that is, In terms of (39), we have

Equation (41) can be seen as the dynamical stiffness equation of the th structural cell, and 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 and , respectively, and we have where describes the relationship between and and is a unit matrix. Combining (38) and (40) with (43) yields where Substituting (44) into (41) and multiplying (41) with yield In terms of (46), we have 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), can be computed. Once is obtained, the vector can be computed using (44). Then, can also be computed by applying to (42). When and 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 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 of the interpolation points for approximating 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.

4. 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  (kg) and  (kg) and  (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 th 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.

(a)

(b)

(a)
(b)

Figure 2 

The spring-mass model.

To study the effect of the number of interpolation points on the accuracy of SPIM, and 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  s. The time step for SPIM is  s, two time steps of  s and  s are used for the Newmark method, and 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 and are the displacement and velocity vectors computed using the R-K method, respectively. and 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 and for two different time steps and 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  s,  s,  s, and  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  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.

(a) Displacements of the nodes at  s

(b) Velocities of the nodes at  s

(a) Displacements of the nodes at  s
(b) Velocities of the nodes at  s

Figure 3 

The displacements and velocities of all nodes at  s.

(a) Relative errors of the displacement vectors

(b) Relative errors of the velocity vectors

(a) Relative errors of the displacement vectors
(b) Relative errors of the velocity vectors

Figure 4 

The relative errors for Example 1. SPIM with and time step (solid line), SPIM with and time step (), Newmark method with time step  s (□), and Newmark method with time step  s ().

Figure 5 

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

Figures 4(a) and 4(b) demonstrate that the results obtained using the SPIM with are more accurate than those obtained using the SPIM with when the time steps used for the two cases are both  s. From Figure 5 and Table 1, it can be concluded that the order of accuracy of the SPIM is approximately . However, in terms of Table 1, the SPIM with requires slightly more CPU time than the SPIM with . 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 and , 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 , 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 are  Pa,  Pa,  kg/m³,  kg/m³, , and . 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 structural cells. The bottom of the phononic crystal is fixed. On top of the phononic crystal acts the uniform load  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 interval s. The time step for SPIM is  s, while , 2, and 3 are used to study the effect of the number of interpolation points on the accuracy of SPIM. Four time steps  s,  s,  s, and  s are selected to determine the order of the SPIM. For the Newmark method, two time steps of  s and  s are used, and 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).

Figures 9(a) and 9(b) present the displacement distributions at  s in the and 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.