The boundary element formulation for multiparameter structural shape optimization

The boundary element formulation for multiparameter structural shape optimization T. Burczyfiski and T. A d a m c z y k Institute of Mechanics and Fundarn en tals of Machine Design, Silesian Technical University, Konarskiego 18a, 44-101 Gliwice, Poland (Received July 1983)

A problem of optimal shaping of a free and loaded boundary has been formulated in terms of the boundary element method. A maximal stiffness criterion together with the limitation of volume of a body has been used. The problem has been solved with the application of integral optimality conditions being defined on an optimal boundary. Discretization of a boundary and optimality conditions respectively, with the aid of boundary elements, leads to an iterative procedure, from which one can determine the shape parameters generating the position of the optimal boundary, Key words: mathematical model, boundary element method, shape optimization

The determination of shape of two- and three-dimensional bodies plays an essential role in machine element design. The only universal method of obtaining a solution for arbitrary complex design systems is an application of numerical methods. The boundary element method has found broad application to analysis problems in which, as as a boundary shape, a boundary condition is known. ~-3 In such problems one should determine a distribution of statical or dynamical fields due to given boundary conditions and body forces. The problem of shape synthesis is much more complicated, because only a proportion of the parameters defining a boundary of a body is given. The remaining parameters which explicitly define the shape of a free and loaded boundary of a body should be defined in such a way that constraints have to be satisfied and a functional of shape synthesis ought to reach an extremum. Two optimization techniques, 4,s optimality conditions and dynamical programming, are most often applied to problems formulated in such a way. Optimal shaping has been solved with the aid of the finite element method by Zienkiewicz and Campbell, 6 Ramakrishnan and Francavilla 7 and Dems and Mr6z. 1° The boundary eleme at and linear programming method was described by Fu agamis through its application to the optimal control of aeat con0307-904X/85/03195-06/$03.00 @ 1985 Butterworth & Co. (Publishers) Ltd

duction phenomena. A formulation of optimal shaping of design problems with the aid of boundary elements has been presented by Burczyfiski and Adamczyk. 9 The boundary element method is a remarkably convenient and natural way of modelling optimal shaping in design. Optimality conditions is in this case an especially suitable optimization technique, because one can formulate the conditions in terms of integral expressions determined on an unknown boundary being optimized. Discretization of these conditions with the aid of boundary elements leads to a set of nonlinear algebraic equatibns, which together with a set of linear algebraic equations, obtained by discretization of boundary integral equations, enables determination of unknown optimal shape parameters. The present authors have concentrated on an application of the boundary element method to a multiparameter optimal shaping using a maximum stiffness of design as a criterion. Integral o p t i m a l i t y c o n d i t i o n s f o r l o a d e d and free boundary Consider an elastic body with domain V in Figure 1 bounded by S. Let the boundary consist of three parts: S t

Appl. Math. Modelling, 1985, Vol. 9, June

195

when tile boundary S ' = So is optimized: d

5u.

~

(6)

cb = W - - - (u . t ° ) + 2 H u . t ° - 6 o dn

.

when the boundary S' = S t is optimized; and:

Z

f I"

)",

(7)

d V = Vo

in equation (6) H is mean curvature of the S t surface, while 5g, in equation (4), is a variation of the g vector. It should be noted that the integral optimality condition (4) is defined only on boundaries which are optimized.

X1

Figure I

g vector which generates an optimal boundary of a body

on which the tractions t i = t o = a i j n j are given, S u on which the displacements II i tl 0 are given and the surface So being free from tractions: evidently S = S t U S u u S o . Components of displacement function g ( x ) = ( g i ( x ) ) (i = I, 2, 3 ; F i g u r e 1) are in the required functions, which generate the optimal boundary, having a certain reference boundary as a basis. An optimization problem is defined as follows: minimize or maximize a certain functional of synthesis of shape rr' in such a way, that the volume V of the domain in question should not exceed a certain value V0 stated previously:

Discretization of optimum conditions with the aid of boundary elements An optimal boundary is required among a family of polyhedrons with triangular sides. These sides are built up by the boundary elements ( F i g u r e 2 a ) . A local coordinate system is used with a unitary basis e i ( F i g u r e 2 b ) on every

=

rain or max rr' for V~< Vo

S 1

(1)

q

The problem may be reduced to finding a stationary condition functional: rr = r r ' - c o ( V - Vo)

(2)

where co denotes the Lagrange multiplier. Depending on the optimization criteria being used one can consider different kinds of synthesis functionals rr'. Here it is required to limit the maxixnunl design stiffness. This criterion, which is also termed the minimal strain or minimal compliance criterion, may be regarded as an expansion of energy formulations of mechanics on problems of shaping mechanical and geometrical characteristics of design. Taking potential energy as a measure of total stiffness:

.'= fw(e) dV- fu.t o ds, V

x2

a

i

(3)

03= • 3 .

St

where W(e) is a strain potential. Applying the principle of virtual work and the complementary virtual work respectively for the variable boundary stationary conditions for the functional are obtained, which can be treated as the necessary optimality conditions for the boundary: 1°'11

~57r= ,!qbn • 5g dS' = 0

(4)

s

/ x I

b

where: ~ = W--co

196

Appl. Math. Modelling, 1985, Vol. 9, June

(5)

Figure 2

(a) Distribution of surface using boundary elements; (b) local coordinate system on a triangular boundary element

boundary element. A set o f S c' (a = i, 2 . . . . . M') elements is described in the form of a sum with N components equal to the number of nodes, i.e.: N

{s

integral on application of the Ostrogradsky-Gauss theorem: f d V = - 3I f r . n d S = - 1

V

= u p=l

Figure2b).

The function which generates the optimum boundary may now be written for every S~ as follows: 1

1

S

r~.n~ dS~

~ f

=1

(13)

Sa

Taking into account that the position of an arbitrarily chosen point on the element S ~' can be presented by radiusvector r~:

where {S~} is a set of all these elements which join the same node p (p = i , 2 . . . . , N:

/

3

ra(~l, [z) = r~ + ~iei

(14)

Equation (13) may be written in the following form:

\

1

M

3 y- r 'n lS t =vo

(15)

0~=I

1

1

+ L----~~'ap + 1 ~,a=p

(8)

L2 a = 1,2 . . . . . Mo ; p = 1,2 . . . . . No. Vector shape parameters ap, a~, and a~ are defined by:

%

~

= ~2 = 0

1 when

~l=Lt,~2

2

~

ap

0

(9)

= O, ~= = L~

a ={aip } (i = 1,2, 3;p = 1,2 . . . . ,No). From equation (4) it follows that if, on the surface being optimized, there are Mo elements and No nodes, the optimal boundary will be generated by 3No shape parameters aip. The optimality condition (4) may be written as follows:

where [Sal is an area of the S a element and r~ is a radiusvector of the node p where the local coordinate system of the S '~ element is placed. Equations (12) and (15) create a set of 3No+ I nonlinear algebraical equations, from which one can determine shape parameters aip generating the optimal boundary and Lagrange multiplier co respectively. In a case of twodimensional problems one requires the optimal boundary among a family of polygons. Sides of this polygon are built by rectilinear segments of length LR which are in fact the boundary elements (Figure 3). The optimality condition after discretization with the aid of the boundary elements can then be written in the form: l;

~(Lp -- ~) np d~

0

Lp_ I

8it= f cbn . Og 8ai, dS '

1 f

+-Lp-I

S'

dP(Lp_l-~)np_ld~=O

(16)

o

. . . . . No. The condition of maintaining constant area of polygon can be expressed by:

p = 1,2

= ? lqSbf n,,' =

(1o)

Baip~g6aip dS,~ = O

1 N

By proper assembly of expressions in equation (10) one obtains:

bTr=

y. ~ p = 1 i =1 g~i

~Pn. ~g 8 a i p d S ~ i = O Baip

(11)

In equation (11) the summation convention is only valid with regard to indexes, ap denotes a number of boundary elements which join the node p. Putting the optimal boundary generating function (8) into equation (11) and taking a free choice of 8aip variation into account, the following stationary conditions for the functional 7r are obtained:

(bni(l--~ll ~ , - 5 ~ 2 ) s i n T d ~ l d ~ a = 0

(17)

2 S" r p . n p L p = A o p=l

No es

. ........

__S,ooounoory n0

/,-

%o ,,

(12) rp %

where the relation dS~i = sin 3' d~1 d/j2 has been used, while 3' is an angle between base vectors of the,local coordinate system which are situated on the boundary element S'ai and which join the node.p. The constraint (7) can be transformed into a surface

,,1 ~

,

~.,I

b()undary

Figure 3 Discretization of two-dimensional problems using boundary elements

Appl. Math. Modelling, 1985, Vol. 9, June 197

where Ao denotes the given area of the polygon. Equations (16) and (17) create a set of 2No + 1 nonlinear algebraic equations with regard to unknown shape parameters. Nonlinear algebraic equations (12) and (15) (for three-dimensional case) or (16) and (17) (for two-dimensional case) can be written in uniform matrix form:

Fm(bq) = 0

Discretizing the surface with the aid of boundary elements S a yields:

o[f

2 u(y) = ~.

~=t

(18)

U ( x , y ) t(x) dSa(x)

Sa

- f T(x,y)u(x)

where

(25)

S~

b3p - 2 = alp, b3p - I = a2p, b3p = a3p, b3No+ l = w (m, q = 1, 2 . . . . . 3No + 1 ) for three-dimensional case;

b2p = a2p, b2p- 1 = alp (Figure 3), b2No+ 1 = t..o (m, q = 1,2 . . . . . 2No + 1 ) for two-dimensional case. The nonlinear functions F,,, depend on the strain potential W. E v a l u a t i o n o f strain p o t e n t i a l using B E M Using an orthonormal local coordinate system {el, ) (Figure 2b); the vector e3' is normal to the surface of the boundary element S~. The strain potential can be written in this coordinate system as:

W(e) = pei'i,ei,/, + 2 ek'k'el'l'

( 19)

where p and X are the Lame constants. Components of the strain tensor eft/, on the boundary element S '~ can be expressed by components of the tensor el~ in the system with basis (el} as follows:

ei,], = Ai,kAi'tekt

(20)

where:

[Ai'k]=

l l

--c°ty

(sin ~,)-i

0

0

Components of the strain tensor ei/on the surface S are determined using BEM. For this, a boundary integral equation is taken into account:

The displacement vector u and the traction vector t can be approximated using the interpolation functions:

(26)

u(}) = Na(}) un t(}) = N~(}) tn

} = ( } , , }2)

where N~(}) is the matrix of the interpolation functions, u n and t n are the nodal vector displacements and tractions. As a result, from equation (25), one has a set of linear algebraic equations with regard to unknown displacements at nodes on the surface S t tO So and to nodal tractions on the boundary S u. In matrix form this set can be written as follows:

KX = Y

(27)

where K is the matrix of coefficients, X is the vector of unknown tractions and displacements at nodes on the boundary, Y is the vector of known tractions and displacements at nodes on the boundary. Solving the set of algebraic equations (27) the strain vector e can be defined in terms of the nodal displacement un as follows:

e = Bu n

(28)

where B is the gradient matrix relating the strains on the boundary element and the displacement. In the case of the free surface So there is no need to determine all of the'components of the strain tensor according to the above procedure. For the surface free from loading, So are only four of the strain tensor components: q T , e 2 ' 2 ' , e 3 ' 3' and ef2, are non-zero elements, while ea'a' one can determine for Hooke's law putting 03, 3, = 0: e3' 3, =

-?,(e,' l' + e2'2')/(2u + X). For the surface S t components of the strain tensor e f f, e2,=', el, =, are also unknown, while the rest one can determine from the relation:

= f u(x, y) t(x) dS(x) s

t o, = o f 3,= 2 p q , 3, +

-- f T(x, ), ) . ( x ) aS(x)

e l , 3, =

(21)

s where components of the fundamental solution tensors U and T are given for three-dimensional problems:

t o, = o 2 ' 3 ' = 2 U e z ' a ' + e = ' 3 ' =

1

- - t o, 2/a

lut°'

(29)

t o, = o3, 3, = 2pc3, 3, + Xek,u'

Uii(x'Y)=

167rp(1- u) l--2v

(1'~

[(3--4u)Si/+r.ir,/1 Or

3 1 + 2u r'ir'

(22)

,)

--nir, i + nir, j ]

(23)

yi)] 1/2

(24)

where: r = [(x i --yi)(xi

--

~, is Poisson's ratio.

198

Appl. Math. Modelling, 1985, Vol. 9, June

] e3'3 ' = ~

2U+ X

[t30 ' -

~.(el'l'k

-I- e2,2,)]

For optimal shaping of the boundary St a method of evaluating d/dn(u, t °) should be chosen (see equation (6)), while the evaluation of du/dn is especially difficult. Working in the local coordinate system {el'}, one has: du --=.

dn

g i,tli,----t~,3,

It follows from this equality, that one needs to determine u l ' y , u2',y and u3' 3,, while U3',3'= e3'3'. The method of evaluation is shown in the example:

= 2 e f 3 ' - - Afke3"u,k(e3' = e3) For two-dimensional problems the method of determination of the strain potential is almost the same as for three-dimensional problems. When a boundary being considered is free from loading and inside the body there is a two-dimensional strain state, only el'l' and e2, 2, are not equal to zero; e2'2' = -- X/(2/~ + ~)e1'1' (for the loaded boundary e2'2' = -- (t o, -- her' f)/(2/a + ?,) and el,2' = t°,/21a). Only el' 1' = u i', 1' = du i'/d~ is determined by straight differentiation of displacements, which are obtained from the solution of the boundary integral equation for the twodimensional problem. A similar method can be applied in the case of the two-dimensional stress state.

A l g o r i t h m o f o p t i m a l shaping A solution of equation (18) is rather difficult due to the nonlinear functions Fro, which are implicit functions of the parameters b/. These difficulties can be overcome by applying the Newton-Raphson iterative procedure in the form of:

q (m, q = 1,2 . . . . . 3No+ 1) for three-dimensional problem; (m, q = 1,2 . . . . . 2No + 1) for two-dimensional problem where Abq are increments of the parameters bq. New values of the parameters bq in (i + 1)th iteration are defined from the relation: (31)

The matrix of derivatives [~Fm/abq], which occurs in equation. (30) may be defined by differentiation of proper optimality conditions or using the direct method in terms of finite differences. In the first case the problem reduces to determination of derivatives of nodal displacements and the radius-vectors of nodes of the boundary elements r = ( x , y , z). The radiusvector of the pth node in the (i + 1)th iteration can be expressed by the relation: r~ +1)= tp-(i)+ a~')

(32)

Differentiating relation (32) with regard to the parameters bq and taking expressions which connect them with the p a r a m e t e r s ap into account relations for coordinates of the radius-vector can be obtained as follows:

x(i+p,ql)____-q~3p--2' YP,q" (i+1) =

-qt~3P- 1,

Zp,q(i+l)= 53p

(33)

where 5~i denotes the Kronecker delta. Derivatives of nodal displacements are obtained by differentiating equation (27). Then by solution of the following set of equations values of the derivatives can be found: K X . q = Y q -- K, q X

~F m

-- Fm(b 1. . . . , bq -- ~Sbq. . . . )

abq

u(, 3, = 2 e l y -- uy, t' = 2el' 3' -- A (kU3',k

by + l) = b 0") + Abq

Fm(bl . . . . . bq_ 1, bq+ 5bq . . . . )

(34)

For evaluating derivatives aFm/Ob q a simpler direct method can be applied instead:

25bq

(35)

where 5bq denotes a given small increment of the shape parameter bq. Application of such a method of evaluation requires solution of the set of equations (27) 6N times (for threedimensional problems) and 4N timeg (for two-dimensional problems) in every iteration step. It follows from equation (30) that a certain initial reference boundary is taken after which the increments Abq are determined. On this basis a shape of the new boundary is defined using the new iteration. The optimization process ends when in two succeeding steps values of the parameters bq are sufficiently close. A multiparameter procedure for optimal shaping which has been described above on the basis of integral optimality conditions and the boundary element method is shown schematically below: (1) Set up an initial shape of boundary being optimized. Perform a discretization of the boundary S in terms of boundary elements. (2) Compute: the matrices K, Y; the nodal displacements of boundary elements using equation (27). (3) Compute: the strain potential W for the boundary elements of varying shape; the matrix {Fro). (4) Compute: the derivatives of radius-vectors of nodes of the boundary elements; the derivatives of matrices K q and Y,q; the derivatives of nodal displacements of the boundary elements; the derivatives of m a t r i x [aFm/abq] or the derivatives of m a t r i x [~Frn/abq]using equation (35). (5) Compute the shape parameter increments Abq from equation (30). (6) Check whether values of Abq being obtained are greater than given evaluation accuracy, i.e.

/

q~ (Abq)2 ~< e b

If yes, go to step (7); if no, go to step (2). (7) Print the nodal coordinates of the optimal boundary elements (OBE).

Conclusions Problems of optimal shaping are much more complicated than problems of analysis. Despite the formulation of a problem for a body with linear elasticity, after discretization of the optimality conditions the set of nonlinear algebraic equations has been obtained to solve. Application of an iterative algorithm allows for determination of the shape parameters, which generate a configuration of a body being determined by positions of nodes of optimal boundary elements (OBE). Computer implementation of the iterative algorithm will be discussed in a future paper. Because of the OBE technique, in a whole shape synthesis process, it is only necessary to operate on values defined on a varying surface of a body. This effective feature of the method may be regarded as its main advantage.

Appl. Math. Modelling, 1985, Vol. 9, June 199

Acknowledgement This wotk is part of the Crucial Problem 05.12 coordinated by the Institute of F u n d a m e n t a l Technological Research o f the Polish A c a d e m y o f Sciences.

References 1 Cruse, T. A. and Rizzo, F. J. (eds.), 'Boundary integral equation method', ACME, New York, 1975 2 Brebbia, C. A. and Walker, S. 'Boundary element tecnniques in engineering', Newnes-Butterworths, London, 1980 3 Banerjee, P. K. and Butterfield, R. 'Boundary element methods in engineering science', McGraw-Hill, London, 1981 4 Kelly, D. W., Moriss, A. J., Bartholomew, P. and Stafford, R. O. 'A review of techniques for automated structural design', Comput. Meth. in Appl. Engng. 1977, 12 (2), 219 5 Venkayya, V. B. 'Structural optimization: a review and some recommendations', hzt. J. Numer. Meth. in Engng, 1978, 13 (2), 205

200 Appl. Math. Modelling, 1985, Vol. 9, June

6

Zienkiewicz, O. C. and Campbell, J. S. 'Shape optimization and sequential linear programming', Chapter 7 of 'Optimum structural design' (ed. R. H. Gallagher and O. C. Zienkiewicz), Wiley. New York, 1973 7 Ramarkrishnan, C. V. and Francavilla, A. 'Structural shape optimization using penalty functions', J. Struct. Mech. 1975, 3 (4), 403 8 Futagami, T. 'Boundary element and linear programming method in optimization of partial differential equation system', in 'Boundary element methods' (C. A. Brebbia), Proc. 3rd Int. Seminar, Irvin, California, 1981, 457--471, Springer-Verlag, Berlin, 1981 9 Burczyfiski, T. and Adamczyk, T. 'The application of the boundary element methods to optimal design of shape of the structure', Proc. Conf. Meth. and hTstrumentation o f Computer Aided Design, Warsaw, 1983 (in Polish) 10 Dems, K. and Mr6z, Z. 'Multiparameter structural shape optimization by the finite element method', Int. J. Numer. Meth. in Engng 1978, 13 (2), 247 11 Dems, K. 'Multi-parameter structural shape optimization', Scientific Bulletin, L6d~ Technical University, No. 29, L6d~, 1980 (in Polish)