Sensitivity analysis for free vibration of rectangular plate

Sensitivity analysis for free vibration of rectangular plate

Journal of Sound and Vibration 332 (2013) 1610–1625 Contents lists available at SciVerse ScienceDirect Journal of Sound and Vibration journal homepa...

3MB Sizes 0 Downloads 64 Views

Journal of Sound and Vibration 332 (2013) 1610–1625

Contents lists available at SciVerse ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Sensitivity analysis for free vibration of rectangular plate Myung-Soo Choi a,n, Jung-Hwan Byun b a b

Department of Maritime Police Science, Chonnam National University, San 96-1, Dundeok-dong, Yeosu, Jeonnam 550-749, South Korea Faculty of Marine Technology, Chonnam National University, San 96-1, Dundeok-dong, Yeosu, Jeonnam 550-749, South Korea

a r t i c l e in f o

abstract

Article history: Received 20 September 2011 Received in revised form 14 September 2012 Accepted 17 September 2012 Handling Editor: S. Ilanko Available online 1 December 2012

This paper presents an effective sensitivity analysis algorithm for free vibration of a rectangular plate structure by using the finite element-transfer stiffness coefficient method (FE-TSCM). The basic concept of FE-TSCM combines both the modeling technique of the finite element method (FEM) and the transfer technique of the transfer stiffness coefficient method (TSCM) in order to benefit from the merits of both FEM and TSCM in the dynamic analysis of a structure. From the results computing the sensitivities of eigenvalues and eigenvectors for a simply supported rectangular plate structure by FE-TSCM and Fox’s method, we can confirm that FE-TSCM has merits from the viewpoint of computational time and memory in sensitivity analysis for free vibration of the rectangular plate structure with a large number of degrees-of-freedom. In addition, when the design variables are the thicknesses of the plate elements, the results computing the sensitivities of the first eigenvalues for the rectangular plate structures with clamped edges are illustrated in detail by using the computation program introducing the concept of FE-TSCM. & 2012 Elsevier Ltd. All rights reserved.

1. Introduction When designing mechanical or structural systems, it is important to understand the characteristics of the systems. The sensitivity analysis provides designers with reasonable and effective data for structural modification, optimization, and so on [1–2]. Since the 1960s, many studies for obtaining the dynamic sensitivity of systems have been carried out [3–5]. Most of the studies have focused on obtaining the derivatives of eigenvalues and eigenvectors with respect to design variables by using large matrices with total degrees-of-freedom of the discrete system; for example, the system stiffness and mass matrices of the finite element method (FEM) [6]. However, large matrices in dynamic analysis for a structure with a lot of degrees-of-freedom may result in many problems from the viewpoint of computational time and memory [7]. In order to overcome these problems, many researchers suggested a variety of methods [8–11] using the transfer matrix method [12]. As a result of studying computational algorithms to analyze effectively the free vibration of plate structures, the finite element-transfer stiffness coefficient method (FE-TSCM) was suggested by Choi [13]. The basic concept of FE-TSCM combines both the modeling technique of FEM and the transfer technique of the transfer stiffness coefficient method (TSCM) [14,15] in order to benefit from the merits of both FEM and TSCM in the free vibration analysis of the rectangular plate structure.

n

Corresponding author. Tel.: þ 82 61 659 7183; fax: þ 82 61 659 7189. E-mail address: [email protected] (M.-S. Choi).

0022-460X/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jsv.2012.09.045

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

1611

The sensitivity analysis algorithm for free vibration of a straight-line beam structure by using the transfer stiffness coefficient method (TSCM) was suggested in a previous paper [16]. However, TSCM in the previous study was limited to the sensitivity analysis of a simple one-dimensional structure such as beam and shaft. Therefore, it was difficult to apply the previous TSCM to the sensitivity analysis of a two-dimensional structure such as plate, panel and shell. In addition, the previous TSCM has drawback, which requires much computational time for sensitivity analysis. In this paper, the authors apply FE-TSCM to the free vibration sensitivity analysis algorithm for a rectangular plate structure. In order to confirm the effectiveness of FE-TSCM, a simply supported rectangular plate structure is given as a computation model for obtaining sensitivities of eigenvalues and eigenvectors. The computational results of FE-TSCM are compared with those of Fox’s method [3] in terms of computational accuracy, time, and memory. In addition, when the design variables are the thicknesses of the plate elements, the results computing the sensitivities of the first eigenvalues for the rectangular plate structures with clamped edges are illustrated in detail by using the computation program introducing the concept of FE-TSCM. 2. Algorithm 2.1. Modeling of rectangular plate structure The rectangular plate structure of Fig. 1 consists of a rectangular plate and elastic springs supporting the plate from the base. The rectangular plate is divided into m strips, and each strip is subdivided into n  1 rectangular plate elements, as shown in Fig. 2. The sections between the strips are called nodal lines, which are designated as the nodal line 1, the nodal line 2, y, and the nodal line m þ1 consecutively from the left-hand edge of the plate to the right-hand edge.

Fig. 1. Rectangular plate structure.

Fig. 2. Nodal lines, plate elements and nodes of strip i and the definition of positive direction of state variables.

1612

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

Each nodal line has n nodes, and a node has 3 degrees-of-freedom for analyzing the bending vibration of the plate T

structure. The displacement vector of the node jðuj ¼ fw, f, ygj Þ is composed of a z-directional displacement (w), an xdirectional angular displacement (f) and a y-directional angular displacement (y). The force vector of the node jðf j ¼ fq,M,NgTj Þ consists of a z-directional force (q), an x-directional moment (M) and a y-directional moment (N). T

Therefore, each nodal line has 3n degrees-of-freedom. The displacement vector of the nodal line iðUi ¼ fuT1 ,uT2 ,. . .,uTn gi Þ is T

T

T T

comprised of the displacement vectors of all nodes in the nodal line i. The force vector of the nodal line iðFi ¼ ff 1 ,f 2 ,. . .,f n gi Þ consists of the force vectors of all nodes in the nodal line i. The strip i represents the ith strip of the plate, which exists between the nodal line i and the nodal line iþ1. Each strip has 2n nodes. The relationship between the force and the displacement vectors at both sections of the strip i can be written as 8 L9 2 38 L 9 ~ = < F~ = A~ i B~ i < U i i 4 5 ¼ , (1) T R : F~ ; ~R; C~ i : U B~ i

i

i

~ and displacement vector ðUÞ ~ represent the left- and rightwhere the subscript i and superscripts L and R of the force vector ðFÞ ~ ~ hand sides of the strip i. The submatrices of the dynamic stiffness matrix of the strip i, A i , B i , and C~ i , are as follows: 2 3 A~ i B~ i ~ i o2 M ~ i ¼K ~ i, ~ i lM 4 T 5¼K (2) C~ i B~ i

~ i and K ~ i are where o is the natural circular frequency of free vibration and the eigenvalue l is the square of o. The matrices M the mass and stiffness matrices for the strip i, which can be obtained respectively by assembling the mass and stiffness matrices ~ i and M ~ i can be also obtained by assembling the triangular of the rectangular plate elements in the strip i [13]. The matrices K plate elements. If there are parts that support some nodes of the plate from the base, they are modeled as three elastic springs per node. The springs consist of one linear and two rotational springs of which constants are kzj, Kxj and Kyj in the case of the node j. The point stiffness matrix Pi is a diagonal matrix and consists of linear and rotational spring constants of the nodal line i as follows:   (3) Pi ¼ diag kz1 ,K x1 ,K y1 ,kz2 ,K x2 ,K y2 ,. . .,kzn ,K xn ,K yn i , where ‘diag’ denotes the diagonal matrix. In the present method, the boundary conditions of structures are modeled as the elastic springs. For example, in the case of the free condition at the left-hand edge of the plate, we consider the point stiffness matrix as P1 ¼ diagð0,0,0,. . .,0Þi , and in the case of the clamped condition as P1 ¼ diagð1,1,1,. . .,1Þi , where N is considered as 1.0E20 in the numerical calculation. 2.2. Free vibration analysis algorithm of rectangular plate structure Because the free vibration analysis precedes the sensitivity analysis, the algorithm for analyzing the free vibration of the rectangular plate structure using FE-TSCM [13] is described briefly. In the case of being some elastic springs supporting a nodal line, the force vector of the left-hand side of the nodal line is not equal to the force vector of the right-hand side of the nodal line. Therefore, a nodal line is analytically divided into the left- and right-hand sides of the nodal line to explain easily the present method. We denote the left-hand side of a nodal line with the head mark, , on the symbols. The right-hand side of the nodal line is denoted with the head mark, ^, on the same symbols. In the case of being some elastic springs supporting a nodal line, the displacement vector of the left-hand side of the nodal line is equal to the displacement vector of the right-hand side of the nodal line. The relationships between the force and the displacement vectors at the left- and right-hand sides of the nodal line i are defined by using the stiffness coefficient matrices (size: 3n  3n) as follows: T

Fi ¼ Si Ui ,Si ¼ Si ,

(4)

T F^ i ¼ S^ i Ui , S^ i ¼ S^ i ,

(5)

where T

Fi ¼ fq1 ,M 1 ,N 1 ,q2 ,M2 ,N 2 ,. . .,qn ,Mn ,N n gi , ^ 1,N ^ 1 , q^ , M ^ 2,N ^ 2 ,. . ., q^ , M ^ n,N ^ n gT , F^ i ¼ fq^ 1 , M 2 n i T

Ui ¼ fw1 , f1 , y1 ,w2 , f2 , y2 ,. . .,wn , fn , yn gi :

(6)

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

1613

The matrices Si and S^ i are the symmetric matrices (see Appendix A). The positive directions of the state variables are defined by the arrows in Fig. 2. If there are springs at the nodal line i, the force balance at the nodal line i can be expressed as F^ i ¼ Fi þ Pi Ui ,

(7)

where the point stiffness matrix Pi is defined in Eq. (3). Consider the transfer of the stiffness coefficient from the left-hand side of the nodal line i to the right-hand side of the nodal line i, across the spring elements of the nodal line i. We can derive the matrix S^ i from Eqs. (4), (5) and (7) as follows: S^ i ¼ Si þ Pi

ði ¼ 1,2,3,. . .,m þ 1Þ:

(8)

Eq. (8) is called the point transmission rule of the stiffness coefficient matrix in the present method. The left-hand section of the strip i is equal to the right-hand side of the nodal line i, and the right-hand section of the strip i is equal to the left-hand side of the nodal line iþ 1. Therefore, Eq. (1) can be changed into 3( ( ) 2 ~ ) A i B~ i Ui F^ i 5 ¼4 T , (9) Ui þ 1 Fi þ 1 C~ i B~ i

where L R ~ L ,Ui þ 1 ¼ U ~ R: F^ i ¼ F~ i ,Fi þ 1 ¼ F~ i ,Ui ¼ U i i

(10)

Consider the transfer of the stiffness coefficient from the right-hand side of the nodal line i to the left-hand side of the nodal line i þ1, across the strip i. We can derive the matrix Si þ 1 from the equation substituted with i¼ iþ1 in Eq. (4) and from Eqs. (5) and (9) as follows: T Si þ 1 ¼ C~ i þ B~ i Vi

ði ¼ 1,2,3,. . .,mÞ,

(11)

where ~ Gi ¼ GTi ¼ S^ i þ A~ i ,Vi ¼ G1 i Bi:

(12)

Eq. (11) is called the field transmission rule of the stiffness coefficient matrix. In the present method, the boundary condition of the structure is modeled by springs at the nodal lines. Therefore, the force vectors at the left-hand side of the nodal line 1 and at the right-hand side of the nodal line mþ1 can be considered as null vectors as follows: F1 ¼ 0,

(13)

F^ m þ 1 ¼ 0:

(14)

From Eq. (13) and the equations substituted with i¼1 in Eqs. (5) and (7), the matrix S^ 1 can be derived as S^ 1 ¼ P1 :

(15)

After first obtaining the matrix S^ 1 from Eq. (15), the matrix S^ m þ 1 , the stiffness coefficient matrix at the right-hand side of the last nodal line, can be finally obtained by applying Eqs. (8) and (11) consecutively. From Eq. (14) and the equation substituted with i¼mþ1 in Eq. (5), we obtain S^ m þ 1 Um þ 1 ¼ 0:

(16)

Therefore, the frequency equation can be derived from Eq. (16) as follows: det S^ m þ 1 ðlÞ ¼ 0:

(17)

The natural mode corresponding to each natural frequency is computed by the following procedure. First, the displacement vector (Um þ 1) at the last nodal line is computed from Eq. (16). The displacement vectors of the other nodal lines can then be successively computed by using Eq. (18) obtained from Eqs. (5), (9) and (12). Ui ¼ V i Ui þ 1

ði ¼ m,m1, m2,. . .,1Þ:

(18)

2.3. Free vibration sensitivity analysis algorithm of rectangular plate structure 2.3.1. Definition of sensitivity stiffness coefficient matrix From differentiating Eqs. (4) and (5) with respect to a design variable d, we get n

n

Fi ¼ Si Ui þ Si Uni ,

(19)

n n F^ i ¼ S^ i Ui þ S^ i Uni ,

(20)

1614

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

n n where the symbol * denotes @/@d. The matrices Si and S^ i are the derivatives of the stiffness coefficient matrices with respect to the design variable d; they are defined as the sensitivity stiffness coefficient matrices (size: 3n  3n) at the leftn n and right-hand sides of the nodal line i. The matrices Si and S^ i are the symmetric matrices, because the matrices Si and S^ i are the symmetric matrices.

2.3.2. Transmission rule of sensitivity stiffness coefficient matrix Differentiation of Eq. (7) with respect to the design variable d gives n

Fni ¼ Fi þ Pni Ui þ Pi Uni :

(21)

n

If we know the sensitivity stiffness coefficient matrix Si at the left-hand side of the nodal line i, we can obtain the n sensitivity stiffness coefficient matrix S^ i at the right-hand side of the nodal line i from Eqs. (8), (19)–(21) as follows: n n S^ i ¼ Si þ Pni

ði ¼ 1,2,3,. . .,m þ1Þ:

(22)

Eq. (22) is called the point transmission rule of the sensitivity stiffness coefficient matrix in the present method. From differentiating Eq. (9) with respect to the design variable d, we obtain 8 n 9 2 3 3( ) 2 ~ ) n n ( < F^ = A i B~ i A~ i B~ i Ui Uni i 4 5 4 5 ¼ þ : (23) T n Tn Ui þ 1 Uni þ 1 : Fni þ 1 ; C~ i B~ i C~ i B~ i n n If we take the matrix S^ i at the right-hand side of the nodal line i, we can obtain the matrix Si þ 1 at the left-hand side of the nodal line iþ1 from the equation substituted with i ¼iþ1 in Eq. (19) and from Eqs. (12), (20) and (23). n

n

Tn

T

Si þ 1 ¼ C~ i þ B~ i Vi þ B~ i Vni

ði ¼ 1,2,3,. . .,mÞ,

(24)

where  n  n n B~ i þ Gni Vi : Gni ¼ S^ i þ A~ i ,Vni ¼ G1 i

(25)

Eq. (24) is called the field transmission rule of the sensitivity stiffness coefficient matrix. 2.3.3. Derivative of eigenvalue From differentiating Eqs. (13) and (14) with respect to a design variable d, we obtain n

F1 ¼ 0,

(26)

n F^ m þ 1 ¼ 0:

(27)

The matrix Sn1 can be derived from Eqs. (15) and (26) and the equations substituted with i ¼1 in Eqs. (20) and (21) as follows: n S^ 1 ¼ Pn1 :

(28)

n n After first obtaining the matrix S^ 1 from Eq. (28), the matrix S^ m þ 1 , the sensitivity stiffness coefficient matrix at the righthand side of the last nodal line, can be finally obtained by applying Eqs. (22) and (24) consecutively. From Eq. (27) and the equation substituted with i¼mþ1 in Eq. (20), we obtain n S^ m þ 1 Um þ 1 þ S^ m þ 1 Unm þ 1 ¼ 0:

(29)

Premultiplication of Eq. (29) by UTm þ 1 yields n UTm þ 1 S^ m þ 1 Um þ 1 þ UTm þ 1 S^ m þ 1 Unm þ 1 ¼ 0:

UTm þ 1 S^ m þ 1 Unm þ 1

UTm þ 1 S^ m þ 1

(30)

T

in Eq. (30) becomes zero because ¼ 0 from the transposition of Eq. (16) and the symmetry of the stiffness coefficient matrix. Therefore, we can reduce Eq. (30) shortly as follows:  n  n n f l ¼ UTm þ 1 S^ m þ 1 l Um þ 1 ¼ 0: (31) Because the function f(l*) of Eq. (31) is the linear function of unknown eigenvalue derivative l* in the numerical calculation finding l* (see Appendix B), we can obtain the eigenvalue derivative l* from Eq. (32) after calculating f(a) and f(b) of two trial values (a and b) for the eigenvalue derivative l* in Eq. (31),

ln ¼

af ðbÞbf ðaÞ : f ðbÞf ðaÞ

(32)

If one chooses the two trial values (a and b) as very close values, for example, a ¼0 and b ¼0.1, computational result may be bad. Therefore, two trial values are selected as a ¼0 and b ¼1.0e7 in computational program.

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

1615

2.3.4. Derivative of eigenvector After determining the derivative of the eigenvalue, the derivative of the eigenvector is obtained. Eq. (29) is expressed by S^ m þ 1 Unm þ 1 ¼ gm þ 1 ,

(33)

n gm þ 1 ¼ S^ m þ 1 Um þ 1 :

(34)

where a known vector gm þ 1 is

Because the matrix S^ m þ 1 is the singular matrix, we cannot find out the vector Unm þ 1 directly from Eq. (33). Because the displacement vectors of the free vibration have already been determined at the stage of the sensitivity analysis, the value of Eq. (35) becomes a constant, UTm þ 1 Um þ 1 ¼ g:

(35)

From differentiating Eq. (35) with respect to the design variable d, we get UTm þ 1 Unm þ 1 ¼ 0:

(36)

If we eliminate an equation from the simultaneous equations, Eq. (33), the remainder simultaneous equations become 0 S^ m þ 1 Unm þ 1 ¼ g0m þ 1 :

Eqs. (36) and (37) can be written together as 2 4

0 S^ m þ 1

UTm þ 1

3 5Un mþ1 ¼



g0m þ 1 0

(37)

 :

Therefore, the derivative of displacement vector at the nodal line m þ1 can be obtained as 2 0 31  0  S^ m þ 1 gm þ 1 n 4 5 : Um þ 1 ¼ 0 UTm þ 1

(38)

(39)

The relationship of the derivatives of displacement vectors at both sides of the strip i is obtained from Eqs. (12), (20), (23) and (25) as follows: Uni ¼ Vni Ui þ 1 þVi Uni þ 1

ði ¼ m,m1, m2,. . .,1Þ:

(40)

After obtaining the derivative of the displacement vector of the nodal line m þ1 from Eq. (39), the derivatives of the displacement vectors at the other nodal lines can be obtained by applying Eq. (40) successively from the last nodal line to the first nodal line of the structure.

2.3.5. Computation procedure of sensitivity analysis The following computation procedure is performed for obtaining the derivatives of eigenvalues and eigenvectors of the rectangular plate structure shown in Fig. 1 in the computer program using FE-TSCM. (1) We find out the natural frequency (o) of the rectangular plate structure from the free vibration analysis using FE-TSCM [13], as described in Section 2.2. The eigenvalue (l) is calculated from the square of o. (2) After we compute the displacement vectors ðU1 ,U2 ,U3 ,. . .,Um þ 1 Þ of all the nodal lines for the rectangular plate structure from the free vibration analysis using FE-TSCM, the displacement vectors are normalized by considering the largest value among the components in all the displacement vectors as one. The eigenvector ðX ¼ T

(3) (4) (5) (6) (7)

fUT1 ,UT2 ,UT3 ,. . .,UTm þ 1 g Þ consists of the displacement vectors of all the nodal lines. n n After we compute S^ 1 ðaÞ from Eq. (28), S^ m þ 1 ðaÞ is computed from the transmission rules, Eqs. (22) and (24), where a is the first trial value for the eigenvalue derivative (l*). n S^ m þ 1 ðbÞ is computed from the same manner as the above procedure (3), where b is the second trial value for the eigenvalue derivative. After computing f(a) and f(b), respectively, from Eq. (31), the eigenvalue derivative (l*) is calculated from Eq. (32). n n S^ m þ 1 ðl Þ is computed from the same manner as the above procedure (3) with the exact eigenvalue derivative (l*) obtained in the above procedure (5). After the derivative of the displacement vector of the last nodal line, Unm þ 1 , is computed from Eq. (39), the derivatives   of the displacement vectors of the other nodal lines Un1 ,Un2 ,Un3 ,. . .,Unm are successively computed from Eq. (40). T

nT The eigenvector derivative ðXn ¼ fU1nT ,U2nT ,U3nT ,. . .,Um þ 1 g Þ consists of the derivatives of the displacement vectors of all the nodal lines.

1616

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

3. Computation results and discussion In order to confirm the effectiveness of the present method, FE-TSCM, in the sensitivity analysis of a rectangular plate, the authors compare the results of FE-TSCM with those of Fox’s method [3]. In Fox’s method, there is one formulation obtaining the sensitivities of eigenvalues and two formulations obtaining the sensitivities of eigenvectors. Therefore, we made three computational programs by MATLAB [17], which are called FE-TSCM, Fox 1, and Fox 2 in this study. The computational program, Fox 1, uses the formulation 1 defined by Fox for calculating the sensitivities of eigenvectors. Fox 2 uses the formulation 2 (see Appendix C). In the procedure calculating the sensitivities of eigenvalues and eigenvectors, we need the results of the free vibration analysis, which are eigenvalues and eigenvectors. In this study, FE-TSCM uses the bisection method when solving the frequency equation, Eq. (17), for obtaining eigenvalues. Fox 1 uses the determinant search method introducing the bisection method when solving the frequency equation of Eq. (C.1), for obtaining eigenvalues. When calculating the sensitivity of the ith eigenvector, FE-TSCM and Fox 1 need only the ith eigenvector among many eigenvectors. However, Fox 2 need many eigenvectors from the first eigenvector to the rth eigenvector for the computational accuracy of the sensitivity of the ith eigenvector, as shown in Eq. (C.6). In this study, the authors made the computer program, Fox 2, obtain and use all eigenvectors in calculating the sensitivity of the ith eigenvector. Fox 2 uses the function of MATLAB for eigenvalue analysis, EIG [17], for obtaining all eigenvalues and eigenvectors. Because the matrices used in Fox’s method have a lot of zero elements, the authors made effective computational program considering the point. If we do not consider the point, the computational program using Fox’s method requires very much computational time and storage. In Section 3.1, results calculating the eigenvalues, eigenvectors, and their derivatives for a simply supported rectangular plate by using FE-TSCM are compared to those obtained using Fox’s method. In Section 3.2, when the design variables are the thicknesses of the plate elements, the results computing the sensitivities of the first eigenvalues for the rectangular plate with clamped edges are calculated and illustrated in detail by using FE-TSCM. 3.1. Simply supported rectangular plate Computation model 1 is a simply supported rectangular plate. The physical and geometrical properties of the plate are indicated in Table 1. The rectangular plate is divided into 9  6, 15  10, 21  14, 30  20, and 45  30 mesh patterns, respectively. For example, Fig. 3 shows that the plate is divided into a 9  6 mesh pattern, which has a total of 54 plate elements, 70 nodes, and total 210 degrees-of-freedom. Eigenvalues and eigenvectors of computation model 1 calculated by FE-TSCM coincided well with those of the finite element method, as in the previous study [13]. When the design variables are the thicknesses of each plate element, the derivatives of eigenvalues calculated by using FE-TSCM coincided well with those of Fox’s method. Table 2 shows the derivatives of the first three eigenvalues (l1 ¼3.0525E3, l2 ¼1.1084E4, l3 ¼2.8821E4) obtained by both methods for computation model 1 when the plate is divided into a 9  6 mesh pattern, as shown in Fig. 3. In Table 2, ti denotes the thickness of the plate element i. From Table 2, we know that modifying the thickness of the plate elements 27 and 28 is more effective than the other plate elements in order to change the first eigenvalue significantly. To change the second eigenvalue significantly, modifying the thicknesses of the plate elements 15, 16, 39, and 40 is effective. For the third eigenvalue, modifying the thicknesses of the plate elements 26 and 29 is effective. In order to verify the accuracy of the derivative values of eigenvalues, when the thickness of a plate element changes, the reanalyzed eigenvalues are compared with the estimated eigenvalues obtained by using the first-order derivatives of the eigenvalues. Table 3 shows the original eigenvalues, the reanalyzed eigenvalues, and the estimated eigenvalues based on the sensitivity analysis results when the thickness of the plate element 27 is increased by 10%, 20%, and 30%, respectively. Because the estimated eigenvalues are very similar to the reanalyzed eigenvalues in spite of using only the first-order eigenvalue derivatives, we can confirm the accuracy of the sensitivity analysis for the eigenvalue of the rectangular plate using FE-TSCM. We calculated the sensitivities of the eigenvectors by using FE-TSCM and Fox’s method, respectively. In order to verify the accuracy of the sensitivities of the eigenvectors obtained by both methods, we compared the reanalyzed eigenvectors Table 1 Physical and geometrical properties of rectangular plate. Item

Value

Unit

Length (a) Width (b) Thickness (t) Young’s modulus (E) Density (r) Poisson’s ratio (n)

3.0 2.0 0.01 73.0 2770 0.32

m m m GPa kg/m3 Dimensionless

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

1617

Fig. 3. Rectangular plate with 9  6 mesh pattern. Table 2 Derivatives of the first three eigenvalues of computation model 1 with respect to thickness of plate element. Plate element (i)

@l1/@ti

@l2/@ti

@l3/@ti

1, 6, 49, 53 2, 5, 50, 53 3, 4, 51, 52 7, 12, 43, 48 8, 11, 44, 47 9, 10, 45, 46 13, 18, 37, 42 14, 17, 38, 41 15, 16, 39, 40 19, 24, 31, 36 20, 23, 32, 35 21, 22, 33, 34 25, 30 26, 29 27, 28

1.7650E4 1.0109E4 0.2568E4 1.4179E4 1.0669E4 0.7159E4 0.8862E4 1.1527E4 1.4192E4 0.4185E4 1.2281E4 2.0377E4 0.2338E4 1.2579E4 2.2820E4

6.3545E4 4.0364E4 1.7182E4 2.5012E4 4.1499E4 5.7986E4 1.1630E4 4.1894E4 7.2157E4 4.5515E4 4.0895E4 3.6275E4 7.0666E4 4.0154E4 0.9642E4

0.5849E5 0.1882E5 0.5849E5 0.6818E5 0.8574E5 0.6818E5 0.8302E5 1.8827E5 0.8302E5 0.9608E5 2.7843E5 0.9608E5 1.0123E5 3.1404E5 1.0123E5

Table 3 Change of the eigenvalues of computation model 1 for variation of thickness of plate element 27. Order

Original

Reanalysis

Estimation

(a) 10% increase 1 2 3

3.0525E3 1.1084E4 2.8821E4

3.0740E3 1.1093E4 2.8919E4

3.0754E3 1.1093E4 2.8922E4

(b) 20% increase 1 2 3

3.0525E3 1.1084E4 2.8821E4

3.0922E3 1.1102E4 2.9010E4

3.0982E3 1.1103E4 2.9023E4

(c) 30% increase 1 2 3

3.0525E3 1.1084E4 2.8821E4

3.1068E3 1.1111E4 2.9091E4

3.1210E3 1.1113E4 2.9125E4

with the estimated eigenvectors using the first-order sensitivities calculated by FE-TSCM and Fox’s method, respectively. When the thickness of the plate element 27 is increased by 30%, Table 4 shows the z-directional displacement components in the original first eigenvector, the reanalyzed first eigenvector, and the estimated first eigenvectors obtained by FE-TSCM, Fox 1, and Fox 2, respectively. The number of eigenvectors selected in Fox 2 for calculating the derivative of the first eigenvector is equal to that of the system degrees-of-freedom. For effective comparison, the data of Table 4 were

1618

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

Table 4 Change of the z-directional displacement components in the first eigenvector of computation model 1 when the thickness of plate element 27 is increased by 30%. Node

Original

Reanalysis

FE-TSCM

Fox 1

Fox 2

1, 64 2, 65 3, 66 4, 67 5, 68 6, 69 7, 70 8, 57 9, 58 10, 59 11, 60 12, 61 13, 62 14, 63 15, 50 16, 51 17, 52 18, 53 19, 54 20, 55 21, 56 22, 43 23, 44 24, 45 25, 46 26, 47 27, 48 28, 49 29, 36 30, 37 31, 38 32, 39 33, 40 34, 41 35, 42

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.1736 0.3008 0.3473 0.3008 0.1736 0.0000 0.0000 0.3264 0.5653 0.6527 0.5653 0.3264 0.0000 0.0000 0.4397 0.7616 0.8794 0.7616 0.4397 0.0000 0.0000 0.5000 0.8660 1.0000 0.8660 0.5000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.1782 0.3089 0.3573 0.3100 0.1793 0.0000 0.0000 0.3338 0.5783 0.6692 0.5815 0.3367 0.0000 0.0000 0.4476 0.7724 0.8948 0.7812 0.4533 0.0000 0.0000 0.5070 0.8617 1.0000 0.8862 0.5153 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.1794 0.3109 0.3596 0.3121 0.1805 0.0000 0.0000 0.3357 0.5816 0.6731 0.5850 0.3389 0.0000 0.0000 0.4493 0.7753 0.8985 0.7850 0.4558 0.0000 0.0000 0.5079 0.8610 1.0000 0.8894 0.5178 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.1793 0.3109 0.3596 0.3120 0.1805 0.0000 0.0000 0.3356 0.5815 0.6729 0.5849 0.3388 0.0000 0.0000 0.4492 0.7752 0.8983 0.7848 0.4557 0.0000 0.0000 0.5078 0.8610 1.0000 0.8892 0.5177 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.1793 0.3109 0.3595 0.3120 0.1805 0.0000 0.0000 0.3356 0.5815 0.6729 0.5849 0.3388 0.0000 0.0000 0.4492 0.7752 0.8983 0.7848 0.4557 0.0000 0.0000 0.5078 0.8610 1.000 0.8892 0.5177 0.0000

Table 5 Total computational times for obtaining the sensitivities of the first eigenvalue and the first eigenvector according to mesh patterns (s). Mesh

Fox 1

Fox 2

FE-TSCM

96 15  10 21  14 30  20 45  30

0.110 1.890 14.250 113.187 –

0.109 2.062 5.484 175.531 –

0.219 1.125 4.485 23.985 200.14

normalized by considering the maximum value among the z-directional displacement components of the first eigenvector as one. In Table 4, the results of Fox 1 coincided well with those of Fox 2, except for node 11 and node 60; the results of FETSCM are very similar to those of Fox 1 and Fox 2. Because the estimated eigenvectors obtained by using FE-TSCM are similar to the reanalyzed eigenvectors in spite of using only the first-order eigenvector derivatives, we can verify the accuracy of the sensitivity analysis for the eigenvector of the rectangular plate using FE-TSCM. According to various mesh patterns, Table 5 shows total computational times taken by FE-TSCM and the Fox’s methods for obtaining the sensitivities of the first eigenvalue and the first eigenvector through the free vibration analysis and the sensitivity analysis. The computational times shown in Table 5 are the sum of both the free vibration analysis and the sensitivity analysis. Because design variables are the thicknesses of all plate elements, the number of design variables equals the number of plate element.

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

1619

Table 5 show FE-TSCM is superior to Fox 1 and Fox 2 from the viewpoint of the computational time, when the partition number of plate elements is increased. The reasons why the sensitivity analysis using FE-TSCM is fast are as follows: (i) The derivative of the eigenvector can be obtained by computing once with the obtained eigenvalue derivative, like procedures (6) and (7) of Section 2.3.5. (ii) FE-TSCM does not need the calculation of large-size matrices in comparison with the standard FEM and Fox’s method, because it uses the transmission rules of small-size matrices. (iii) FE-TSCM does not need the other eigenvalues and eigenvectors except the first eigenvalue and eigenvector, for obtaining derivatives of the first eigenvalue and eigenvector, like Fox 1. Table 6 shows detail computational times taken by FE-TSCM and Fox’s methods for obtaining the sensitivities of the first eigenvalue and the first eigenvector in mesh pattern 30  20. When obtaining only the sensitivity of eigenvalue, Fox’s method is superior to FE-TSCM from the viewpoint of the computational time, because Fox’s method can reduce effectively the size of the problem from the total degrees-of-freedom of the plate system to the degrees-of-freedom of a plate element when a design variable is the thickness of a plate element. For example, if a design variable is the thickness of the jth plate element, the elements of the matrices KnS and MnS of Eq. (C.2) are zeros except those differentiating the jth element stiffness and mass matrices with respect to the design variable. However, when obtaining only the sensitivity of eigenvector, FE-TSCM is superior to Fox’s method from the viewpoint of the computational time. In general, when a dynamic system is modeled by the discrete system, if the number of elements increases, we can enhance the accuracy of the dynamic analysis. Therefore, when the number of rectangular plate elements increases, we can obtain more accurate eigenvalues and eigenvectors and their derivatives. In order to examine in detail the derivative of the eigenvalue with respect to the thickness of the plate element, when the plate is divided into a 45  30 mesh pattern, the derivatives of the first three eigenvalues (l1 ¼3.1055E3, l2 ¼ 1.1476E4, l3 ¼2.9395E4) were computed by using FE-TSCM, as illustrated in Figs. 4–6. The height of the bar illustrated at the position of the plate element in Figs. 4–6 indicates the derivative value of the eigenvalue with respect to the thickness of the corresponding plate element. Table 6 Detail computational times for obtaining eigenvalue, eigenvector and their sensitivities in mesh pattern 30  20 (s). Method

Fox 1

Computational time for free vibration analysis Computational time for sensitivity analysis Sensitivity of the first eigenvalue Sensitivity of the first eigenvector Total computational time

Fox 2

FE-TSCM

13.797

149.485

1.110

0.109 99.281 113.187

0.094 25.952 175.531

14.047 8.828 23.985

Eigenvalue derivative

1200 1000 800 600 400 200 0

45

30

40 35

25 30

20

25 15

20 15

10 5 Y

10 5

X

Fig. 4. Derivative of the first eigenvalue for computation model 1 by FE-TSCM.

1620

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

Eigenvalue derivative

4000 3000 2000 1000 0

45

30

40 35

25 30

20

25 15

20 15

10 5

10 5

X

Y

Fig. 5. Derivative of the second eigenvalue for computation model 1 by FE-TSCM.

Eigenvalue derivative

15000

10000

5000

0

45

30

40 35

25 30

20

25 15

20 15

10 5 Y

10 5

X

Fig. 6. Derivative of the third eigenvalue for computation model 1 by FE-TSCM.

In the case of mesh pattern 9  6 as shown in Fig. 3, the sizes of the system mass matrix (Ms) and the system stiffness matrix (Ks) are 210  210, respectively. The maximum sizes of the derivative matrices of the system mass matrix and the system stiffness matrix with respect to a design variable (Mns and Kns ) are 210  210, too. When considering the symmetry band matrix for compact computational memory in Fox’s method, the sizes of the matrices Ms and Ks become 210  27, respectively. However, in the case of the same mesh pattern, the sizes of the stiffness coefficient matrices ðS^ i Þ and the n sensitivity stiffness coefficient matrices ðS^ i Þ are 21  21 in FE-TSCM. Because of transfer calculation of the matrices with the small size, FE-TSCM does not require large computational memory simultaneously. However, much computational memory is required in Fox’s method, because Fox’s method use simultaneously large matrices with the total degrees-of-freedom of system. Therefore, we know that FE-TSCM is superior to Fox’s method in the view point of the computational memory. As shown in Table 5, we cannot achieve the sensitivity analysis in 45  30 mesh patterns by Fox’s method, because of the excessive use of the computational memory.

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

1621

3.2. Rectangular plate with clamped edges 3.2.1. Rectangular plate with one clamped edge Computation model 2 is a cantilever rectangular plate with one clamped edge, as shown in Fig. 7(a). The physical and geometrical properties of the plate are equal to those of computation model 1. When the rectangular plate is divided into a 45  30 mesh pattern and the design variables are the thicknesses of the 1350 plate elements, results computing the derivatives of the first eigenvalue (l1 ¼ 35.8287) of computation model 2 using

Eigenvalue derivative

Fig. 7. Rectangular plate structures with clamped edges: (a) computation model 2; (b) computation model 3; (c) computation model 4; and (d) computation model 5.

30 20 10 0 −10 45

30

40 35

25 30

20

25 15 Y

20 15

10 5

X

10 5

Fig. 8. Derivative of the first eigenvalue for computation model 2 by FE-TSCM.

1622

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

FE-TSCM are illustrated in Fig. 8. The height of the bar illustrated at the position of the plate element in Fig. 8 indicates the derivative value of the first eigenvalue with respect to the thickness of the corresponding plate element. Fig. 8 shows that the derivative values of the first eigenvalue for computation model 2 are large in the vicinity of the clamped edge. If the thickness of a plate element in the vicinity of the clamped edge increases, the first eigenvalue becomes much larger. On the other hand, the derivative values of the first eigenvalue in the vicinity of the free edge are negative. Therefore, if the thickness of a plate element in the vicinity of the free edge increases, the first eigenvalue of computation model 2 becomes lower.

3.2.2. Rectangular plate with both clamped edges Computation model 3 is a rectangular plate with both clamped edges, as shown in Fig. 7(b). The physical and geometrical properties of the plate are equal to those of computation model 1.

Eigenvalue derivative

1500

1000

500

0 45

30

40 35

25 30

20

25 15

20 15

10

Y

5

X

10 5

Fig. 9. Derivative of the first eigenvalue for computation model 3 by FE-TSCM.

Eigenvalue derivative

5000 4000 3000 2000 1000 0 45

30

40 35

25 30

20 Y

25 15

20 15

10 5

X

10 5

Fig. 10. Derivative of the first eigenvalue for computation model 4 by FE-TSCM.

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

1623

Eigenvalue derivative

15000

10000

5000

0 45

30

40 35

25 30

20

25 15 Y

20 15

10 5

X

10 5

Fig. 11. Derivative of the first eigenvalue for computation model 5 by FE-TSCM.

When the rectangular plate is divided into a 45  30 mesh pattern and the design variables are the thicknesses of the 1350 plate elements, results computing the derivatives of the first eigenvalue (l1 ¼ 1.4686E3) of computation model 3 using FE-TSCM are illustrated in Fig. 9. Fig. 9 shows that the derivative values of the first eigenvalue for computation model 3 are large in the vicinity of both clamped edges. Therefore, if the thickness of a plate element in the vicinity of both clamped edges increases, the first eigenvalue of computation model 3 becomes larger. 3.2.3. Rectangular plate with three clamped edges Computation model 4 is a rectangular plate with three clamped edges, as shown in Fig. 7(c). The physical and geometrical properties of the plate are equal to those of computation model 1. When the rectangular plate is divided into a 45  30 mesh pattern and the design variables are the thicknesses of the 1350 plate elements, results computing the derivatives of the first eigenvalue (l1 ¼ 2.1229E3) of computation model 4 using FE-TSCM are illustrated in Fig. 10. We can know that the derivative values of the first eigenvalue for computation model 4 are large in the vicinity of the corners connecting the clamped edge and the free edge. 3.2.4. Rectangular plate with four clamped edges Computation model 5 is a rectangular plate with four clamped edges, as shown in Fig. 7(d). The physical and geometrical properties of the plate are equal to those of computation model 1. When the rectangular plate is divided into a 45  30 mesh pattern, results computing the derivatives of the first eigenvalue (l1 ¼1.1137E4) of computation model 5 using FE-TSCM are illustrated in Fig. 11. Fig. 11 shows that the derivative values of the first eigenvalue for computation model 5 are large in the vicinity of the long clamped edge. 4. Conclusions This paper presents an effective algorithm using FE-TSCM in order to achieve the sensitivity analysis for free vibration of the rectangular plate structure. The effectiveness of FE-TSCM was confirmed through the sensitivity analysis of the eigenvalues and the eigenvectors for the simply supported the rectangular plate structure by FE-TSCM and Fox’s method. In addition, when the design variables were the thicknesses of the plate elements, the results computing the sensitivities of the first eigenvalues for the rectangular plate structures with clamped edges were illustrated in detail by using the computation program introducing the concept of FE-TSCM. The present study determined that FE-TSCM has merits from the viewpoint of computational time and memory in the sensitivity analysis for free vibration of the rectangular plate structure. Therefore, the authors recommend using FE-TSCM when carrying out the free vibration sensitivity analysis of the rectangular plate structure with a large number of plate elements. Appendix A. The stiffness coefficient matrices are the symmetric matrices ~ i and M ~ i are the mass and stiffness matrices for the strip i, which can be obtained respectively by The matrices K assembling the mass and stiffness matrices of the rectangular plate elements in the strip i. Because the mass and stiffness ~ i and M ~ i become the symmetric matrices of the rectangular plate elements are symmetric matrices, the matrices K

1624

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

matrices. Because the strip i, A~ i and If the matrix S^

i

~ i and M ~ i are the symmetric matrices, the submatrices of the dynamic stiffness matrix of the matrices K C~ i , are the symmetric matrices. is the symmetric matrix, the matrix Gi becomes the symmetric matrix, because G ¼ S^ þ A~ . The matrix i

i

i

T T T 1 ~ T ~ ~ T 1 T ~ T T T ~T ~T T B~ i Vi is the symmetric matrix, because B~ i Vi ¼ ðB~ i ÞðG1 i ÞðB i Þ ¼ fðB i Þ ðGi Þ ðB i Þ g ¼ fðB i ÞðGi ÞðB i Þg ¼ ðB i V i Þ . The T matrix Si þ 1 becomes the symmetric matrix, because Si þ 1 ¼ C~ i þ B~ Vi . Therefore, the symmetry of the stiffness coefficient i

matrix is always kept in the field transmission rule. ^ i and the point If the matrix Si is the symmetric matrix, the matrix S^ i becomes the symmetric matrix, because S^ i ¼ Si þ P stiffness matrix Pi is a diagonal matrix. Therefore, the symmetry of the stiffness coefficient matrix is always kept in the point transmission rule, too. Because the point stiffness matrix P1 is a diagonal matrix, the stiffness coefficient matrix S^ 1 becomes the symmetric matrix, because S^ 1 ¼ P1 . After first obtaining the matrix S^ 1 from the equation, S^ 1 ¼ P1 , the stiffness coefficient matrices, Si and S^ i , are obtained by applying in the field and point transmission rules consecutively. Therefore, all the stiffness coefficient matrices of FE-TSCM become the symmetric matrices. Appendix B. Linear function f(k*) If a design variable is the thickness of a plate element in the strip k, the derivatives of the submatrices of the dynamic stiffness matrix of the strip k with respect to the design variable are as follows: 2 n 3 n A~ i B~ i ~ n lM ~ n ln M ~ k ði ¼ kÞ: 4 nT 5¼K (B.1) n k k B~ i C~ i The derivatives of the submatrices of the dynamic stiffness matrix of the other strips except the strip k become 2 n 3 n A~ i B~ i ~ i ðiak,i ¼ 1,2,. . .,mÞ: 4 nT 5 ¼ ln M (B.2) n B~ i C~ i n n n From Eqs. (B.1) and (B.2), we know that the elements in the matrices A~ i , B~ i , and C~ i are always the first degree * polynomial of l . When the design variable is the thickness of the plate element in the strip k, the derivative of the point stiffness matrix Pi becomes

Pni ¼ 0

ði ¼ 1,2,3,. . .,mþ 1Þ:

n

n

(B.3) n

n

n

From Eqs. (28) and (B.3), the matrix S^ 1 is null matrix. Because S^ 1 ¼ 0 and the elements of the matrices A~ 1 , B~ 1 , and C~ 1 are the first degree polynomial of l*, the elements of the matrices Gn1 and Vn1 of Eq. (25) become the first degree polynomial of n

l*. Then, the elements of the matrix S2 are the first degree polynomial of l*, from Eq. (24). Because the matrix Pn2 is null n

matrix, the elements of the matrix S^ 2 are the first degree polynomial of l*, from Eq. (22). n From Eqs. (22) and (24), the elements of the matrix S^ 3 are the first degree polynomial of l*, because the matrix Pn3 is null n n n matrix and the elements of the matrices A~ , B~ , and C~ are the first degree polynomial of l*. Then, the elements of the 2

2

2

n n n n matrices S^ 4 , S^ 5 , . . ., S^ m , and S^ m þ 1 are the first degree polynomial of l*, because the matrix Pni is null matrix and the n n n elements of the matrices A~ i , B~ i , and C~ i are the first degree polynomial of l*. Therefore, the function f(l*) of Eq. (31) is the first degree polynomial and the linear function of unknown eigenvalue derivative l* in the numerical calculation finding l*.

Appendix C. Derivatives of eigenvalue and eigenvector proposed by Fox The eigenvalue problem for an undamped system in structural dynamics is KS X ¼ lMS X,

(C.1)

where l and X are the eigenvalue and the eigenvector of the system, respectively. The N  N matrices KS and MS are the system stiffness matrix and the system mass matrix, respectively. Their order N corresponds to the number of total degrees-of-freedom of the system; for example, N is equal to 3  n  (mþ1) in the rectangular plate structure of Fig. 1. n Fox proposed Eq. (C.2) in order to obtain the derivative of the ith eigenvalue, li , for a dynamic system in the previous paper [3],

lni ¼ XTi ½KnS li MnS Xi , n

(C.2)

n

where the N  N matrices KS and MS are differentiating the system stiffness and mass matrices with respect to the design variable d. li and Xi are the ith eigenvalue and eigenvector of the system, respectively. The ith eigenvector is normalized as XTi MS Xi ¼ 1:

(C.3)

M.-S. Choi, J.-H. Byun / Journal of Sound and Vibration 332 (2013) 1610–1625

1625

Fox proposed two types of formulations in order to obtain the derivative of the eigenvector of the system. In formulation 1, the derivative of the eigenvector is computed by n T n Xni ¼ Z1 S ½DS DS þMS Xi Xi MS Xi ,

(C.4)

where ZS ¼ DS DS þ 2MS Xi XTi MS , DS ¼ KS li MS , n

DnS ¼ KnS li MnS li MS :

(C.5)

In formulation 2, the derivative of the eigenvector is computed by Xni ffi

r X

aij Xj ,

(C.6)

j¼1

where r is the number of modes used in calculation,   aij ¼ XTj ½KnS li MnS Xi = li lj   aii ¼  XTi ½MnS Xi =2

ðiajÞ, ði ¼ jÞ:

(C.7)

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]

¨ R.T. Haftka, Z. Gurdal, Elements of Structural Optimization, Kluwer Academic Publishers, Dordrecht, 1992. S.S. Rao, Engineering Optimization, Wiley, New York, 1996. R.L. Fox, M.P. Kapoor, Rates of change of eigenvalues and eigenvectors, AIAA Journal 6 (12) (1968) 2426–2429. R.B. Nelson, Simplified calculation of eigenvector derivatives, AIAA Journal 14 (9) (1976) 1201–1205. B.P. Wang, Improved approximate methods for computing eigenvector derivatives in structural dynamics, AIAA Journal 29 (6) (1991) 1018–1020. K.J. Bathe, Finite Element Procedures, Prentice-Hall, New Jersey, 1996. N.S. Sehmi, Large Order Structural Eigenanalysis Techniques (Algorithms for Finite Element Systems), Ellis Horwood Limited, Chichester, 1989. M.A. Dokainish, A new approach for plate vibration: combination of transfer matrix and finite element technique, Transactions of the ASME Journal of Engineering for Industry 94 (1972) 526–530. H. Xue, A combined finite element-stiffness equation transfer method for steady state vibration response analysis of structures, Journal of Sound and Vibration 265 (2003) 783–793. L.K. Abbas, M. Lei, X. Rui, Natural vibrations of open-variable thickness circular cylindrical shells in high temperature field, Journal of Aerospace Engineering 23 (3) (2010) 205–212. B. Rong, X. Rui, G. Wang, Modified finite element transfer matrix method for eigenvalue problem of flexible structures, Journal of Applied Mechanics 78 (2011). 021016-1-7. E.C. Pestel, F.A. Leckie, Matrix Methods in Elastomechanics, McGraw-Hill, New York, 1963. M.S. Choi, Free vibration analysis of plate structures using finite element-transfer stiffness coefficient method, Journal of Mechanical Science and Technology 17 (6) (2003) 805–815. T. Kondou, T. Ayabe, A. Sueoka, Transfer stiffness coefficient method combined with concept of substructure synthesis method (linear free and forced vibration analyses of a straight-line beam structure), JSME International Journal (Series C) 40 (2) (1997) 187–196. M.S. Choi, T. Kondou, Y. Bonkobara, Development of free vibration analysis algorithm for beam structures by combining Sylvester’s inertia theorem and transfer stiffness coefficient method, Journal of Mechanical Science and Technology 17 (6) (2003) 805–815. D.H. Moon, M.S. Choi, Development of sensitivity analysis algorithm for a straight-line beam structure by the transfer stiffness coefficient method, JSME International Journal (Series C) 46 (1) (2003) 138–144. /http://www.mathworks.comS.