An efficient algorithm based on group theory and the Woodbury formula for the dynamic responses of periodic structures

Computers and Structures 182 (2017) 238–251

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

An efficient algorithm based on group theory and the Woodbury formula for the dynamic responses of periodic structures X.Q. Liang, Q. Gao ⇑, W.A. Yao State Key Laboratory of Structural Analysis for Industrial Equipment, International Center for Computational Mechanics, Department of Engineering Mechanics, Faculty of Vehicle Engineering and Mechanics, Dalian University of Technology, Dalian 116023, PR China

a r t i c l e

i n f o

Article history: Received 18 August 2016 Accepted 5 December 2016

Keywords: Periodic structures Woodbury formula Group theory Newmark method Dynamic responses

a b s t r a c t An efficient numerical method for computing the dynamic responses of periodic structures is proposed. Efficiently solving a system of linear equations is a key issue for computing the dynamic response of a large-scale structure. Based on the periodic properties of a structure and using condensation technology, the system of linear equations that corresponds to the periodic structure is reduced to a small-scale system of linear equations. Based on the Woodbury formula, the solution for the small-scale system of linear equations is obtained by solving a new system of linear equations whose coefficient matrix corresponds to a cycle periodic structure. Using group theory, the new system of linear equations is efficiently solved. Superior efficiency is achieved because the scale of the system of linear equations for the entire structure is significantly reduced and the coefficient matrix for the new system of linear equations can be converted into a block-diagonal matrix using group theory. Numerical examples are presented to illustrate the high efficiency of the proposed method. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction A periodic structure consists of a number of identical structural components (unit cells) that are repeated in one, two or three directions to form a complete system, such as a beam lattice [1,2], periodic beam structures [3,4], cellular solids [5], periodic composite structures [6–8], and photonic [9] and phononic [10,11] crystals. The theoretical studies and applications of the periodic structures have aroused substantial concern among research institutions and scholars all over the world. The wave motion and propagation are very important aspects for the periodic structures. In this area, it is mainly concentrated on the dispersion relations and the scattering analysis of the periodic structures. The dispersion relation has attracted the attention of many scholars. Tassilly [12] derived the dispersion relation between wavenumber and frequency for a class of non-uniform periodic elements. Gavric [13] used the concepts of cross-section modes and finite element method to calculate the dispersion relation for an free rail. Using the transfer matrix method and Floquet’s theorem, Shen and Cao [14] derived a dispersion relation for acoustic wave propagation in a periodic layered structure. For the infi⇑ Corresponding author. E-mail addresses: [email protected] (X.Q. Liang), [email protected] (Q. Gao), [email protected] (W.A. Yao). http://dx.doi.org/10.1016/j.compstruc.2016.12.002 0045-7949/Ó 2016 Elsevier Ltd. All rights reserved.

nite two-dimensional periodic lattices, the dispersion curves were obtained by solving the eigenvalue problem for wave propagation [15]. By studying the wave behavior in periodically simply supported beam and periodic frame structures, Sonekar and Mitra [16] analyzed the dispersion relation for the periodic structures. Through analysis of the elastic wave propagation in onedimensional periodic elastic rod structure and two-dimensional periodic elastic beam structure, Tian et al. [17] derived the dispersion relation between Bloch wave vectors and eigenfrequencies. Trainiti et al. [18] investigated the effects of periodic geometric undulations on the dispersion properties of one-dimensional and two-dimensional elastic structures. Meanwhile, the scattering analysis has also received much attention. Vonflotow [19] showed the scattering behavior of the junction of the Timoshenko bending model and the periodic torsion model. Cai and Lin [20] analyzed wave scattering in generic structural networks. By studying the propagation properties of flexural wave in the periodic beam on elastic foundations, Yu et al. [21] designed and calculated Bragg scattering gap. For analyzing the scattering in infinite periodic structures, Petersson and Jin [22] proposed a new twodimensional time domain finite element method formulation. Dynamic vibration is another important aspect for the periodic structures. In this area, it is mainly concentrated on the energy band and dynamic responses of the periodic structures. Because the band gap of periodic structure could be exploited in devices

239

X.Q. Liang et al. / Computers and Structures 182 (2017) 238–251

with widespread application in engineering, it has been received much attention. Schmidt and Kauf [23] computed the band structure of two-dimensional photonic crystals. Olhoff et al. [24] computed the maximizing frequency gaps of beam structures. Wang et al. [11] investigated the wave band gaps in two-dimensional piezoelectric and piezomagnetic phononic crystals. By employing the Bloch-Floquet theorem, Xiang and Shi [25] determined the flexural vibration band gaps in periodic beams. The vibration band-gap characteristics of periodic rigid frame structures composed of Timoshenko beams [26] and periodic Mindlin plate structure with two simply supported opposite edges [27] were studied. Meanwhile, solving the dynamic responses also have been attracted many attentions. Engels [28] investigated the dynamic responses of infinite and semi-infinite periodic structures to harmonic loads. Wu et al. [29] analyzed the dynamic behavior of periodic plate structures. Cai et al. [30] applied a U-transformation method to obtain dynamic solutions for periodic structures. Zhou et al. proposed efficient numerical approaches for studying the vibration of complex one-dimensional [31] and two-dimensional [32] periodic structures. Hawreliak et al. [33] investigated the dynamic responses of additively manufactured engineered lattice materials. Luongo and Romeo [34] presented a modified version of the traditional wave vector computational scheme for the dynamic analysis of long undamped periodic structures. Using a combination of wave and finite element approaches, Duhamel et al. [35] presented a method for calculating the force responses of periodic structures. Mencik [36] investigated a wave finite element method for computing the low- and mid-frequency forced responses of straight elastic structures. Based on a precise integration method, Gao et al. proposed efficient algorithms for computing the dynamic responses of one-dimensional [37] and two-dimensional [38] periodic structures. Gao et al. [39] derived an exact analytical solution for the dynamic response of a periodic structure for which the unit cell consists of one mass and one spring. Wu et al. [40] developed a subdomain precise integration method for the dynamic responses of periodic structures. For large scale periodic structures, improving the computational efficiency of the numerical methods for solving the dynamic responses is an essential issue. It is well known that the group theory is an effective tool for improving the computational efficiency for symmetrical structures and so has been widely used in science and engineering [41–48]. Using the group-theoretic methods and substructuring technique, Zhong and Qiu [41] analyzed the symmetric or partially symmetric structures subjected to arbitrary loads. Based on the symmetry groups and representation theory, Zingoni [42] proposed an efficient computational scheme for the linear vibration analysis of high tension cable nets. Based on the full symmetry properties of a structure, group representation theory and Fourier method, the full stiffness matrix of a structure was converted into the form of block diagonalization [43], and the similar techniques was used to block-diagonalise the equilibrium matrix of a symmetric structure [44]. The unsymmetrical (or weakly symmetrical) spring-mass dynamic systems was transformed into equivalent dynamic systems featuring maximum possible symmetry, and then the group theoretic procedure employed to calculate all the eigenvalues of the systems [45]. Using the group theory, the eigenvalue problems were decomposed into independent subproblems due to the block diagonalized [46,47]. Zingoni [48] showed that the group-theoretic approach enables considerable simplifications and reductions in computational effort. However, although the infinite periodic structure is translational symmetry, but the periodic structure composed of finite unit cells has no symmetric property, so the group theory cannot be applied directly to the finite periodic structures. This paper proposed an efficient algorithm for solving the dynamic responses of one-dimensional and two-dimensional peri-

odic structures. The basic idea of the numerical algorithms for computing the dynamic responses is to transform the dynamic equations into a system of linear equations. Therefore, efficiently solving the system of linear equations is a key issue for computing the dynamic responses. Firstly, based on the periodic properties of a structure and using condensation technology, the system of linear equations for the entire periodic structure is reduced to a small-scale system of linear equations in this paper. Then, based on the Woodbury formula, the solution of the small-scale system of linear equations is obtained by solving a new system of linear equations whose coefficient matrix corresponds to a cycle periodic structure. Finally, using group theory, the new system of linear equations is efficiently solved. Using the three technologies enables an efficient algorithm to be proposed for solving the system of linear equations. This algorithm leads to an efficient method for computing the dynamic responses of periodic structures.

2. The basic equations Assume that a structure is considered and that its stiffness, mass and damping matrices are K, M and C， respectively. The dynamic equation can be written as

€ þ Cu_ þ Ku ¼ f Mu

ð1Þ

€ denote the where f denotes the external force vector, and u, u_ and u displacement vector, velocity vector and acceleration vector, respectively. Many numerical algorithms, such as the Newmark method [49], the generalized a method [50], the Bathe method [51,52], and the central difference method [53], can be applied for solving Eq. (1). For all methods, the basic idea is to transform Eq. (1) into the following system of linear equations

 t þDt ¼ f Ku 0

ð2Þ

 and f denote the equivalent stiffness matrix and the equivwhere K alent external force vector. Solving Eq. (2) gives the dynamic responses of Eq. (1). The main computational effort for solving the structural dynamic responses comprises repeatedly performing the system of linear Eq. (2). For periodic structures, the degrees of freedom (DOF) of the entire structure is significant when it contains many identical unit cells. Thus, it is very time-consuming for solving the system of linear equations. Based on the characteristics of a periodic structure, an efficient method is proposed in this paper to solve the system of linear equations (2) using the Woodbury formula and group theory. Because all previously mentioned numerical methods are required to solve Eq. (2), the Newmark method is used as an example in this paper to explain the idea of the proposed method. For 

the Newmark method, the equivalent stiffness matrix K and equiv

alent external force vector f are [49] 

K ¼ K þ c 0 M þ c1 C f ¼ f t þDt þ Mðc0 ut þ c2 u_ t þ c3 u € t0 Þ þ Cðc1 ut0 þ c4 u_ t0 þ c5 u € t0 Þ 0 0 0 ð3Þ € t0 are the displacement vector, velocity vector where ut0 , u_ t0 and u and acceleration vector at t0 , f t0 þDt denotes the external force vector at t 0 þ Dt and c0 –c5 in Eq. (3) are the parameters of the Newmark method [49]. Then, the displacement vector ut0 þDt at t 0 þ Dt can be obtained by solving Eq. (2), and the velocity and acceleration vectors at t 0 þ Dt can be obtained [49].

240

X.Q. Liang et al. / Computers and Structures 182 (2017) 238–251

The detailed procedures of the proposed method for onedimensional and two-dimensional periodic structures are given in Sections 3 and 4, respectively.

boundary of the entire structure. When the boundary conditions are applied, the DOFs of u1L and uNþ1 may be changed. R 3.1. The application of condensation technology

3. The efficient method for one-dimensional periodic structures This section discusses a finite one-dimensional periodic structure that consists of N þ 1 unit cells. The arbitrary boundary conditions are applied on the left-hand and right-hand ends of the onedimensional periodic structure, as shown in Fig. 1. For the n-th unit cell, as shown in Fig. 2, assume that unL , unR and unI denote the displacement vectors that correspond to the DOFs on its left boundary, right boundary and internal DOFs, respectively, n n n that f L , f R and f I denote the external force vectors that correspond to the DOFs on its left boundary, right boundary and internal DOFs, respectively, and that pnL and pnR denote the internal force vectors that correspond to the DOFs on its left boundary and right boundary, respectively. The number of DOFs on its left and right boundaries is p. Based on these notations, the displacement and external force vectors, which correspond to the n-th unit cell, can be written in block form as


un ¼ ðunL ÞT

ðunI ÞT

h n f ¼ ðf nL ÞT

n T

ðf I Þ

ðunR ÞT n T

ðf R Þ

T

iT

ð4Þ ð5Þ

It’s clear that

unR ¼ unþ1 L

n ¼ 1; 2;    ; N

ð6Þ

According to the form of the displacement vectors in Eq. (4), the stiffness matrix Ke , the mass matrix Me and the damping matrix Ce of the unit cell can also be written in block form. Then, the displacement vector u and external force vector f for the entire structure can be given according to Eqs. (4)–(6), and the stiffness, mass and damping matrices for the entire structure can be obtained by the superposition of the stiffness, mass and damping matrices for all unit cells. Note that u1L and uNþ1 are the disR placement vectors that correspond to the left boundary and right

When a structure includes many unit cells, the scale of Eq. (2) is very large. Thus, solving Eq. (2) is time-consuming. To improve the computational efficiency, the condensation technology is used to reduce the size of Eq. (2), based on the periodic properties of the structure. For the one-dimensional periodic structure in Fig. 1, Eq. (2) can be described as follows: First, the stiffness, mass, damping matrices and force vector of the entire structure are formed. Second, the boundary constraints are applied. Last, the equivalent stiffness matrix and equivalent force vector of the entire structure can be obtained. Eq. (2) can also be obtained in the following manner. First, in a unit cell, the equivalent stiffness matrix and equivalent force vector are computed from the stiffness, mass, damping matrices and the force vector. Second, they are superposed to form the equivalent stiffness matrix and equivalent force that correspond to the entire structure. Last, the boundary conditions are applied to the structure. The benefit of using the latter method is that we can condensate the internal DOFs in a unit cell; then, a small system of linear equation can be obtained for the entire structure. Thus, the system of linear equation can be solved more efficiently. Because every unit cell in the periodic structure is identical, only one calculation is needed. The specific process is given as follows:  e and equivalent force vector f n The equivalent stiffness matrix K of the n-th unit cell can be obtained from Eq. (3). Then, the system of linear equations of the n-th unit cell can be expressed as n n  e un K t 0 þDt ¼ f þ p

ð7Þ

unt0 þDt

where denotes the displacement vector of the n-th unit cell at t0 þ Dt, and pn denotes the internal force vector, which can be written in the following block form as


pn ¼ ðpnL ÞT

0 ðpnR ÞT

T

ð8Þ

According to Eqs. (4), (5) and (8), Eq. (7) can be written in the block form, then using the condensation technology, the internal DOFs that corresponds the displacement vector unI can be eliminated, which produces



Fig. 1. One-dimensional periodic structure.

WLL

WLR

WRL

WRR



unL unR



 ¼ t0 þDt

  n pL rnL þ rnR pnR

where the coefficient matrix denotes the equivalent stiffness matrix of condensation unit cell, rnL and rnR denote the equivalent external force vectors of condensation unit cell. Because all unit cells of the entire structure are identical, the computational effort of the condensation is very small. The system of linear equations of the entire structure can be obtained by combining the condensation Eq. (9) of all unit cells. Let uA denote the displacement vector that corresponds to the DOFs on the left-edge and right-edge of the entire structure, and let uB denote the displacement vector that corresponds to remaining DOFs, i.e.,

h i T T T uA ¼ ðu1L Þ ; ðuNþ1 Þ ; R h i T T T T T Þ ; ðuNR Þ uB ¼ ðu1R Þ ; ðu2R Þ ; . . . ; ðuN1 R

Fig. 2. Unit cell for one-dimensional periodic structure.

ð9Þ

ð10Þ

Then, the system of linear equations of the entire structure can be given by

241

X.Q. Liang et al. / Computers and Structures 182 (2017) 238–251



RAA RBA

RAB RBB



uA uB



 ¼

rA rB

 ð11Þ

2W

where

rA ¼

r1L rNþ1 R

# ;

r1R þ r2L

3

7 6 2 6 rR þ r3L 7 7 6 7 6 7 6 .. rB ¼ 6 7 . 7 6 7 6 6 rN1 þ rN 7 4 R L 5

ð12Þ

and RBB is a N  N block matrix in the following form

RBB

6 6 6 6 6 ¼6 6 6 6 6 4

WRR þ WLL

WLR

0

0

0

WRL

WRR þ WLL

WLR

0

0

0

WRL

..

..

0

0

0

..

0

0

0

. .

.

3 7 7 7 7 7 7 7 7 7 7 5

WRR þ WLL

WLR

WRL

WRR þ WLL ð13Þ

Once Eq. (11) is solved, the displacement vectors unL and unR for all unit cells can be obtained, then the displacement vector unI for the internal DOFs can be obtained. Using condensation technology, Eq. (2) is transformed into a smaller system of linear Eq. (11) in this subsection, which can reduce the computational effort.

0

WRL

WRL

WRR þ WLL WRL

0 .. .

0

0

0

0 WLR

0 0

WLR .. . .. . 0

WRR þ WLL WRL

WLR WRR þ WLL

7 7 7 7 7 7 7 5

Comparing the matrix RBB with the matrix R shows that the block matrices WLR and WRL are added in the bottom-left and top-right corners of the matrix R, which forms

RBB ¼ R þ UV  U¼

1 0  0 0 0 0  0 1

ð15Þ

where uA denotes the displacement vector that corresponds to the DOFs of two ends of the structure, and uB denotes the displacement vector that corresponds to the remaining DOFs. The dimension of uA is much smaller than the dimension of uB . Because the matrix RAA  RAB R1 BB RBA is a small matrix, its inverse can be easily com1 1 puted. Thus, computing R1 BB RBA , R BB rB and RBB ðrB  R BA uA Þ is the key issue when solving Eqs. (14) and (15). They can be computed by solving the following equation:

ð16Þ

Note that R1 BB RBA is only computed once. In this subsection, the Woodbury formula is used to transform Eq. (16) into another system of linear equations that can be solved more efficiently. The basic idea of the Woodbury formula is as follows: assume that we want to solve a system of linear equation with the coefficient matrix A but it is time-consuming. If we can obtain the matrix B by making a ‘‘small” change in matrix A and if we have a fast method for solving a system of linear equations with the coefficient matrix B, the Woodbury formula enables us to efficiently solve a system of linear equations with the coefficient matrix A [54].

 0 V¼ WLR

0    0 WRL 0  0

ðR þ UVÞw ¼ b



0

ð20Þ

According to the Woodbury formula [54], the solution of Eq. (20) can be given by

where

uB ¼ R1 BB ðrB  R BA uA Þ

 Ip ;

where  denotes the Kronecker product [55], and Ip is the identity matrix of p dimensions. U is an Np  2p matrix, and V is a 2p  Np matrix. Using Eq. (18), Eq. (16) is equivalent to the following equation:

In the last section, the scale of the system of linear equations is reduced using condensation technology. Solving Eq. (11) remains time-consuming when the structure includes many unit cells. Therefore, based on the Woodbury formula, Eq. (11) is transformed into a system of linear equations, which can be solved more efficiently. Eq. (11) can be solved in the following manner

ð14Þ

T

ð19Þ

w ¼ y  Z½H1 ðVyÞ

  1 RAA  RAB R1 BB R BA uA ¼ rA  R AB R BB rB

ð18Þ

where

3.2. The transformation based on the Woodbury formula

RBB w ¼ b

3

0

ð17Þ

rNR þ rNþ1 L 2

6 6 6 6 R¼6 6 6 4

þ WLL

WLR

RR

2 "

Based on the idea of Woodbury formula, we construct the matrix R by making a ‘‘small” change in RBB , i.e.,

ð21Þ

H ¼ I þ VZ

ð22Þ

RZ ¼ U

ð23Þ

Ry ¼ b

ð24Þ

Therefore, solving Eq. (20) is converted into solving Eqs. (23) and (24) and computing Eqs. (21) and (22). Eqs. (19) and (23) show that the dimensions of the matrix H is 2p  2p; thus, the computational effort of Eq. (21) is small. Therefore, the major computational effort is to solve Eq. (24). If Eq. (24) can be efficiently solved, Eq. (16) can be efficiently solved, and the dynamic responses of the structure can be given by efficiently solving Eqs. (14) and (15). Therefore, efficiently solving Eq. (24) is the key for computing the dynamic responses of the structure. 3.3. The group theory for one-dimensional periodic structures Section 3.2 shows that improving the computational efficiency of Eq. (24) is a key point for solving the dynamic responses for a periodic structure. The structure of the coefficient matrix R in Eq. (24) is given by Eq. (17), which shows that the matrix R has the following properties, namely, the second row of blocks is the permutation of the first row of blocks, the third row of blocks is the permutation of the second row of blocks, and so on, and the first row of blocks is the permutation of the last row of blocks. The special form of the equivalent stiffness matrix R is identical to the special form of the cycle periodic structure. It is well known that group theory can be applied to cycle periodic structures, which can reduce the computational effort and reduce the usage of internal storage. Based on this argument, group theory can also be applied

242

X.Q. Liang et al. / Computers and Structures 182 (2017) 238–251

 1 ¼ WLL þ WRR þ WLR þ WRL R

to Eq. (24) to improve the computational efficiency. In this subsection, we introduce the method. The basic idea of solving Eq. (24) is as follows: based on the symmetric properties of the matrix R, a transformation matrix can be computed using the irreducible representation of a symmetry group, which can transform the matrix R into a block-diagonal form. Because the transformation matrix can be easily obtained and solving a linear equation with a block-diagonal coefficient matrix is more efficient [41], Eq. (24) can be efficiently solved. The details are as follows: Based on the irreducible representation of the C N symmetry group, for any integer N (the number of unit cells), the transformation matrix TN that corresponds to the matrix R can be given as

(
TN ¼




TI

TIII;1

TI

TII

TIV;1 TIII;1

TIII;2 TIV;1

TIV;2 TIII;2

   TIII;ððN1Þ=2Þ TIV;2

TIV;ððN1Þ=2Þ

   TIII;ðN=21Þ

 i ¼ I2  ðWLL þ WRR Þ þ Pi  WLR þ ðPi ÞT  WRL R N N 2; 3; . . . ; ðN þ 1Þ=2; when N is an odd integer i¼ 1; 2; 3; . . . ; N=2; when N is an even integer

ð31Þ

ð32Þ

the coefficient matrix PN has the following form

 P1N

T

TIV;ðN=21Þ

when N is an odd integer

T

¼

1

"



0

0 1

;

PiN

¼

Pi;1;1 N

PNi;1;2

Pi;2;1 N

PNi;2;2

# i ¼ 2; 3; . . .

ð33Þ

when N is an odd integer

ð25Þ

when N is an even integer

in which



1 sinðð2N  2Þði  1ÞhÞ ðN  1Þ cosðði  1ÞhÞ þ 2 cosððN  1Þði  1ÞhÞ þ N 2 sinðði  1ÞhÞ

1 1  cosðð2N  2Þði  1ÞhÞ ¼ ðN  1Þ sinðði  1ÞhÞ þ N 2 sinðði  1ÞhÞ

1 1  cosðð2N  2Þði  1ÞhÞ 2 sinððN  1Þði  1ÞhÞ  ðN  1Þ sinðði  1ÞhÞ þ ¼ N 2 sinðði  1ÞhÞ

1 sinðð2N  2Þði  1ÞhÞ ðN  1Þ cosðði  1ÞhÞ  ¼ N 2 sinðði  1ÞhÞ

P i;1;1 ¼ N P i;1;2 N P i;2;1 N P i;2;2 N

where

TI ¼ TII ¼

h

Then, the solution of Eq. (24) can be given by p1ffiffiffi N

h

p1ffiffiffi N

h TIII;n ¼ 1 h TIV;n ¼ 0

p1ffiffiffi N



1 p ffiffiffi N



pffiffiffiffiffiffi

cosðnhÞ N=2

pffiffiffiffiffiffi

sinðnhÞ N=2

p1ffiffiffi N p1ffiffiffi N

p1ffiffiffi N

iT

1 p ffiffiffi N

 1 ½ðTN  Ip Þbg y ¼ ðTTN  Ip ÞfR

iT

pffiffiffiffiffiffi



cosððN2ÞnhÞ



sinððN2ÞnhÞ

N=2

pffiffiffiffiffiffi

N=2

ð26Þ

i cosððN1ÞnhÞ T pffiffiffiffiffiffi

N=2

i sinððN1ÞnhÞ T pffiffiffiffiffiffi

N=2

easy to verify that TN TTN ¼ I. Based on the transformation matrix TN , we introduce a new vec, which satisfies the following equation tor y

 y ¼ ðTTN  Ip Þy

ð27Þ

Because R is a block-diagonal matrix, Eq. (35) can be efficiently computed. Because the scale of matrix TN  Ip is large, directly computing ðTN  Ip Þb is time-consuming. Fortunately, we can compute ðTN  Ip Þb in the following manner: According to the block form of matrix R in Eq. (17), vector b can be given in block form as


b ¼ ðb1 ÞT

ðb2 Þ

T

T

   ðbN1 Þ

 T T

ð36Þ

ðbN Þ

Then


ðTN  Ip Þb ¼ ðc1 ÞT

ðc2 ÞT

   ðcN1 ÞT

ðcN ÞT

T

ð37Þ

in which

Then Eq. (24) can be written as

y  ¼ ðTN  Ip Þb R

ð28Þ

ci ¼ T i1 b1 þ T i2 b2 þ    þ T iN bN

ð38Þ

Eq. (38) can be written in the following form as

where

2

 ¼ ðTN  Ip ÞRðTT  Ip Þ R N

ð29Þ

Because the transformation matrix TN is obtained by the irreducible 

representation of the symmetry group, we know that R is a blockdiagonal matrix with the following form

 ¼ diagfR  ig R

ð35Þ



in which h ¼ 2p=N, and n ¼ 1; 2; . . . ; ðN  1Þ=2 when N is an odd integer and n ¼ 1; 2; . . . ; N=2  1 when N is an even integer. It is

where

ð34Þ

ð30Þ

½ c1

c2

   cN  ¼ ½ b1

b2

3 T N1 T N2 7 7 .. 7 7 . 5

T 11 6T 6 12    bN 6 6 .. 4 .

T 21 T 22 .. .

  .. .

T 1N

T 2N

   T NN ð39Þ

According to Eq. (37), the value of ðTN  Ip Þb can be obtained by Eq. (39).

X.Q. Liang et al. / Computers and Structures 182 (2017) 238–251

This subsection explains the application of group theory in onedimensional periodic structures. Based on the irreducible representation of the symmetry group, the coefficient matrix in Eq. (24) is converted to a block-diagonal matrix; thus, Eq. (24) can be easily solved.

243

25. Based on uA , uB and unI , the displacement vector ut0 þDt is obtained. Then, the velocity vector u_ t0 þDt and the accelera€ t0 þDt are obtained. tion vector u 26. Go to step 15 for the next time step.

3.4. Algorithm summarization

4. The efficient method for two-dimensional periodic structures

Based on condensation technology, the Woodbury formula and group theory, an efficient algorithm for computing the dynamic responses of one-dimensional periodic structures is proposed. Using condensation technology, the scale of the system of linear equations is reduced. Based on the Woodbury formula, the solution of the system of linear equations can be obtained by solving another system of linear equations with a special coefficient matrix, which can be efficiently by solved using group theory. Based on these three technologies, the algorithm that is proposed in this paper is highly efficient and uses significantly much less internal storage. The procedure for solving the dynamic responses of a onedimensional periodic structure is described in the following description.

Based on similar ideas for solving the dynamic responses of one-dimensional periodic structures, i.e., condensation technology, the Woodbury formula and group theory, an efficient algorithm for two-dimensional periodic structures is presented. This section addresses a finite two-dimensional periodic structure that consists of ðM þ 1Þ  ðN þ 1Þ unit cells; the arbitrary boundary conditions are applied on the four edges of the twodimensional periodic structure, as shown in Fig. 3. For convenience, each unit cell is numbered, as shown in Fig. 3; for example, the unit cell in the southwest corner is denoted ð1; 1Þ. For the ðm; nÞ-th unit cell as shown in Fig. 4, assume that um;n b , m;n m;n , u and u denote the displacement vectors that correspond um;n d f h

1. Given the number of unit cells N þ 1. 2. The initial displacement vector, velocity vector and acceler€ t0 , respectively. ation vector are ut0 , u_ t0 and u 3. The time step is Dt and the parameters for Newmark method are given. 4. The stiffness matrix Ke , mass matrix Me and damping matrix Ce for the unit cell are formed, and the external force vectors n f (n ¼ 1; 2; . . . ; N þ 1) in Eq. (5) for all unit cells are computed.  e can be obtained from Eq. 5. The equivalent stiffness matrix K (3). 6. Based on the condensation formula, the blocks of the equivalent stiffness matrix WLL , WLR , WRL and WRR of the unit cell can be obtained. 7. Based on the information of DOFs in Eq. (10), the matrices RAA , RAB and RBA in Eq. (11) can be obtained. 8. Matrices U and V are formed, and the transformation matrix TN is computed.  can be given by Eqs. (30)–(32), 9. The diagonal-block matrix R  1 can be obtained. and the matrix R  1 ½ðTN  Ip ÞUg. 10. Compute Z ¼ ðTT  Ip ÞfR

to the DOFs on the four corners, the number of DOFs on the four m;n m;n corners is pp ; um;n and um;n denote the displacement veca , uc , ue g tors that correspond to the DOFs on the four boundaries. The number of DOFs that correspond to um;n and um;n is px , and the number a e

N

11. Compute H ¼ I þ VZ.  1 ½ðTN  Ip ÞRBA g. 12. Compute Y ¼ ðTTN  Ip ÞfR 1 13. Compute R1 BB R BA ¼ Y  Z½H ðVYÞ.

14. Compute RAA  RAB R1 BB R BA .

15. The equivalent external forces f n for all unit cells are obtained by Eq. (3). 16. The equivalent force vectors rnL and rnR for all unit cells are obtained. 17. Based on the DOF information in Eq. (10), the force vectors rA and rB can be obtained.  1 ½ðTN  Ip ÞrB g. 18. Compute y ¼ ðTT  Ip ÞfR 1

Fig. 3. Two-dimensional periodic structure.

N

1 19. Compute R1 BB rB ¼ y 1  Z½H ðVy 1 Þ.

20. Compute rA  RAB R1 BB rB . 1

1 21. Compute uA ¼ ðRAA  RAB R1 BB R BA Þ ðrA  R AB R BB rB Þ. T 1  22. Compute y ¼ ðT  Ip ÞfR ½ðTN  Ip ÞðrB  RBA uA Þg. 2

N

23. Compute uB ¼ y2  Z½H1 ðVy2 Þ. 24. The displacement vectors unI that correspond to the internal DOFs can be obtained.

Fig. 4. Unit cell for two-dimensional periodic structure.

244

X.Q. Liang et al. / Computers and Structures 182 (2017) 238–251

of DOFs that correspond to um;n and um;n is py ; um;n denotes the disc g I m;n

m;n

placement vector that corresponds to the internal DOFs; f b , f d , m;n m;n f f and f h denote the external force vectors that correspond to m;n

m;n

m;n

m;n

the DOFs on the four corners; f a , f c , f e and f g denote the external force vectors that correspond to the DOFs on the four m;n boundaries; f I denotes the external force vector that corresponds m;n m;n and pm;n denote the internal to the internal DOFs; pm;n b , pd , pf h force vectors that correspond to the DOFs on the four corners; m;n m;n and pm;n and pm;n denote the internal force vectors that a , pc , pe g correspond to the DOFs on the four boundaries. For the convenience of subsequent derivation, let p ¼ pp þ px þ py .Based on these notations, the displacement vector, external force vector and internal force vector, which corresponding to the ðm; nÞ-th unit cell, can be written in block form as h iT T T m;n T m;n T m;n T m;n T m;n T m;n T m;n T um;n ¼ ðum;n Þ ðum;n Þ a Þ ðub Þ ðuc Þ ðud Þ ðue Þ ðuf g Þ ðuh Þ ðuI ð40Þ

f

m;n

h iT T m;n T m;n T m;n T m;n T m;n T m;n T m;n T m;n T ¼ ðf m;n a Þ ðf b Þ ðf c Þ ðf d Þ ðf e Þ ðf f Þ ðf g Þ ðf h Þ ðf I Þ ð41Þ

h iT T m;n T m;n T m;n T m;n T m;n T m;n T m;n T pm;n ¼ ðpm;n a Þ ðpb Þ ðpc Þ ðpd Þ ðpe Þ ðpf Þ ðpg Þ ðph Þ 0 ð42Þ

And it is clear that

um;n ¼ um1;n a e um;n ¼ um1;n ¼ um;n1 ¼ um1;n1 b d h f

ð43Þ

um;n ¼ um;n1 c g

2

Waa Wab Wac Wad Wae Waf Wag Wah

6 6 Wba 6 6 6W 6 ca 6 6W 6 da 6 6 6 Wea 6 6 6 Wfa 6 6 6 Wga 4 Wha

uam;n

3

76 m;n 7 7 6 Wbb Wbc Wbd Wbe Wbf Wbg Wbh 7 76 ub 7 76 7 m;n 7 7 6 Wcb Wcc Wcd Wce Wcf Wcg Wch 76 uc 7 76 7 6 m;n 7 Wdb Wdc Wdd Wde Wdf Wdg Wdh 7 76 ud 7 76 7 76 7 Web Wec Wed Wee Wef Weg Weh 76 uem;n 7 76 7 76 m;n 7 Wfb Wfc Wfd Wfe Wff Wfg Wfh 76 uf 7 76 7 76 m;n 7 6 ug 7 Wgb Wgc Wgd Wge Wgf Wgg Wgh 7 54 5 Whb Whc Whd Whe Whf Whg Whh uhm;n t

2

rm;n a

3 2

pm;n a

3

6 m;n 7 6 m;n 7 6r 7 6p 7 6 b 7 6 b 7 7 6 7 6 6 rm;n 7 6 pm;n 7 6 c 7 6 c 7 7 6 7 6 6 rm;n 7 6 pm;n 7 6 d 7 6 d 7 7 6 6 ¼ 6 m;n 7 þ 6 m;n 7 7 6 re 7 6 pe 7 7 6 7 6 6 m;n 7 6 m;n 7 6 rf 7 6 pf 7 7 6 7 6 6 m;n 7 6 m;n 7 6 rg 7 6 pg 7 5 4 5 4 0 þDt

rm;n h

pm;n h ð45Þ

where the coefficient matrix denotes the equivalent stiffness matrix (b ¼ a; b; c; d; e; f ; g; h) denote the of condensation unit cell, rm;n b equivalent external force vectors of condensation unit cell. Because all unit cells of the entire structure are identical, the computational effort of the condensation is very small. Then, the system of linear equations that correspond to the entire structure can be obtained by combining the condensation equations (45) of all unit cells. Let uB denotes the displacement vector that corresponds to the DOFs enclosed by the dashed lines shown in Fig. 5, and uA denotes the displacement vector that corresponds to the remaining DOFs. The displacement vector uB can be written in the following form as

3 u1j 6 27 6 uj 7 6 7 uB ¼ 6 . 7 ; 6 .. 7 4 5 uNj 2

2

u1;j j 6 2;j 6 uj 6 ujj ¼ 6 . 6 .. 4

3

7 7 7 7; 7 5

2

umþ1;nþ1 a

3

6 mþ1;nþ1 7 um;n 5 j ¼ 4 ub

M;j

uj

m ¼ 1; 2; . . . ; M;

where m ¼ 2; 3; . . . ; M þ 1 and n ¼ 2; 3; . . . ; N þ 1. According to the form of the displacement vectors in Eq. (40), the stiffness matrix Ke , the mass matrix Me and the damping matrix Ce of the unit cell can also be written in block form. The stiffness, mass and damping matrices for the entire structure can be obtained by the superposition of the stiffness, mass and damping matrices for all unit cells, and the force and displacement vectors for the entire structure can be given by Eqs. (40) and (41), respectively. For the two-dimensional periodic structure with many unit cells, the scale of Eq. (2) is very large and timeconsuming. Therefore, similar methods for one-dimensional periodic structures, i.e., condensation technology, the Woodbury formula and group theory, are used to solve the dynamic responses for two-dimensional periodic structures.

32

umþ1;nþ1 c

ð46Þ

n; j ¼ 1; 2; . . . ; N

Then, the system of linear equations of the entire structure can be given by



RAA

RAB

RBA

RBB



uA uB





¼

rA rB



ð47Þ

4.1. The application of condensation technology To improve the computational efficiency for solving Eq. (2), which corresponds to two-dimensional periodic structures, the similar condensation technology given in Section 3 is used to reduce the scale of Eq. (2). According to Eq. (3), we can obtain the equivalent stiffness  e and equivalent force vector f m;n of the ðm; nÞ-th unit cell. matrix K Then, the system of linear equations for the ðm; nÞ-th unit cell can be expressed as

 e um;n ¼ f m;n þ pm;n K t 0 þDt

ð44Þ

where um;n t 0 þDt denotes the displacement vector of the ðm; nÞ-th unit cell at t 0 þ Dt, and pm;n is given in Eq. (42). According to Eqs. (40)–(42), Eq. (44) can be written in the block form, then using condensation technology, the internal DOFs that corresponds the displacement vector unI can be eliminated, which produces

Fig. 5. DOFs of a two-dimensional periodic structure, where the DOFs in the dashed line box corresponds to the displacement vector ub .

245

X.Q. Liang et al. / Computers and Structures 182 (2017) 238–251

where 2

4.2. The transformation based on the Woodbury formula 3

2

r1j 6 27 6 rj 7 6 7 rB ¼ 6 . 7 ; 6 .. 7 4 5 rNj

r1;j j 6 2;j 6 rj 6 rjj ¼ 6 . 6 . 4 .

3

2

7 7 7 7; 7 5

rm;n j

rM;j j

m ¼ 1; 2; . . . ; M;

3

þ rm;nþ1 rmþ1;nþ1 a e

7 6 mþ1;nþ1 þ rm;nþ1 þ rmþ1;n þ rfm;n 7 ¼6 d h 5 4 rb rmþ1;nþ1 þ rmþ1;n c g

n; j ¼ 1; 2; . . . ; N ð48Þ

and RBB is a N  N block matrix in the following form

2

RBB

W0A þ W0B

6 6 0 T 6 ðWC Þ 6 6 ¼6 0 6 6 6 6 0 4 0

W0C

0

0

0

W0A þ W0B

W0C .. . .. .

0 .. .

0

T

ðW0C Þ 0 0

W0A

0

þ W0B 0 T

W0C W0A

ðWC Þ

0

3

þ

7 7 7 7 7 7 7 7 7 7 5

ð49Þ

RBB w ¼ b

W0B

where W0c (c ¼ A; B; C) is an M  M block matrix in the following form

2 6 6 6 6 6 W0c ¼ 6 6 6 6 6 4

W1c þ W2c

W3c

0

0

0

W4c

W1c þ W2c

W3c

0

0

0

W4c

..

..

0

.

0

0

..

0

0

0

.

.

1

2

3

3

Wc þ Wc

Wc

W4c

W1c þ W2c

7 7 7 7 7 7 7 7 7 7 5

ð50Þ

Waa

6 W1A ¼ 6 4 Wba

Wab Wbb

Wac

2

3

7 Wbc 7 5;

6 W2A ¼ 6 4 Wde

Wca

Wcb Wcc 2 Wae Wad 6 3 4 T 6 WA ¼ ðWA Þ ¼ 4 Wbe Wbd Wce

2

2

0

Wch

0 Waf

6 W3C ¼ 6 4 0 Waf 0 Waf

Wcg 3 0 7 07 5; 0

0

Wed

0

0 W2B 0

3

2

0 0 6 ¼6 4 0 Wff 0 0

0

0

3

7 07 5; 0

ð51Þ

7 07 5 0

2

0

Wef

0

6 ðW ÞT 6 C 6 6 0 R ¼6 0 6 6 6 4 0 where

3

7 07 5;

Wdd

0

3

0

WA þ WB

WC

0

0

0

WA þ WB

0 .. .

0

0

WC .. . .. .

WA þ WB

WC

0

0

ðWC ÞT

WA þ WB

ðWC ÞT

0

7 07 5

Wcd 3

0 0 0 7 6 7 ¼6 4 0 Whh Whg 5; 0 Wgh Wgg 2 0 0 6 3 4 T WB ¼ ðWB Þ ¼ 6 4 0 Whf 0 Wgf 2 3 0 Wah Wag 6 7 7 W1C ¼ 6 4 0 Wbh Wbg 5;

W1B

Wee

ð52Þ

where RBB is given in Eqs. (49) and (50). Based on the idea of the Woodbury formula used for onedimensional periodic structure, we construct the matrix R by making a ‘‘small” change in the matrix RBB . First, we construct the matrix R0 by making a ‘‘small” change in the matrix RBB , i.e.,

2

in which

2

Based on the idea given in Section 3.2 for a one-dimensional periodic structure, Eq. (47) for a two-dimensional periodic structure is transformed into a new system of linear equations using the Woodbury formula, which can be efficiently solved using group theory. Eq. (47) can be solved in the similar manner as that for solving Eq. (11) for one-dimensional periodic structure. In Eq. (47), uB denotes the displacement vector that corresponds to the DOFs enclosed by the dashed lines shown in Fig. 5, and uA denotes the displacement vector that corresponds to the remaining DOFs. Since the dimension of uA is much smaller than the dimension of uB , according to the similar analysis for Eqs. (14) and (15) for onedimensional periodic structure, solving the following equation efficiently is the key for solving Eq. (47)

3

6 7 7 ¼6 4 0 Wdf 0 5 0 0 0 3 2 0 Weh Weg 7 6 7 W4C ¼ 6 4 0 Wdh Wdg 5 0 0 0 W2C

Once Eq. (47) is solved, the displacement vector um;n for the internal I DOFs can be obtained. Using condensation technology, the system of linear Eq. (2) that correspond to two-dimensional periodic structures is transformed into a smaller system of linear Eq. (47) that can be solved using few computational resources.

2

0

0

0

W4c

0

0

3 7 7 7 7 7 7 7 7 5

ð53Þ

3

7 6 6 0 0 0 0 0 7 7 6 7 6 .. .. Wc ¼ W0c þ 6 . . 0 7 7 6 0 0 7 6 6 0 0 ... 0 0 7 5 4 0 W3c 0 0 0 2 1 Wc þ W2c W3c 0 0 6 4 1 2 3 6 W W þ W W 0 c c c c 6 6 6 .. .. ¼6 0 W4c . . 6 6 6 .. 6 0 0 . W1c þ W2c 4 W3c 0 0 W4c

W4c 0 0 W3c

3 7 7 7 7 7 7 7 7 7 7 5

ð54Þ

W1c þ W2c

Then, we construct the matrix R by making a ‘‘small” change in the matrix R0 , i.e.,

2

0 6 0 6 6 6 R ¼ R0 þ 6 0 6 6 4 0 2

WC WA þ WB

6 6 ðWC ÞT 6 6 ¼6 0 6 6 6 4 0 WC

0

0

0

0

0 0 .. .. 0 . . .. 0 . 0 0

0

0

ðWC ÞT 0 0 0

3 7 7 7 7 7 7 7 5

0

WC

0

0

ðWC ÞT

WA þ WB

0 .. .

0

0

WC .. . .. .

WA þ WB

WC

0

0

ðWC ÞT

WA þ WB

ðWC ÞT

0

3 7 7 7 7 7 7 7 7 5

ð55Þ

246

X.Q. Liang et al. / Computers and Structures 182 (2017) 238–251

According to Eqs. (54) and (55), the relationship between the matrix RBB and the matrix R can be given

RBB ¼ R þ UV

ð56Þ

where

2

IM

0

60 6 6 6 U¼6 0 6 6 40

Iu

0

0 0 .. .. . . .. . Iu

07 7 7 7 0 7  Ip 7 7 05

0

0

0

IM

0

0

3

0

0

0

ð57Þ

where IM and Ip are the identity matrices of M dimensions and p dimensions, respectively, and the Iu is a matrix of M  2 dimensions, which is defined as

 Iu ¼

1 0  0 0 0

T

ð58Þ

0  0 1

of the last row of blocks. Moreover, every block Wc (c ¼ A; B; C) of R can be divided into M  M small blocks, which has the following properties, namely, the second row of blocks is the permutation of the first row of blocks, and the third row of blocks is the permutation of the second row of blocks, and so on, and the first row of blocks is the permutation of the last row of blocks. The special form of the equivalent stiffness matrix R is the same as the special form of a structure that has the cyclic periodicity in two directions. In Eq. (55), the equivalent stiffness matrix R is written as a N  N blocks form, which corresponds to a cycle periodic structure. Based on group theory for the one-dimensional periodic structure in Section 3.3, the transformation matrix TN  IMp can be given and the equivalent stiffness matrix R can be transformed into the block-diagonal matrix, i.e.,

R ¼ ðTN  IMp ÞRðTTN  IMp Þ

ð62Þ

where R is a block-diagonal matrix, and

R1 ¼ WA þ WB þ WC þ ðWC ÞT

when N is an odd integer

ð63Þ

and the matrix V can be written in the following block form as

2 6 6 6 6 6 V ¼ 6 6 6 6 6 4

GA þ GB

GC

0

0

0

ðWC ÞT

VD

V A þ VB

0

VD

0 .. .

0

0

0

0

0

0

VC .. . .. .

V A þ VB

VC

0

0

0

0

VD

V A þ VB

VC

WC

0

0

0

ðGC ÞT

GA þ GB

T

3

Rj ¼ I2  ðWA þ WB Þ þ PNj  WC þ ðPNj Þ  ðWC ÞT 2; 3; . . . ; ðN þ 1Þ=2 when N is an odd integer j¼ 1; 2; 3; . . . ; N=2 when N is an even integer

7 7 7 7 7 7 7 7 7 7 5 ð59Þ

where

2

0

0

6 6 Gc ¼ 6 0 . . . 4 W3c 0 2 0 0 VD ¼ 4 4 T ðWC Þ 0

W4c

3

7 7 ; 0 7 5 0

" Vc ¼

0 W3c

0    0 W4c 0  0

#

0

3 T    0 ðW3C Þ 5  0 0

ð60Þ

lc ¼ A; B; C Based on the Woodbury formula [54] and using Eq. (56), the solution of Eq. (52) can be obtained by using Eqs. (21)–(24). According to Eqs. (21)–(24), the major computational effort is to solve the following equation

Ry ¼ b

ð61Þ

where R is given in Eq. (55). If Eq. (61) can be efficiently solved, Eq. (52) can be efficiently solved, and the dynamic responses of the two-dimensional periodic structures can be given by efficiently solving Eq. (47) Therefore, efficiently solving Eq. (61) is the key for computing the dynamic responses of the structure.

in which PN are given in Eq. (33). According to Eqs (63) and (64), the matrix R1 can be written in block form, as shown in Fig. 6, when N is an odd integer, and matrix Rj in Eq. (64) can be written in block form, as shown in Fig. 7(a). R1 can be regarded as an equivalent stiffness matrix for a structure that has cyclic periodicity in one direction. Although the Rj shown in Fig. 7(a) does not directly correspond to a cyclic structure, using the following procedure, Rj can be transformed into a matrix that corresponds to a cyclic structure. Interchanging the i-th (i ¼ 1; 2; . . . ; M) row blocks and the ð2i  1Þ-th row blocks, the ði þ MÞ-th row blocks and the ð2iÞ-th row blocks, the i-th column blocks and the ð2i  1Þ-th column blocks, and the ði þ MÞ-th column blocks and the ð2iÞ-th column blocks, Rj can be transformed into the matrix Xj as shown in Fig. 7(b), which can be regarded as an equivalent stiffness matrix for a cyclic structure. This procedure is equivalent to the adjustment of the DOFs. Thus, the relationship between Rj and Xj can be given by

Xj ¼ ðX  Ip ÞRj ðX  Ip ÞT

ð65Þ

where

     1 0 X ¼ IM  IM  0 1

ð66Þ

Because R1 and Xj can be considered as an equivalent stiffness matrix for a cyclic structure, group theory can be reapplied. Based on the group theory for one-dimensional periodic structures in Section 3.3, the matrix Xj can be converted into a block-diagonal  j using the transformation matrix TM  I2p , i.e., matrix R

4.3. The group theory for two-dimensional periodic structures In Section 4.2, we have shown that improving the computational efficiency of Eq. (61) is the key point for solving the dynamic responses of two-dimensional periodic structures. Based on the similar idea used for one-dimensional periodic structures, group theory is used to improve the computational efficiency for twodimensional periodic structures. The coefficient matrix of Eq. (61) has the following properties, namely, the matrix R can be divided into N  N blocks, and the second row of blocks is the permutation of the first row of blocks, and the third row of blocks is the permutation of the second row of blocks, and so on, and the first row of blocks is the permutation

ð64Þ

Fig. 6. Form of the matrix R1 .

X.Q. Liang et al. / Computers and Structures 182 (2017) 238–251

247

Fig. 7. Form of the matrices: (a) Rj and (b) Xj .

 j ¼ ðTM  I2p ÞXj ðTM  I2p ÞT R

ð67Þ

In a similar manner, R1 can be converted into the block-diagonal  1 , i.e., matrix R

 1 ¼ ðTM  Ip ÞR1 ðTM  Ip ÞT R

ð68Þ

Substituting Eq. (65) into Eq. (67) gives

 j ¼ ðT  M  Ip ÞRj ðT  T  Ip Þ R M

ð69Þ

where

     0  M ¼ ðTM  I2 ÞX ¼ ðTM  I2 Þ IM  1 T IM  0 1

ð70Þ

When N is an odd integer, the transformation matrix for R1 is TM  M . When N is an even inteand the transformation matrix for Rj is T  M . Therefore, let ger, the transformation matrix for Rj is T

8  0 > < TM when N is an odd integer M 0 IðN1Þ=2  T Tx ¼ > : M IN=2  T when N is an even integer

ð71Þ

then

 j Þ ¼ ðTx  Ip Þ½diagðRj ÞðTx  Ip ÞT diagðR

5. Numerical examples

 ¼ ðT  Ip ÞRðTT  Ip Þ R

ð73Þ

in which

T ¼ Tx  ðTN  IM Þ

ð74Þ





and R is block-diagonal matrix. When M and N are even integers, R consists of MN=4 blocks, i.e.,

 ¼ diagðR  kÞ R

ð75Þ "



Rk ¼ I4  ðW1A þ W2A þ W1B þ W2B Þ þ 2 þ

i 4 ðPM Þ

0

3

T

0 ðPiM Þ

T

PiM

0

0

PiM

#  ðW3A þ W3B Þ

5  ðW4 þ W4 Þ þ ðP j  I2 Þ  ðW1 þ W2 Þ A B C C N

h i h i T T T þ ðPNj  PiM Þ  W3C þ PNj  ðPiM Þ  W4C þ ðPNj Þ  I2  ðW1C þ W2C Þ h i h i T T T T T þ ðPNj Þ  PiM  ðW4C Þ þ ðPNj Þ  ðPiM Þ  ðW3C Þ ð76Þ

k¼

M ðj  1Þ þ i; 2

i ¼ 1; 2; . . . ; M=2;

(i ¼ 1; 2; . . . ; M=2) and PNj (j ¼ 1; 2; . . . ; N=2) are given in Eq. (33). If either M or N is an odd integer, a similar formula can be given. Based on transformation matrix given in Eq. (74) and the similar formula given in Eqs. (27)–(29), the equivalent stiffness matrix R can be converted to a block-diagonal matrix. Then, the solution of Eq. (61) can be easily obtained by using the same way for onedimensional periodic structures. This subsection explains the application of group theory to twodimensional periodic structures. Based on the irreducible representation of the symmetry group, the coefficient matrix in Eq. (61) is converted to a block-diagonal matrix; thus, Eq. (61) can be easily solved. Based on the periodic properties of the structure, condensation technology and group theory were used for two-dimensional periodic structures, which produce an efficient method for computing the dynamic responses of two-dimensional periodic structures. The algorithm for two-dimensional periodic structures is based on similar ideas for one-dimensional periodic structures. Thus, the numerical procedures for computing the dynamic responses of one-dimensional and two-dimensional periodic structures are similar.

ð72Þ

Combining Eq. (62) into Eq. (72) gives

where

where I4 is the identity matrix of four dimensions, and PiM

j ¼ 1; 2; . . . ; N=2

As mentioned in Section 2, efficiently solving a system of linear equations is the key issue for computing the dynamic responses. Usually, the system of linear equations is solved by a direct method that is based on the Cholesky factorization or the preconditioned conjugate gradient method (PCGM). Based on the periodic properties of the structure, the three technologies, i.e., condensation technology, the Woodbury formula and group theory, were applied in this paper to efficiently compute the system of linear equations. Therefore, the computational efficiency and accuracy of the proposed method are discussed by the following three examples. All examples were performed on a computer which has four cores, 12 GB of internal storage and a 2.26 GHz processor. The parameters of the Newmark method are a ¼ 0:6 and d ¼ 0:5. Example 1. This example concerns the dynamic response of a perfect periodic structure that consists of masses and springs, as shown in Fig. 8(a). It consists of N unit cells and is fixed at both ends. For the unit cell shown in Fig. 8(b), the mass and stiffnesses are m ¼ 2  103 kg, ka ¼ 2  105 N=m and kb ¼ 3  105 N=m. The structure has 2N  1 DOFs, the initial displacements and velocities

248

X.Q. Liang et al. / Computers and Structures 182 (2017) 238–251

Fig. 8. (a) Periodic structure of Example 1 and (b) the unit cell.

are zero for all DOFs, and the external force 105 sin½np=ð2Nþ 1Þ sinðtÞN is applied at the n-th node in the horizontal direction. The Rayleigh damping is C ¼ 0:05K. The numerical results of the direct method were employed as the reference solution to compute the relative errors of the proposed method from

e¼

ku  uD k2 kuD k2

in which e represents the relative errors for the displacement vector, u and uD are the displacement vectors computed using the proposed method and the direct method, respectively.

In this example, the two time steps 0.2 s and 2 s were employed, and the number of time steps was 105 . Figs. 9 and 10 give the relative errors for the displacement responses for the time steps 0.2 s and 2 s, respectively. In each figure, subfigure (a) and subfigure (b) correspond to the periodic structures, which consist of 49 unit cells and 149 unit cells, respectively. Figs. 9 and 10 show that the proposed method has very high precision in solving the system of linear equations and is independent of the time steps and the number of unit cells. Example 2. This example concerns the dynamic response of a onedimensional periodic structure (refer to Fig. 11), which consists of N unit cells. The unit cell in Fig. 12(a) is a square plate with a

Fig. 9. Relative errors of the proposed method when time step is 0.2 s: (a) N ¼ 49 and (b) N ¼ 149.

Fig. 10. Relative errors of the proposed method when time step is 2 s: (a) N ¼ 49 and (b) N ¼ 149.

249

X.Q. Liang et al. / Computers and Structures 182 (2017) 238–251

Fig. 11. One-dimensional periodic structure for Example 2.

Fig. 12. Unit cell: (a) geometry size and (b) mesh of FEM.

circular hole in the center. The geometric sizes of the unit cell were L ¼ 0:1 m, r ¼ 0:025 m and the thickness h ¼ 0:025 m. The physical properties were the density q ¼ 9:8  103 kg=m3 , Young’s modulus E ¼ 2:1  1011 N=m2 and Poisson’s ratio l ¼ 0:3. The unit cell was discretized using a FEM with triangular elements, as shown in Fig. 12(b), with 192 triangular elements and 120 nodes (240 DOFs). The structure is clamped at both ends of structure and at the boundaries of the circular hole of every unit cell. The initial displacements and velocities of all DOFs were zero. The external forces were applied to the nodes that correspond to the top of the structure, had the form 105 sinðlp=NLÞ sinðtÞN and were oriented in the y-direction, in which l denotes the distance from the left end. The Rayleigh damping of this example is C ¼ 0:05K þ 0:06M. To compare the computational efficiency, the system of linear equations for the Newmark method is solved by the following three methods, i.e., the direct method, the PCGM and the proposed method. The time step and the interval of integration were 0.002 s and 0–2 s. Fig. 13 shows the central processing unit (CPU) times for

the three methods, where the solid line and the solid line with solid circles represent the CPU times by the direct method and the proposed method, respectively, and the solid lines with triangles, squares and circles represent the CPU times used by the PCGM with the three different convergence tolerances s ¼ 104 ; 106 and

108 , respectively. Fig. 13 shows that the proposed method is much more efficient than the remaining two methods, and the efficiency becomes more significant with an increase in the number of unit cells. When the number of unit cells is 1499 (approximately 300 thousand DOFs), the efficiency of the proposed method is 2.3 times better than the efficiency of the direct method, and is 15, 21 and 25 times better than the efficiency of the PCGM for the three different convergence tolerances

s ¼ 104 ; 106 and 108 , respectively.

Example 3. This example concerns the dynamic responses of a two-dimensional periodic structure (refer to Fig. 14) that consists of N  N unit cells. The unit cell is identical to the unit cell in Example 2. The structure is clamped at the left end and at the

Fig. 13. CPU times for one-dimensional periodic structures.

250

X.Q. Liang et al. / Computers and Structures 182 (2017) 238–251

6. Conclusion

Fig. 14. Two-dimensional periodic structure for Example 3.

boundaries of the circular hole of every unit cell. The external forces were applied to the nodes that correspond to the right end of the structure, had the form 103 sinðtÞ N and were oriented in the xdirection. The initial displacements and velocities were zero for all DOFs. The Rayleigh damping is C ¼ 0:05K þ 0:06M. The time step and the interval of integration were 0.002 s and 0–2 s, respectively. Fig. 15 shows the CPU times of the integration procedure for the three methods, where the solid line and the solid line with solid circles give the CPU times by the direct method and the proposed method, respectively, and the solid lines with triangles, squares and circles give the CPU times used by the PCGM with the three different convergence tolerances

s ¼ 104 ; 106 and 108 , respectively. Fig. 15 shows that the proposed method is much more efficient than the remaining two methods, and the efficiency becomes more significant as the number of unit cells increases. When the number of unit cells is 127  127 (approximately 3.25 million DOFs), the efficiency of the proposed method is approximately five times better than the efficiency of the direct method, and is 12.3, 16.7 and 20.1 times better than the efficiency of the PCGM for the three convergence tolerances, respectively.

Based on the efficient numerical method for computing a system of linear equations, an efficient algorithm has been presented for solving the dynamic responses of periodic structures. For computing the system of linear equations with high efficiency, the three technologies, i.e., the condensation technology, the Woodbury formula and group theory, were applied. Firstly, based on the periodic properties of the structure and condensation technology, the system of linear equations that corresponds to the periodic structures is reduced to a small-scale system of linear equations. Then, based on the Woodbury formula, the solution of the small-scale system of linear equations is obtained by solving a new system of linear equations whose coefficient matrix corresponds to a cycle periodic structure. Finally, using group theory, the new system of linear equations is efficiently solved. Thus, an efficient algorithm for computing the dynamic responses of one-dimensional and two-dimensional periodic structures is generated. The proposed method for computing the dynamic responses of the periodic structures has the following advantages: (1) Using the periodic properties of the structure, only the mass, damping and stiffness matrices of one unit cell are needed in the procedure of the condensation. Thus, it is very convenient and has superior computational efficiency. (2) Based on the Woodbury formula and group theory, the coefficient matrix of the system of linear equations can be converted into a block diagonal matrix, which can be solved with high computational efficiency and few computational resources. (3) The numerical examples demonstrate that the proposed method is more efficient than the direct method and the preconditioned conjugate gradient method. Acknowledgements The authors are grateful for the support of the Natural Science Foundation of China (Nos. 11272076 and 11572076); the 973 program (No. 2014CB049000); Program for New Century Excellent Talents in University (No. NCET-13-0072); and the Fundamental Research Funds for the Central Universities (Nos. DUT13LK12 and DUT14YQ202).

References [1] Lebee A, Sab K. Homogenization of a space frame as a thick plate: application of the bending-gradient theory to a beam lattice. Comput Struct 2013;127:88–101. [2] Nukala PKVV, Barai P, Zapperi S, Alava MJ, Simunovic S. Fracture roughness in three-dimensional beam lattice systems. Phys Rev E 2010;82.

Fig. 15. CPU times of the integration procedure for two-dimensional periodic structures.

X.Q. Liang et al. / Computers and Structures 182 (2017) 238–251 [3] Yi SN, Cheng GD, Xu L. Stiffness design of heterogeneous periodic beam by topology optimization with integration of commercial software. Comput Struct 2016;172:71–80. [4] Domagalski L, Jedrysiak J. Geometrically nonlinear vibrations of slender mesoperiodic beams. The tolerance modeling approach. Compos Struct 2016;136:270–7. [5] Su WZ, Liu ST. Vibration analysis of periodic cellular solids based on an effective couple-stress continuum model. Int J Solids Struct 2014;51:2676–86. [6] Kushwaha MS, Halevi P, Dobrzynski L, Djafari-Rouhani B. Acoustic band structure of periodic elastic composites. Phys Rev Lett 1993;71:2022–5. [7] Michel JC, Moulinec H, Suquet P. Effective properties of composite materials with periodic microstructure: a computational approach. Comput Meth Appl Mech Eng 1999;172:109–43. [8] Dizy J, Palacios R, Pinho ST. Homogenisation of slender periodic composite structures. Int J Solids Struct 2013;50:1473–81. [9] Klos JW, Krawczyk M, Dadoenkova YS, Dadoenkova NN, Lyubchanskii IL. Photonic-magnonic crystals: multifunctional periodic structures for magnonic and photonic applications. J Appl Phys 2014;115. [10] Liu Y, Gao LT. Explicit dynamic finite element method for band-structure calculations of 2D phononic crystals. Solid State Commun 2007;144:89–93. [11] Wang YZ, Li FM, Huang WH, Jiang X, Wang YS, Kishimoto K. Wave band gaps in two-dimensional piezoelectric/piezomagnetic phononic crystals. Int J Solids Struct 2008;45:4203–10. [12] Tassilly E. Propagation of bending waves in a periodic beam. Int J Eng Sci 1987;25:85–94. [13] Gavric L. Computation of propagative waves in free rail using a finite element technique. J Sound Vib 1995;185:531–43. [14] Shen MR, Cao WW. Acoustic bandgap formation in a periodic structure with multilayer unit cells. J Phys D: Appl Phys 2000;33:1150–4. [15] Phani AS, Woodhouse J, Fleck NA. Wave propagation in two-dimensional periodic lattices. J Acoust Soc Am 2006;119:1995–2005. [16] Sonekar P, Mitra M. A wavelet-based model of one-dimensional periodic structure for wave-propagation analysis. P Roy Soc A-Math Phy 2010;466:263–81. [17] Tian BY, Tie B, Aubry D, Su XY. Elastic wave propagation in periodic cellular structures. CMES-Comp Model Eng 2011;76:217–33. [18] Trainiti G, Rimoli JJ, Ruzzene M. Wave propagation in periodically undulated beams and plates. Int J Solids Struct 2015;75–76:260–76. [19] Vonflotow AH. Disturbance propagation in structural networks. J Sound Vib 1986;106:433–50. [20] Cai GQ, Lin YK. Wave propagation and scattering in structural networks. J Eng Mech-ASCE 1991;117:1555–74. [21] Yu DL, Wen JH, Shen HJ, Xiao Y, Wen XS. Propagation of flexural wave in periodic beam on elastic foundations. Phys Lett A 2012;376:626–30. [22] Petersson LER, Jin JM. A two-dimensional time-domain finite element formulation for periodic structures. IEEE T Antenn Propag 2005;53:1480–8. [23] Schmidt K, Kauf P. Computation of the band structure of two-dimensional photonic crystals with hp finite elements. Comput Meth Appl Mech Eng 2009;198:1249–59. [24] Olhoff N, Niu B, Cheng G. Optimum design of band-gap beam structures. Int J Solids Struct 2012;49:3158–69. [25] Xiang HJ, Shi ZF. Analysis of flexural vibration band gaps in periodic beams using differential quadrature method. Comput Struct 2009;87:1559–66. [26] Zuo SL, Li FM, Zhang CZ. Numerical and experimental investigations on the vibration band-gap properties of periodic rigid frame structures. Acta Mech 2016;227:1653–69. [27] Wu ZJ, Li FM, Wang YZ. Vibration band gap properties of periodic Mindlin plate structure using the spectral element method. Meccanica 2014;49:725–37. [28] Engels RC. Response of infinite periodic structures. J Sound Vib 1980;69:181–97. [29] Wu ZJ, Li FM, Wang YZ. Study on vibration characteristics in periodic plate structures using the spectral element method. Acta Mech 2013;224:1089–101.

251

[30] Cai CW, Cheung YK, Chan HC. Uncoupling of dynamic equations for periodic structures. J Sound Vib 1990;193:253–63. [31] Zhou CW, Laine JP, Ichchou MN, Zine AM. Wave finite element method based on reduced model for one-dimensional periodic structures. Int J Appl Mech 2015;7:1550018. [32] Zhou CW, Laine JP, Ichchou MN, Zine AM. Multi-scale modelling for twodimensional periodic structures using a combined mode/wave based approach. Comput Struct 2015;154:145–62. [33] Hawreliak JA, Lind J, Maddox B, Barham M, Messner M, Barton N, et al. Dynamic behavior of engineered lattice materials. Sci Rep 2016;6. [34] Luongo A, Romeo F. Real wave vectors for dynamic analysis of periodic structures. J Sound Vib 2005;279:309–25. [35] Duhamel D, Mace BR, Brennan MJ. Finite element analysis of the vibrations of waveguides and periodic structures. J Sound Vib 2006;294:205–20. [36] Mencik JM. On the low- and mid-frequency forced response of elastic structures using wave finite elements with one-dimensional propagation. Comput Struct 2010;88:674–89. [37] Gao Q, Yao WA, Wu F, Zhang HW, Lin JH, Zhong WX, et al. An efficient algorithm for computing the dynamic responses of one-dimensional periodic structures and periodic structures with defects. Comput Mech 2013;52:525–34. [38] Gao Q, Zhang HW, Zhong WX, Howson WP, Williams FW. An accurate and efficient method for dynamic analysis of two-dimensional periodic structures. Int J Appl Mech 2016;08:1650013. [39] Gao Q, Wu F, Zhang HW, Zhong WX, Howson WP, Williams FW. Exact solutions for dynamic response of a periodic spring and mass structure. J Sound Vib 2012;331:1183–90. [40] Wu F, Gao Q, Zhong WX. Subdomain precise integration method for periodic structures. Shock Vib 2014;2014:1–11. [41] Zhong WX, Qiu CH. Analysis of symmetric or partially symmetric structures. Comput Meth Appl Mech Eng 1983;38:1–18. [42] Zingoni A. An efficient computational scheme for the vibration analysis of high tension cable nets. J Sound Vib 1996;189:55–79. [43] Kangwai RD, Guest SD, Pellegrino S. An introduction to the analysis of symmetric structures. Comput Struct 1999;71:671–88. [44] Kangwai RD, Guest SD. Symmetry-adapted equilibrium matrices. Int J Solids Struct 2000;37:1525–48. [45] Zingoni A. On group-theoretic computation of natural frequencies for springmass dynamic systems with rectilinear motion. Commun Numer Meth En 2008;24:973–87. [46] Mohan SJ, Pratap R. A group theoretic approach to the linear free vibration analysis of shells with dihedral symmetry. J Sound Vib 2002;252:317–41. [47] Kaveh A, Nikbakht M. Improved group-theoretical method for eigenvalue problems of special symmetric structures, using graph theory. Adv Eng Softw 2010;41:22–31. [48] Zingoni A. Group-theoretic exploitations of symmetry in computational solid and structural mechanics. Int J Numer Meth Eng 2009;79:253–89. [49] Newmark NM. A method of computation for structural dynamics. ASCE J Eng Mech 1959;85:67–94. [50] Chung J, Hulbert GM. A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-a method. J Appl MechTrans ASME 1993;60:371–5. [51] Bathe KJ. Conserving energy and momentum in nonlinear dynamics: A simple implicit time integration scheme. Comput Struct 2007;85:437–45. [52] Bathe KJ, Noh G. Insight into an implicit time integration scheme for structural dynamics. Comput Struct 2012;98–99:1–6. [53] Park KC, Underwood PG. A variable-step central difference method for structural dynamics analysis, 1: theoretical aspects. Comput Meth Appl Mech Eng 1980;22:241–58. [54] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes in C: the art of scientific computing: The Press Syndicate of the University of Cambridge; 1988–1992. [55] Jacob HG. A theorem on Kronecker products. Bull Am Math Soc 1953;59:235.