A numerical method of calculating first and second derivatives of dynamic response based on Gauss precise time step integration method

European Journal of Mechanics A/Solids 29 (2010) 370–377

Contents lists available at ScienceDirect

European Journal of Mechanics A/Solids journal homepage: www.elsevier.com/locate/ejmsol

A numerical method of calculating first and second derivatives of dynamic response based on Gauss precise time step integration method Qimao Liu a, b, *, Jing Zhang c, Liubin Yan b a

Department of Civil Engineering, Guangxi University of Technology, Liuzhou 545006, PR China College of Civil and Architecture Engineering, Guangxi University, Nanning 530004, PR China c Department of Mechanical Engineering, University of Alaska Fairbanks, Fairbanks AK 99775, United States b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 10 June 2009 Accepted 18 November 2009 Available online 26 November 2009

This paper describes an accurate and efficient method for calculating the first and second derivatives of dynamic response with respect to design variables for linear structural systems subjected to transient loads. An efficient algorithm to calculate the dynamic responses and their first and second derivatives is formulated based on Gauss precise time step integration method. The algorithm is achieved by direct differentiation and only a single dynamic analysis is required. Several numerical examples are comparatively demonstrated using the new developed method, analytical method, and central difference method. The results show that the new method is highly accurate compared with the analytical approach and is more efficient than the central difference method. Ó 2009 Elsevier Masson SAS. All rights reserved.

Keywords: Precise time step integration method Dynamic response first derivative Dynamic response second derivative Sensitivity analysis Hessian matrix

1. Introduction Calculation of first derivatives of structural dynamic response with respect to design variables is the prerequisite when the gradient-based methods are used in optimization design of the structures subjected to transient loads. It is also called sensitivity analysis in structural optimization. Many algorithms have been developed in this regard (Keulen et al., 2005; Kirsch et al., 2005; Bogomolni et al., 2006). There are two broad categories: direct differentiation method (DDM) and adjoint variable method (AVM) (Haftka and Adelman, 1989). In DDM, the direct integration method (Sousa et al., 1997) and the mode superposition method (Dutta and Ramakrishnan, 2000) are often used to calculate the design sensitivity. To do a sensitivity analysis using DDM, the number of second order differential equations to be solved is equal to the number of design variables. Thus, DDM is not efficient when the number of design variables is large. In AVM, an augmented response function is defined to calculate sensitivity of a response. The adjoint variables are introduced to define the augmented response function (Arora and Cardoso, 1992). The AVM is more efficient compared with the DDM when the number of design variables is large and the number of active constraints is small. However, the AVM is not useful when the constraints in a nonintegral form, i.e., pointwise constraints, are

* Corresponding author. Tel.: þ86 772 2686010; fax: þ86 772 2686010. E-mail address: [email protected] (Q. Liu). 0997-7538/$ – see front matter Ó 2009 Elsevier Masson SAS. All rights reserved. doi:10.1016/j.euromechsol.2009.11.006

involved. The AVM is only applied to the constraints in an integral form (Dutta and Ramakrishnan, 1998) and the sensitivity of the equivalent function can be calculated using AVM (Kocer and Arora, 2002). The pointwise constraints can be converted into the constraints in integral form using the Dirac delta function. However, if the number of time grid points to be considered is large, application of AVM to pointwise constraints is not efficient. Generally, the gradient-Hessian matrix-based method is more efficient compared with the gradient-based method when the Hessian matrix can be calculated efficiently. The gradient-Hessian matrix-based methods, e.g., Newton’s Method and Marquardt’s method (Reklaitis et al., 1983), require the information of dynamic responses, their first and second derivatives to construct optimal algorithms. Although many researchers have calculated the first derivatives of dynamic response, there is little work published on dynamic response second derivative analysis. Second derivative analysis is more complicated than first derivative analysis. For dynamic analyses, the precise time step integration method (PTSIM), in theory, can achieve any order accuracy at integration points (Zhong and Williams, 1994). When the equations of motion are homogeneous linear differential equation, the numerical result obtained by PTSIM is the exact solution at integration points. However, when the equations of motion are non-homogeneous linear differential equations, the accuracy of PTSIM depends on the fitted accuracy of non-homogeneous terms. In structural dynamics, the non-homogeneous terms are dynamic loads, such as wind, ocean wave and earthquake, which are often very complicated.

Q. Liu et al. / European Journal of Mechanics A/Solids 29 (2010) 370–377

Generally, it is difficult to achieve the accurate fitting curve of the dynamic loads. Gauss precise time step integration method (GPTSIM) (Wang and Zhou, 2004) was developed to overcome this difficulty because GPTSIM does not require curve fitting of the dynamic loads. A first-order sensitivity analysis method of structural transient linear and non-linear heat conduction was developed using PTSIM by Chen et al. (2003). However, very little work has been published on using GPTSIM to calculate the first and second derivatives of dynamic response of the structures subjected to transient loads. The purpose of this paper is to develop an accurate and efficient method for calculating the first and second derivatives of dynamic response based on GPTSIM for linear structural systems subjected to transient loads. We formulate an efficient algorithm to calculate dynamic responses and their first and second derivatives. The algorithm is achieved by direct differentiation and only a single GPTSIM-based dynamic analysis is required. The detailed computational procedure is provided. Moreover, in order to evaluate the accuracy and efficiency of the new method, we also formulate the analytical method and the central difference method to study the dynamic responses and their first and second derivatives. Using numerical examples, we compare the accuracy of the new method with the analytical solution and the efficiency of the new method with the central difference method. 2. Problem description The structural design variable vector is denoted by d, the number of design variable is m. Consider the equations of motion for a linear system subjected to dynamic forces

€ðtÞ þ C xðtÞ _ Mx þ KxðtÞ ¼ f ðtÞ

(1)

xð0Þ ¼ x0 ; _ xð0Þ ¼ x_ 0

(

vxð0Þ=vdi ¼ x0di ; _ _ 0di vxð0Þ=vd i ¼ x

P ¼ M x_ þ Cx=2

(6)

Eq. (1) becomes:

v_ ¼ Hv þ r

(7)

where

v ¼ ½x

P T

(8)

r ¼ ½0

f T

(9)

H ¼



M 1 C=2 CM 1 C=4  K

M 1 CM 1 =2



(10)

The general solution of the state equation Eq. (7) is

Ht

vðtÞ ¼ e v0 þ

Zt

eHðtsÞ rðsÞds

(11)

0

The dynamic load duration of time is divided evenly into several time intervals with time step Dt. The time station is denoted by tk ¼ kDt (k ¼ 0,1,2,.). The relationship of two adjacent time station is tkþ1 ¼ tk þ Dt. According to Eq. (11), the relationship of v(tk) and v(tkþ1) is

vðtkþ1 Þ ¼ TðDtÞvðtk Þ þ

Ztkþ1

eHðtkþ1 sÞ rðsÞds

(

  v2 xð0Þ= vdi vdj ¼ x0di dj   _ vdi vdj ¼ x_ 0di dj v2 xð0Þ=

ði; j ¼ 1; 2; /; mÞ

ð2Þ

where M, C, and K are mass matrix, damping matrix, and stiffness _ €ðtÞ are unknown nodal matrix, respectively. x(t), xðtÞ and x displacement, velocity and acceleration vectors, and f(t) is the load vector. Rayleigh damping is used to derive the damping matrix in this work.

C ¼ a1 M þ a2 K

(3)

where

  a1 ¼ 2ð21 u2  22 u1 Þu1 u2 = u22  u21

(4)

  a2 ¼ 2ð22 u2  21 u1 Þ= u22  u21

(5)

where u1 and u2 are the first and second natural frequency of the structure, respectively. 21 and 22 are the first and second mode damping ratios, respectively. The problem is to calculate the structural dynamic responses and their first and second derivatives with respect to design variables.

(12)

tk

where

TðDtÞ ¼ eH Dt

with the following initial conditions:



371

(13)

It is critical to calculate the exponential matrix T (Dt) accurately in the PTSIM. We transform Eq. (13) into

TðDtÞ ¼ eH Dt ¼



eH s

2 N

(14)

where s ¼ Dt=2N . According to references (Zhong and Williams, 1994; Wang and Zhou, 2004), the computational results are accurate enough when N ¼ 20. In this work, we use N ¼ 32 in example 1 and N ¼ 20 in examples 2–3. In order to avoid the accuracy loss of the PTSIM in rounding operation of computer, the exponential matrix T (Dt) is calculated as follows: First, the exponential matrix eH s is expanded in Taylor series:

TðsÞ ¼ eH s ¼ I þ T 0 T 0 ¼ H s þ ðH sÞ2 =2! þ / þ ðH sÞL =L!



(15)

where L is truncation order and I is identity matrix. According to references (Zhong and Williams, 1994), the computational results are accurate enough when L ¼ 4. In this paper, we use L ¼ 10 in example 1 and L ¼ 4 in examples 2–3. Second, the following equation is iterated for N times to obtain TN:

T i ¼ 2T i1 þ T i1 ,T i1 ði ¼ 1; 2; /; NÞ

(16)

Finally, the exponential matrix T(Dt) can be calculated by the following equation:

3. Dynamic response analysis with GPTSIM

TðDtÞ ¼ I þ T N

In this section, GPTSIM is presented. We make transformation as follows:

The integration term in Eq. (12) can be achieved by using Gauss integration method:

(17)

372

Ztkþ1

Q. Liu et al. / European Journal of Mechanics A/Solids 29 (2010) 370–377

eHðtkþ1 sÞ rðsÞds ¼

tk

n Dt X

2

  v CM 1 C=4  K wi TðDtð1  xi Þ=2Þ,rðtk þ Dtð1 þ xi Þ=2Þ

n Dt X

2

wi TðDtð1  xi Þ=2Þ,rðtk þ Dtð1 þ xi Þ=2Þ

i¼1

ð19Þ

  v  CM 1 =2

In this section, the formulas for first and second derivatives of dynamic response are derived based on the GPTSIM. Differentiating Eq. (6) with respect to the design variable di, we have

vvðtkþ1 Þ vTðDtÞ vvðtk Þ ¼ vðtk Þ þ TðDtÞ vdi vdi vdi n Dt X vTðDtð1  xi Þ=2Þ þ w ,rðtk þ Dtð1 þ xi Þ=2Þð21Þ 2 i¼1 i vdi where, the first derivatives of exponential matrix with respect to design variable di, i.e., vT (Dt)/vdi and vTðDtð1  xi Þ=2Þ=vdi , can be calculated as follows:

vdi

vH ¼ eH Dt Dt, vdi

Further differentiating Eq. (20) with respect to the design variable dj, we have

v2 P v2 M vM vx_ vM vx_ v2 x_ 1 v2 C ¼ þ þM þ x x_ þ vdi vdj vdi vdj vdi vdj vdj vdi vdi vdj 2 vdi vdj þ

1 vC vx 1 vC vx 1 v2 x þ þ C 2 vdi vdj 2 vdj vdi 2 vdi vdj

,rðtk þ Dtð1 þ xi Þ=2Þ



 v  M 1 C=2

vdi

In Eq. (29), the first derivatives of exponential matrix with respect to design variable have been obtained in Eq. (22). The second derivatives of exponential matrix with respect to design variables can be calculated as follows:

v2 TðDtÞ vH vH v2 H ¼ eH Dt Dt 2 , þ eH Dt Dt, vdi vdj vdj vdi vdi vdj

  v  M 1 C=2

vdi

(30)

where

2

  v2  M 1 C=2

6 6 v2 H vdi vdj   ¼ 6 6 2 vdi vdj 4 v CM 1 C=4  K

  v2 M 1 v2



3

7 7 7 7  CM 1 =2 5

vdi vdj



v M 1



  v2  M 1 C=2

3

7 7 vdi  7 7 v  CM 1 =2 5

vdi vdj

¼ M 1

(23)



vdi

vdi vdj

(31)

vdi vdj

vM 1 C M 1 vC ¼ M 1 M  2 vdi vdi 2

(24)

  v2 M 1 vdi vdj

vM 1 ¼ M 1 M vdi

vM 1 vM 1 C v2 M 1 C M M M þ M 1 vdj vdi 2 vdi vdj 2

M 1 v2 C vM 1 vM 1 C  M 1 M M 2 vdi vdj vdi vdj 2

1 vM 1 vC M þ M 1 2 vdi vdj

where

  v M 1

ð29Þ

In Eq. (31),

6 6 vH vdi   ¼ 6 6 vdi 4 v CM 1 C=4  K

vdi

ð28Þ

v2 vðtkþ1 Þ v2 TðDtÞ vTðDtÞ vvðtk Þ vTðDtÞ vvðtk Þ ¼ vðtk Þ þ þ vdi vdj vdi vdj vdj vdi vdi vdj n 2 2 v vðtk Þ Dt X v TðDtð1  xi Þ=2Þ þTðDtÞ þ w vdi vdj 2 i¼1 i vdi vdj

(22)

In Eq. (22),

2

(27)

(20)

Omitting high-order small quantity in Eq. (19) and differentiating Eq. (19) with respect to the design variable di, we have

vTðDtÞ ¼ vdi

vC M 1 C 1 vM 1 M  M vdi 2 2 vdi

Further differentiating Eq. (21) with respect to the design variable dj, we have

4. Formulas for first and second derivatives of dynamic response

vP vM vx_ 1 vC 1 vx ¼ þ xþ C x_ þ M vdi vdi vdi 2 vdi 2 vdi

¼ 

vdi

Eq. (19) is the iteration format of Gauss precise time step integration.

veH Dt

vC M 1 C vM 1 C CM 1 vC vK M   CM 1 þ 4 vdi vdi vdi 4 vdi 4 (26)

ð18Þ

where n is the number of the Gauss integral points, xi is the coordinates of the integral points, wi is weighting coefficient and T (Dt (1  xi)/2) is exponential matrix. Substitute Eq. (17) and Eq. (18) into Eq. (12), we have

  þ o Dt 2n

vdi

i¼1

  þ o Dt 2n

vðtkþ1 Þ ¼TðDtÞvðtk Þ þ

¼

(25)

1 vM 1 vC þ M 1 M 2 vdj vdi ¼ M 1

ð32Þ

vM 1 vM 1 v2 M 1 M M  M 1 M vdj vdi vdi vdj

þ M 1

vM 1 vM 1 M M vdi vdj

ð33Þ

Q. Liu et al. / European Journal of Mechanics A/Solids 29 (2010) 370–377

  v2 CM 1 C=4  K vdi vdj

Step 6 Calculate the state vector v (tkþ1) at each time station (k ¼ 0,1,2,.) by solving Eq. (19). Step 7 Calculate the first derivatives of state vector vv (tkþ1)/vdi, (i ¼ 1,2,.,m) at each time station (k ¼ 0,1,2,.) by solving Eq. (21). Step 8 Calculate the second derivatives of state vector v2v(tkþ1)/ (vdivdj) (i,j ¼ 1,2,.,m) at each time station (k ¼ 0,1,2,.) by solving Eq. (29).

v2 C M 1 C vC vM 1 C ¼ M 1 M  vdi vdj 4 vdi vdj 4 þ

vC M 1 vC vC 1 vM 1 C  M M vdi 4 vdj vdj vdi 4

þ CM 1

vM 1 vM 1 C M M vdj vdi 4

 CM 1

v2 M 1 C M vdi vdj 4

þ CM 1

vM 1 vM 1 C M M vdi vdj 4

 CM 1

vM 1 1 vC vC M 1 vC M þ vdi 4 vdj vdj 4 vdi

6. Example 1 A 2-bar truss is shown in Fig. 1. Node 2 is subjected to both horizontal dynamic load, Fh(t) ¼ Ph sin ut, and vertical dynamic load, Fv(t) ¼ Pv sin ut. Where Ph ¼ Pv ¼ 1500 N and u ¼ 3p/4. A total duration of 3 s and an incremental time step of 0.01 s are considered in the procedures. We use typical steel material properties, i.e., elastic modulus E ¼ 210 GPa, material density r ¼ 7850 kg/m3 and Poisson’s ratio m ¼ 0.3.

C vM 1 vC CM 1 v2 C  M 1 M þ 4 vdi vdj 4 vdj vdi  v

2



  CM 1 =2 vdi vdj

¼  

v2 K vdi vdj

373

ð34Þ

Fv ( t )

v2 C M 1 vC 1 1 vM 1 M þ M vdi vdj 2 vdi 2 vdj

Fh ( t ) 0

2

x

1 vC 1 vM 1 C 1 vM 1 vM 1 M M þ M M M 2 vdj vdi 2 vdj vdi

C v2 M 1 C 1 vM 1 vM 1  M 1 M þ M M M 2 vdi vdj 2 vdi vdj

2m

y

(35)

3

1 5. Computational procedure of the first and second derivatives of dynamic response Section 4 gives the formulas, i.e., Eq. (21) and Eq. (29), for calculating the first and second derivatives of the dynamic response, i.e., displacement. This section provides the detailed computational procedure. Procedure of calculating first and second derivatives of dynamic response: Step 1 Calculate stiffness matrix K, damping matrix C, mass matrix M and their first and second derivatives with respect to design variables. Calculate the inverse matrix M1. Step 2 Calculate matrix H, its first and second derivatives with respect to design variables by solving Eq. (10), Eq. (23) and Eq. (31), respectively. Step 3 Set time step Dt and the number of Gauss integration points, calculate the exponential matrices, i.e., TðDtÞ; TðDtð1  x1 Þ=2Þ; TðDtð1  x2 Þ=2Þ; /; TðDtð1  xn Þ=2Þ, by solving Eqs. (15)–(17). Step 4 Calculate the first and second derivatives of the exponential matrices by solving Eq. (22) and Eq. (30), respectively. Step 5 Initial calculations: 5.1 Calculate P (0) by solving Eq. (6). Obtain the initial value of state vector, i.e., vð0Þ ¼ ½ xð0Þ Pð0Þ T . 5.2 Calculate vPð0Þ=vdi by solving Eq. (20). Obtain the initial value of the first derivatives of state vector, i.e., vvð0Þ=vdi ¼ ½ vxð0Þ=vdi vPð0Þ=vdi T . 5.3 Calculate v2 Pð0Þ=ðvdi vdj Þ by solving Eq. (28). Obtain the initial value of the second derivatives of state vector, i.e., v2 vð0Þ=ðvdi vdj Þ ¼ ½ v2 xð0Þ=ðvdi vdj Þ v2 Pð0Þ=ðvdi vdj ÞT .

2m Fig. 1. Two-bar truss.

The cross section of the truss element is given in Fig. 2. de and te is the outer diameter and wall thickness, respectively. The design variable vector is d ¼ ½d1 ; d2 ; t 1 ; t 2 ; x2 ; y2 T and the current design point is dc ¼ [90,90,12,12,0,0]T (mm). All the initial conditions in Eq. (2) are zeros. Following the procedures described in Section 5, we can calculate the first and second derivatives of the horizontal and vertical displacements with respect to the design variables at the current design point.

te

de

Fig. 2. Element cross section.

6.1. Analytical method In order to evaluate the accuracy of the new method, we also formulate an analytical method to calculate the dynamic responses and their first and second derivatives with respect to design variables.

374

Q. Liu et al. / European Journal of Mechanics A/Solids 29 (2010) 370–377

In this example, the number of dynamic freedom is 2. The horizontal and vertical displacements of node 2 are unknowns. The dynamic load vector is

f ðtÞ ¼



Ph sinut Pv sinut



¼ P 0 sinut

(36)

v2 B1 =ðvdi vdj Þ and v2B2/(vdivdj) can be obtained by solving Eq. (44). Then the dynamic response second derivative v2x/(vdivdj) is achieved from substituting v2B1/(vdivdj) and v2B2/(vdivdj) into Eq. (45). 6.2. Results and discussion of the new method and the analytical solution

Eq (1) becomes

(37)

_ Because the initial conditions are x(0) ¼ 0 and xð0Þ ¼ 0, we only need to calculate the steady state response. Assuming that the steady state response solution is

x ¼ ½ U21

V21 T sinut þ ½ U22

V22 T cosut

¼ B1 sinut þ B2 cosut

(38)

Substitute Eq. (38) into Eq. (37), we have

  K  u2 M B1  uCB2 ¼ P 0   uCB1 þ K  u2 M B2 ¼ 0

(39)

Eq. (39) can be written as:



K  u2 M uC

uC K  u2 M



B1 B2





P0 0

¼



(40)

B1 and B2 can be obtained by solving Eq. (40). Then the dynamic response x is achieved from substituting B1 and B2 into Eq. (38). Eq. (40) can be written as follows:

~B ~ ¼ P ~ K

(41) "

~ ¼ where K

K  u2 M uC

#

 uC ~ ¼ ,B K  u2 M



B1 B2



~ ¼ and P



 P0 . 0

The calculations were performed on a notebook computer (CPU frequency: 498 MHz, RAM: 128 MB) with Window XP operating system. The computational time of solving the problem is about 11 s with the new method, about 2 s with the analytical method, which indicates that the analytical method is more efficient than the new method. In order to evaluate accuracy of the new method, the results of the two methods are shown in the same figure. The relative error of the results of the new method is also analyzed. One of the first derivative-time curves obtained with the two methods and the relative error-time curve of the new method is shown in Figs. 3 and 4, respectively. One of the second derivative-time curves and the relative error-time curve of the new method is shown in Figs. 5 and 6, respectively. From Figs. 3 and 5, it is observed that the time history curves of first and second derivatives, achieved by the new method and the analytical method, are similar results. From Figs. 4 and 6, the maximum relative error of the numerical results achieved by the new method is w2  103. The relative error 4.E-08

First derivative

€ þ C x_ þ Kx ¼ P 0 sinut Mx

Analytical Method

0.E+00 -2.E-08 -4.E-08 0.0

Differentiating Eq. (41) with respect to design variable di, we have

New Method

2.E-08

0.5

1.0

1.5

2.0

2.5

3.0

Time(s)

~ ~ ~ vB ¼  vK B ~ K vdi vdi

vx vB1 vB ¼ sinut þ 2 cosut vdi vdi vdi

(43)

vB1/vdi and vB2/vdi can be obtained by solving Eq. (42). Then the dynamic response first derivative vx/vdi is achieved from substituting vB1/vdi and vB2/vdi into Eq. (43). Further differentiating Eq. (42) with respect to the design variable dj, we have:

~ v2 K vdi vdj

¼

v K vdi vdj

v M  u2 vd vd i

2 C uvdvi vd j

(44) #

2

j

C uvdv vd i

v2 K vdi vdj

j 2

v M  u2 vd vd i

.

j

Further differentiating Eq. (43) with respect to the design variable dj, we have: 2

2

2

v x v B1 v B2 ¼ sinut þ cosut vdi vdj vdi vdj vdi vdj

2.E-03 1.E-03 0.E+00 -1.E-03 -2.E-03 -3.E-03 0.0

0.5

1.0

1.5

2.0

2.5

3.0

Time(s) Fig. 4. Relative error of first derivative, horizontal displacement of node 2, with respect to d1. (Between new method and analytical method).

(45)

New Method

Second derivative

where

2

3.E-03

1.E-07

2~ 2~ ~ ~ ~ ~ ~ v B ¼  v K B ~  vK vB  vK vB K vdi vdj vdj vdi vdi vdj vdi vdj 2

Relative error of first derivative

vC # " vK  u2 vM u vdi vdi vdi . vC vK vM u  u2 vdi vdi vdi Differentiating Eq. (38) with respect to design variable di, we have ~ vK where ¼ vdi

"

Fig. 3. First derivative, horizontal displacement of node 2, with respect to d1.

(42)

5.E-08

Analytical Method

0.E+00 -5.E-08 -1.E-07 0.0

0.5

1.0

1.5

2.0

2.5

3.0

Time(s)

Fig. 5. Second derivative, horizontal displacement of node 2, with respect to d1, d2.

Relative error of second derivative

Q. Liu et al. / European Journal of Mechanics A/Solids 29 (2010) 370–377

375

Fv ( t )

3.E-03 2.E-03

Fh ( t )

1.E-03

10

x 42

0

0.E+00

y

-1.E-03

41

-2.E-03 -3.E-03 0.0

11

0.5

1.0

1.5

2.0

2.5

3.0

Fig. 6. Relative error of second derivative, horizontal displacement of node 2, with respect to d1, d2. (Between new method and analytical method).

36

40

31

7. Examples 2 and 3

5

In order to evaluate the efficiency of the new method, in this section, we obtain the first and second derivatives of dynamic response using both the new method and the central difference method. In order to obtain the more accurate results, all the stepsizes of design variables are assumed to be 0.001. All simulations thereafter were performed on a desktop PC (CPU frequency: 1 GHz, RAM: 256 MB) with Window XP operating system.

16

6

25

28 15

24

23

17

20

16 19

12

15

14

2

5

13 18

9

10

18 17

7

2 1

29

13.5m

3

7.1. Example 2

33 14

22

4 11

34

30

38 13

27

6 21

39

35

7 26

43 12

32

8

decreases rapidly and does not propagate in the time step integration process. Therefore, the new method has almost the same accuracy as the analytical method. The main reason of high accuracy is that the formulas for the first and second derivatives of dynamic response are derived by direct differentiation and are based on the two reliable theories, i.e., the finite element method and GPTSIM. The new method is not efficient compared to the analytical method. However, the analytical method can only be used in simple structures because the analytical solution of the complex structures is difficult to be achieved. The new method can also calculate the first and second derivatives of dynamic response for both simple and complex structures.

44 37

9

Time(s)

A 45-bar truss system is shown in Fig. 7, which has a much larger degree of dynamic freedom, design variable number and element number than the 2-bar system. Node 10 is subjected to both horizontal and vertical dynamic load, Fh(t) ¼ Fv(t) ¼ 1500 sin (3pt/4) (N). A total duration of 3 s and an incremental time step of 0.01 s are considered in the procedures. The material properties are same as those in example 1. The design variables are cross-sectional areas (i.e., A1wA45) and nodal coordinates (i.e., x2wx19). The truss system is symmetric along y-axis. Therefore, the number of independent design variable is 36. The current design point: Ai ¼ 0.8 m2 (i ¼ 1w45), xi ¼ 1.5 m (i ¼ 11w19), xi ¼ 1.5 m (i ¼ 2w10). We assume that all the initial conditions in Eq. (2) are zeros. We calculate the first and second derivatives of all nodal displacements with respect to the design variables at the current design point with both the new method and the central difference method. The results of the first derivatives of dynamic response are shown in Fig. 8. The results of the second derivatives of dynamic response are shown in Fig. 9. It is observed that the time history curves of first and second derivatives, achieved by the new method and the central difference method, are similar results.

45

8 19

4

3 20

1 3m Fig. 7. 45-bar truss.

7.2. Example 3 A three-storey, two-bay planar concrete frame is shown in Fig. 10. The mass of each storey is assumed to be 30 tonnes. Node 4 is subjected to both horizontal and vertical dynamic load, Fh(t) ¼ Fv(t) ¼ 1500 sin (3pt/4) (N). A total duration of 3 s and an incremental time step of 0.01 s are considered in the procedures.

376

Q. Liu et al. / European Journal of Mechanics A/Solids 29 (2010) 370–377 Table 1 Design variables and current design point of frame (units: m).

4.E-10 New Method

First derivative

2.E-10

Element type

Central Difference Method

Storey level

Member group

0.E+00

Column

0.5

1.0

1.5

2.0

2.5

C5 C6

b9 ¼ 0.9 b8 ¼ 0.9

h9 ¼ 0.9 h8 ¼ 0.9

2nd

C3 C4

b7 ¼ 0.9 b6 ¼ 0.9

h7 ¼ 0.9 h6 ¼ 0.9

1st

C2 C1

b5 ¼ 0.9 b4 ¼ 0.9

h5 ¼ 0.9 h4 ¼ 0.9

3rd

B3

b3 ¼ 0.9

h3 ¼ 0.9

2nd

B2

b2 ¼ 0.9

h2 ¼ 0.9

1st

B1

b1 ¼ 0.9

h1 ¼ 0.9

3.0

Time(s)

Depth he

3rd

-2.E-10 -4.E-10 0.0

Current design point Width be

Fig. 8. First derivative, horizontal displacement of node 10, with respect to A2. Beam

4.E-10

New Method Central Difference Method

2.E-10

1.E-06

0.E+00 -2.E-10 -4.E-10 -6.E-10 0.0

0.5

1.0

1.5

2.0

2.5

3.0

Time(s)

First derivative

Second derivative

6.E-10

Fig. 9. Second derivative, horizontal displacement of node 10, with respect to A1, A2.

5.E-07

New Method Central Difference Method

0.E+00 -5.E-07 -1.E-06 0.0

0.5

1.0

1.5

2.0

2.5

3.0

Time(s)

Fv ( t ) Fh ( t )

Fig. 11. First derivative, horizontal displacement of node 4, with respect to b2.

B3 4

B3 8

12 1.E-06

C5

C6 B2

9

C5

B2 11

3.6m

C3

3

7

2

C4 B1

5.4m

B1 10

6

2 C1

7

C1 9

5 6.6m

5.E-07

New Method Central Difference Method

0.E+00 -5.E-07 -1.E-06 0.0

0.5

1.0

1.5

2.0

2.5

3.0

Time(s)

C2 1

C3

8

Second derivative

3.6m

6.6m

Fig. 12. Second derivative, horizontal displacement of node 4, with respect to b1, b2.

response are shown in Fig. 12. It is observed that the time history curves of first and second derivatives, achieved by the new method and the central difference method, are similar results.

Fig. 10. Three-storey, two-bay planar frame.

7.3. Discussion of the new method and the central difference method Material elastic modulus: E ¼ 20 GPa, material density: r ¼ 2500 kg/m3, and Poisson’s ratio: m ¼ 0.2. The beams and columns are classified into 9 member groups shown in Table 1. The cross section of the beam element is rectangle. The design variables are the width be and depth he of beam element in the eth member group. The design variables and the current design point are also shown in Table 1. The number of design variable is 18. We assume that all the initial conditions in Eq. (2) are zeros. We calculate the first and second derivatives of all nodal displacements with respect to the design variables at the current design point with both the new method and the central difference method. The results of the first derivatives of dynamic response are shown in Fig. 11. The results of the second derivatives of dynamic

In the 45-bar truss example, the computational time is about 1217 s to achieve all the first and second derivative information using the new method. The computational time of a single dynamic analysis is about 0.6 s for the 45-bar truss with GPTSIM. For a problem with 36 design variables, central-difference first and second derivative calculations require total repetition of the analysis for 4  (362  2  36 þ 1) different design points. The central difference method requires about 3218 s to obtain all the first and second derivatives of dynamic response with respect to design variables. In the frame example, the computational time is about 129 s to complete all the first and second derivative calculations using the new method. The computational time of a single dynamic analysis

Q. Liu et al. / European Journal of Mechanics A/Solids 29 (2010) 370–377

is about 0.4 s for the beam frame with GPTSIM. For a problem with 18 design variables, central-difference first and second derivative calculations require total repetition of the analysis for 4  (182  2  18 þ 1) different design points. The central difference method requires about 537 s to obtain all the first and second derivatives of dynamic response with respect to design variables. Therefore, the two numerical examples indicate that the new method is more efficient than the central difference method. The important reason of high efficiency is that the dynamic responses, their first and second derivatives can be calculated simultaneously in only a single dynamic analysis in the new method.

8. Conclusions A new method for first and second derivatives of structural dynamic response is developed. The conclusions summarized as follows: (1) An efficient algorithm is formulated to calculate dynamic responses of linear structural systems and their first and second derivatives with respect to design variables. The algorithm is achieved by direct differentiation and only a single dynamic analysis based on Gauss precise time step integration method is required. (2) The relative error analysis shows that compared with the analytical method, the new method is highly accurate. The main reason of high accuracy is that the formulas for the first and second derivatives of dynamic response are derived by direct differentiation and based on the two reliable theories, i.e., the finite element method and Gauss precise time step integration method. (3) The results show that compared with the central difference method, the new method is more efficient to obtain the first

377

and second derivatives of dynamic response with respect to design variables. The important reason of high efficiency is that the dynamic responses, their first and second derivatives can be calculated simultaneously in only a single dynamic analysis in the new method.

References Arora, J.S., Cardoso, J.B., 1992. Variational principle for shape design sensitivity analysis. AIAA Journal 30, 538–547. Bogomolni, M., Kirsch, U., Sheinman, I., 2006. Efficient design sensitivities of structures subjected to dynamic loading. International Journal of Solids and Structures 43, 5485–5500. Chen, B.S., Gu, Y.X., Zhang, H.W., Zhao, G.Z., 2003. Structural design optimization on thermally induced vibration. International Journal for Numerical Methods in Engineering 58, 1187–1212. Dutta, A., Ramakrishnan, C.V., 1998. Accurate computation of design sensitivities for structures under transient dynamic loads with constraints on stress. Computers and Structures 66, 463–472. Dutta, A., Ramakrishnan, C.V., 2000. Accurate computation of design derivatives for plates and shells under transient dynamic loads. Engineering and Computers 17, 7–27. Haftka, R.T., Adelman, H.M., 1989. Recent developments in structural sensitivity analysis. Structural and Multidisciplinary Optimization 1, 137–151. Keulen, F.V., Haftka, R.T., Kim, N.H., 2005. Review of options for structural design sensitivity analysis. Part 1: linear systems. Computer Methods in Applied Mechanics and Engineering 194, 3213–3243. Kirsch, U., Bogomolni, M., Keulen, F.V., 2005. Efficient finite-difference designsensitivities. AIAA Journal 43, 399–405. Kocer, F.Y., Arora, J.S., 2002. Optimal design of latticed towers subjected to earthquake loading. Journal of Structure Engineering 128, 197–204. Reklaitis, G.V., Ravindran, A., Ragsdell, K.M., 1983. Engineering Optimization Methods and Applications. John Wiley & Sons, New York. Sousa, L.G., Cardoso, J.B., Valido, A.J., 1997. Optimal cross-section and configuration design of elastic-plastic structures subjected to dynamic cyclic loading. Structural Optimization 13, 112–118. Wang, M.F., Zhou, X.Y., 2004. Renewal precise time step integration method of structural dynamic analysis. Acta Mechanica Sinca 36 (2), 191–195 (in Chinese). Zhong, W.X., Williams, F.W., 1994. A precise time step integration method. Proceedings of Institution Mechanical Engineers Part C: Journal of Mechanical Engineering Science 208, 427–430.