Shape design sensitivity analysis using closed-form stiffness matrix for hierarchic triangular elements

Computers d Slructures Vol. 43, No. Printed in Gnat Britain.

I, pp. 69-75,

004%7949192 S5.00 + 0.00 Pergamon Press plc

1992

SHAPE DESIGN SENSITIVITY ANALYSIS USING CLOSED-FORM STIFFNESS MATRIX FOR HIERARCHIC TRIANGULAR ELEMENTS B. P. WANG, D. BABU, R. V. NAMBIARand K. L. LAWRENCE Department of Mechanical Engineering, University of Texas at Arlington, Box 19023, Arlington, TX 76019, U.S.A. (Received 26 February

1991)

Abstract-Shape design sensitivity analysis for finite element analysis using straight-sided hierarchic triangular elements is considered in this paper. Using a recently developed closed-form stiffness matrix, the sensitivity analysis is carried out by direct differentiation of the system equations for finite element analysis. That is, the term dK/aa required in sensitivity analysis is formulated in the closed form. This is in contrast with prior formulations which utilize finite differences (the semianalytical formulation) and

numerical integration to obtain these terms. Both the displacement sensitivity analysis and the strain and stress sensitivity analysis are addressed in this paper.

1NTRODUCTION

entiating eqn (1) with respect to nodal coordinate and rearranging we get

Design sensitivity analysis [l] is a procedure to study the effects of the changes on the overall design. There are two basic methods for efficient sensitivity calculation, namely, analytic and semianalytic methods. The analytic method uses the differentiation of shape functions with respect to the nodes for the explicit differentiation of the stiffness matrix. The semianalytic sensitivity calculation uses the finite difference method for numerical evaluation of the derivatives of the stiffness matrix. The shortcoming of using this method is the uncertainity in the choice of perturbation step size Aai. If the step size is too large, truncation errors may be excessive. These can be thought of as errors due to retention of only the lowest-order terms of a Taylor series representation of a perturbed function. If the step size is too small, condition errors may occur. Condition errors are due to inaccuracy in the calculation of the displacement and roundoff errors in the finite difference calculation. For the analytic sensitivity analysis calculation, numerical integration using Gauss quadrature rule may be used and can become computationally expensive. Sensitivity analysis in the closed form for straight-sided elements is obtained easily using a recently developed closed-form stiffness matrix [2, 31. DISPLACEMENTSENSITIVITY

cc,

The above equation can also be written as

where K, is the element stiffness matrix and F, is the force acting on the element. Note that in eqn (3) e E a means the summation is only over those elements that are affected by a. If the loaded boundary is not affected by the design changes, F, is not a function of a aF

2

aa

co.

(4)

So, eqn (2) becomes

ANALYSIS

The global finite element equation can be written as and the displacement KU=F,

sensitivity is defined as

(1)

where K is the global stiffness matrix, U is the displacement vector and F is the load vector. Differ-

(6) 69

70

B. P. WANG

Usually, the term dK,/da is evaluated numerically. Closed form expressions for the stiffness matrix and load vectors for straight-sided triangular elements have been derived. This paper demonstrates that expressions for shape sensitivity for such elements can also be derived in closed form. For example, for plane stress and plane strain analysis, the element stiffness matrix is given by

(7) where t is the element thickness and XX, XY, YY are matrices whose elements are functions of the nodal coordinates (see Appendix A for details). The closed-form expression for ar<,/&x is

et al.

Note that U, is the nodal element.

STRESS-STRAIN

SENSITIVITY

ANALYSIS

The strain in an element is given by E = PRU,,

(9)

where P is the matrix of inverse of the Jacobian [see eqn (A7)] and R is a matrix of derivatives of shape functions given by N A’

0

N&2 ’

R =

[

1

0

N,L,

0

N42

,

(10)

where N is the shape function in terms of the area coordinates L, , L2 and L,. The strain sensitivity can be computed from

$=ERU.+P~U,+PR$.

(11)

Note that dP/au has been evaluated while computing aK,/ao: (see Appendix A). R is a function of area coordinates and since the strain is evaluated at node point, aRIaa = 0. Therefore eqn (9) becomes

a6 ap au

-=aRU,+PRz.

(12)

of this

Stress sensitivity

The expression for stress is

u

=De

(13)

where D is the elasticity matrix. Therefore, the stress sensitivity is

aa aa

-DE

(14)

which can be calculated using eqn (12).

EXAMPLES

With this information, eqn (6) can be used to compute the displacement sensitivity. The analysis procedure is given in Appendix A.

displacement

AND RESULTS

Displacement sensitivity is calculated in closed form for an equilateral triangle and the results are compared with data found using the finite difference method. Three sets of boundary conditions are applied to a single triangular element and the nodal coordinates are changed in each case. For example, in case 1, node 1 is free to move in both x- and y-directions, while node 2 is fixed in x-direction and node 3 is fixed in both directions. The pressure is applied on the edge opposite to the free node 2 as shown in Fig. l(a). A small change in x- or y-direction is applied to the free node and the displacement sensitivity is calculated for each change. The results calculated by closed form and by the finite difference method are compared and are given in the tables of Fig. 1. Similarly, the boundary conditions and the load conditions are changed as shown in Fig. 2(a) and Fig. 3(a) and the results for changes in the coordinates of node 2 and node 3 are given in Figs 2 and 3. In another example, for the same triangular element, instead of applying a pressure, point loads are applied on the free node and the closed form sensitivity results are compared and tabulated in Figs 4-6. In both the cases the results obtained by closed-form calculation matches the converged finite difference method results.

CONCLUSION

Closed-form formulation for design sensitivity analysis in plane stress analysis using hierarchic triangular elements is derived in this paper. Good results are demonstrated for constant strain elements. In view of the accuracy and efficiency of the closed-form formulations, it is expected that the present sensitivity analysis formulation could be used together with boundary parameterization in shape optimization,

yl

Fig. 1. Comparison

V5yl

Vl,yl

w

*Delta D = 0.1

Ul* xl Vl 1x1 V2A

1

?’

D

3

2

0.0062

0.0017

0.0062

0.0016

0.0017

0.0062

O.tXKH

Method OJD

0.0017 0.0062

0.0017

0.0062

0.0017

0.0062

0.0016

0.0061

of displacement sensitivity analysis results using closed finite difference method for case 1.

0.05 D 0.0000

Method 05D o.oooo

Difference D O.WJO

of node 1.

o.cQ17

0.0062

0.05 D 0.0000

of node 1.

X

Nodal

Finite 10 D 0.0000

Sensitivity for a change in Y coordinate

O.OOMl

Differewe D

O.OOCN3

Finite 10 D

Sensitivity for a change in X coordinate

I

”

Example-l: with Pressun Load:

Fig. 2. Comparison

V3,y2

0.0062 and

V2,p

form analysis

u2,y2 0.0017

*Delta D = 0.1

V3.d

V2N2

u2,fi

Closed Form O.oooO

0.0017

0.0062

Closed Form O.CCHIO

Coordmales

3

1

2

o.fnIO12

0.00009

Difference D O.ooooO

o.ooo12

O.oooo8

Method 05D O.ooooO

-O.OcKlO!3

-0.00009 o.oooo3

-0.OOtM9 o.OooO3

form

o.oOoO3

-O.OOMB

analysis

0.05 D O.tWOOO

of node 2.

0.00012

O.oooo8

0.05 D O.OOf00

of node 2.

of displacement sensitivity results using closed difference method for case 2.

O.oooo3

Method O.SD O.ooooO Ditference D O.OOtWO

Finite 10 D O.owoO

Sensitivity for a change in Y coordinate

0.00012

O.OOUO9

Finite 10 D O.OOfQO

Sensitivity for a change in X coordinate

Roundary Conditions and Loading for Case 2

J7

Ix

Y

and finite

o.oOOO3

-O.OOOOg

Closed Form O.OOtNN

o.OQO12

O.OOOOg

Closed Form O.wooO

-

IT

3

1

\

O.W300

O.ODOOO

0.00016

u3,9

v3xI3

O.ooooO

o.ooM2

u3.y3

v37y3

o.om61

O.WlOO

OSlW45

Difference D

o.ooo61

O.CQOOO

0.00044

Method 0.5 D

form

0.00060

O.OOQOO

0.00044

0.05 D

of node 3.

0.00016

O.OCC00

0.05 D -0.00344

of node 3.

of displacement sensitivity results using closed difference method for case 3.

0.00045

Finite 10 D

VLy3

Fig. 3. Comparison

0.00016

O.OOOOO

Method 0.5 D -0.00044

Sensitivity for a change in Y coordinate

0.00016

D -0.00044

v&3

DitTerence

Sensitivity for a change in X coordinate

Finite 10 D -0.00044

*Delta D = 0.1

2

Boundary Conditions and Laading for Case 3

Lx

Y

analysis

and

O.ooo60

0.0O000

0.00044

Closed Form

OJXQ16

O.ooooO

Closed Form -0.00044

finite

Vl

yl

4.21

3.92

Fig. 4. Percentage

v2, yl

n

“l,yl

2.00 0.00

1.18

3.80

“1x1

Yxl v2, xl

2.00 0.00 I

Boundary Condition

error

I

# 1:

b

3

1 FX

2

I

0.30

0.96

SmP 0.50 0.00

0.15

0.48

Sii 0.25 0.00

between

2.18

2.00

1.00 0.00

closed-form

1.10

1.01

Step 0.50 0.00

0.22

0.20

0.10 0.00

of node 1

0.06

0.19

0.10 0.00

of node 1

X

=

=

Fx Fy

=

L

=

1.0

1.0

10

0.3

method

0.02

0.01

0.01 0.00

0.005

0.018

0.01 0.00

1.0 nu

E

25.00 60.36 73.30

Coordinates 25.00, 60.36. 12.06.

1 2 3

=

Nodal

results and finite difference

0.55

0.51

Size 0.25 0.00

Sensitivity for a change in Y coordinate

0.59

1.91

1.M) 0.00

Sensitivity for a change in X coordinate

IY

Example - 2: With Point Loads

results.

0.01

0.01

0.005 0.00

0.002

0.008

0.005 O.txll

4.61

2.00 0.00 4.22

“3, l0

WY2 V&y2

Fig. 5. Percentage

1.21

1.98

u2,fl

“5x2

“3, y2

2.00 0.00

I

3

f

FY

FX

X

=

t

Fy = 1.0

Fx = 1.0

1.U

0.48

Step 050 0.00

0.56

0.24

si 0.2.5 0.00

between

0.60

1.00 0.00 2.10

closed-form

0.30

SteP 050 0.00 1.04 0.06

0.10 0.00 0.21

of node 2

0.22

0.10

0.10 0.00

0.00

0.01 0.21 0.05

0.023

0.009

0.01 0.016

0.00

0.005 0.54 0.05

0.04

0.02

0.005 0.01

results and finite difference method results.

0.15

Size 0.25 0.00 0.52

Sensitivity for a change in Y coordinate

2.28

0.97

1.00 0.00

1.0

= 0.3

= 1.0

Coordinates

nu

E

Nodal

Sensitivity for a change in X coordinate of node 2

error

Boundary condition # 2-

7.07 Fig. 6. Percentage

y3

y3

0.00 “39

“39

0.06

0.13

“19y3

0.00

0.00

0.03

0.00

Step 0.50 0.81

0.02

0.00

Size 0.25 0.41

error

between

3.48

0.00

1.00 1.66

closed-form

1.72

0.00

Step 0.50 0.82

results

0.86

0.00

Size 0.25 0.41

1.0

X

0.34

0.00

0.10 0.16

of node 3

0.01

OSKI

0.10 0.16

73.30 60.36 I

25.00

method

0.03

0.00

0.01 0.01

0.002

0.001

0.01 0.027

and finite difference

’

Sensitivity for a change in Y coordinate

1.00 1.61

2.00 3.38

=

= 1.0

= 0.3

= 1.0

Fy = 1.0

Fx

t

nu

E

25.00. (12.06. 60.36.

1

Coordinates

32

Nodal

Sensitivity for a change in X coordinate of node 3

FY

2.00 3.16

1

Boundary Condition # 3:

results.

0.02

0.00

O.lXE 0.00

0.002

0.001

0.005 0.027

l% P.

74

WANG

REFERENCES 1. E. J. Haua, K. K. Choi and V. Komkov. Design Sensitivity Analysis of Structural Systems. Academic Press, New York (1986). 2. R. V. Nambiar, Closed form expressions for hierarchic triangular and tetrahedral elements, Doctoral dissertation, University of Texas at Arlington, TX (1989). 3. K. L. Lawrence, R. V. Nambiar and B. Bergmann, Closed form stiffness matrices and error estimators for plane hierarchic triangular elements. Int. J. Numer. Meth. Engng, to be published. 4. Chang So Han, Shape optimization in 2-D continuum structure by using p-version finite element method. Doctoral dissertation, University of Texas at Arlington, TX (1989). 5. I. Babuska, I. N. Katz and B. A. Szabo, Hierarchic families for the p-version of the finite element method. In Advances in Computer Metho& for Partial Differential Equations (Edited by R. Vichnevetsky and R. S. Stapleman), pp. 272-236. IMACS (1979). 6. A. Peano, Hierarchies of conforming finite elements for plane elasticity and plate bending. Comput. Math. Applic. 2, 21 l-224 (1976).

et d.

In the above formulation only the G matrix is a function of nodal coordinates G = PWP, where

(A6)

and

I-=

(A 1) [:I=[;:

where xii = x, - xj, and J is the determinant of the Jacobian, which is twice the area of the triangle, and given by @g)

J = XI3Y23 - X32Y3,.

From eqns (A2) and (A3), we observe that the computation of aK/da involves the calculation of aZ/r?a, where Z = XX, XY or YY. It should be noted that in (A3), only the terms G, and J are functions of a. The term dJ/da can be evaluated using (A8) and aG,,/aa can be evaluated using (A5). Therefore

:][::I?

where N is the shape function in terms of the area coordinates L, , L, and L, (see Appendix B). The nodal displacement vector in the x-direction 6u and the nodal displacement vector in the y-direction 6v are usually interspaced alternately in one vector, but are shown separately here for convenience. The stiffness matrix in closed form can be derived to be

xx

K=t

047) --

APPENDIX A: DETAILS OF THE CLOSED-FORM STIFFNESS MATRIX In the assumed displacement formulation of the stiffness matrices, the displacement vector at any point within the element is given as

(As)

[ XY7

XY

t.42)

YY1 ’

K,a = t

XX,,

XY,,’

XY;

YY,,

VW

and

+ JG,,,A + G&D + 87 + G,,,,C) XY,,= J,,(G,,A + G,,B+ G,,BT+ G,C)

where the matrices XX, XY, YY are obtained from

+ JG,,,J + G,,,J + G,,,X XX = J(G,,A

+ G,,(B + Br) +

XY = J(G,,A

+ G,,B + G&3’+

YY = J(G,,A

+ G,(B

+ G,,,,C)

G,,C) YY,, = J,.(G,,A + G,(E +

B’) + G,C)

G,,C)

+ B’) + GaC)

+ JG,,J

+ G,,,,@ + 8’) + G,,,C).

(A101

(A3)

and A, B, C are matrices that contain only constants and given by [2]

From the above equations stiffness matrix differentiation can be found out in closed form.

APPENDIX B: HIERARCHIC SHAPE FUNCTIONS FOR TRIANGULAR ELEMENTS

B=

A family of hierarchic triangular elements suggested by Babuska et al. [S] and Peano [6], have shape functions as follows. The degree of freedom corresponding to the first three shape functions are the displacement at the vertices

N;, N,Lz dR s

C=

s

N,TL2 N,L2di2.

644)

N(i)=L,,

i=l,2,3.

(Al 1)

Shape design sensitivity analysis Table Bl. Shape functions for hierarchic triangular elements Degree of polynomial 1 2 3 4

Shape functions ‘5 2L, L2

2L2JJ3

L?L,-LiL, LiL,LiL,

LiL,-LiL, L;L, - L;L,

L2

L3

L: L, - L: L, L:L, -LiL,

2L34

L,L,

L, L2L3 L,L,

L,dL,=b.

In addition, three shape functions of the form L$L,-L,(-L,)‘,

L$L,-L,(-L,)‘, L/,L, - L,(-L,)’

(A12)

are added at the mid-sides of the triangles of the order p > 2, for each j = 1,2, . ,p - 1. Additional internal ‘bubble mode’ shape functions, each containing the product L, 4 L,

so that they vanish at the boundaries, are added at the centroid to complete the polynomial basis. Peano [6] has shown that the r entries in N together form a linearly independent set that spans the space of polynomials of order p, where r = (p + 2)(p + 1)/2.

(Al31

For ease of reference, the shape function for polynomial upto degree 4 are explicitly summarized in Table Bl.