An integrated approach to structural shape optimization of linearly elastic structures. Part II: Shape definition and adaptivity

Computing Systems in EngineeringVol. 2, No. 1, pp. 41-56, 1991

0956-0521/91 $3.00+ 0.00 Pergamon Press plc

Printed in Great Britain.

A N I N T E G R A T E D APPROACH TO STRUCTURAL SHAPE OPTIMIZATION OF LINEARLY ELASTIC STRUCTURES. PART II: SHAPE DEFINITION A N D ADAPTIVITY E. HINTON,M. (~ZAK~Aand N. V. R. RAO Department of Civil Engineering, University College of Swansea, Swansea, U.K.

(Received 14 February 1991; revised l May 1991) Almtraet--In this paper we initially study how the number of design variables used affects the final optimum shape of the structure when employing two different types of curves to describe the boundary of the structure, i.e. quadratic Bezier and cubic B-spline curves. The advantage of using better shape definition is highlighted with several examples. An adaptive mesh refinement (AMR) procedure using six-node triangular elements is adopted in the structural shape optimization process. The procedure makes use of an h-version adaptive refinement technique based on error estimates determined from either best-guess stress values or residual terms in the governing equation. An example is presented to illustrate the pertbrmance of these error estimators with respect to their convergence, accuracy and cost of computation. Different strategies for the inclusion of AMR procedures in the shape optimization process are also proposed. Anomalies in predicting the optimum shape due to discretization errors are demonstrated using several examples.

l. INTRODUCTION

1.2. Need for adapth~ity in structural shape optimization

In Part I of this paper, we discussed our overall strategy for the development of an integrated approach to structural shape optimization. Here we concentrate on the need to incorporate better shape definition and A M R procedures into the optimization process.

From an engineering viewpoint the most attractive and yet possibly the most dangerous feature of the finite element (FE) method is the fact that it is approximate. In the hands of an experienced and careful analyst it is a very useful means of obtaining an insight into certain problems for which no other analytical tool is appropriate. However, it is a tool which is often misused especially with respect to the discretization errors arising from poor mesh design. In the FE displacement method the basic unknowns, in other words the displacements, are obtained to a higher degree of accuracy than the stresses and their derivatives which are the input parameters to a structural shape optimizer--in some circumstances such quantities exhibit an oscillatory behaviour. The optimization process uses the results of the FE analysis to predict the optimum shape. In structural shape optimization, the original mesh gradually changes as the iterative optimization method proceeds and the structural geometry changes. Untbrtunately, these changes in the structural geometry can lead to gross element distortions and also to substantial discretization errors. There is even no guarantee that the discretization errors in the original starting mesh are within acceptable limits. Thus, it may be that the resulting optimal structures are based on erroneous FE results and are consequently merely optimal illusions. Recent developments indicate that the discretization error and its distribution in a FE analysis may be improved using adaptivity. For example,

1.1. Need for better shape definition In practical structures the presence of kinks or cusps on the boundary of the structure is usually unwanted as they may lead to singularities and stress concentrations. In structural optimization the smoothness of the final optimum shape of the structure depends, to some extent, on the type and degree of curve used in the shape definition. For example, quadratic Bezier curves, consisting of piecewise parabolic segments (with continuous tangents at the junctions of the segments) are satisfactory when the requirement of smoothness is not strict. Also they are advantageous for decreasing oscillations in the curves. However, in most practical applications there is a need for smoothness of the structural shape. By smoothness of a curve we imply that there should be continuity of curvature. In such applications use of the cubic B-spline curves is recommended as they provide a greater degree of shape control of the structure and in any integrated approach it is necessary to make use of the best available computer aided design (CAD) tools to allow the engineer to concentrate on the more creative aspects of the structural design. 41

42

E. HINTONet al.

Zienkiewicz and Zhu i have introduced a simple error estimator for the h-version of A M R analysis. The error estimator uses a smoothing or projection technique which is based on superconvergent FE stresses. Various authors 2'3 have developed error estimators based on stress jumps around the element boundary and the residual terms in the governing equilibrium equation. The use of the A M R analysis concept in structural shape optimization has been suggested by various researchers 4-7 to increase the accuracy of the analysis as the optimal design is approached and certain constraints become more critical. If the solution is far from the optimum design the A M R may be used to reduce the cost of computation by reducing the total number of degrees of freedom. 1.3. Objectives In this paper we describe the progress of our current research in achieving the main objective which is to combine A M R procedures with FE structural shape optimization to ensure that the FE discretization errors are minimized during the shape optimization process and to obtain efficient and reliable optimum solutions. In the preliminary section we study how the number of design variables affects the final optimum shape of the structure when using quadratic Bezier or cubic B-spline curves. The next section provides an introduction to the concept and basic algorithm of the A M R procedure. In the following sections, the main ingredients of the A M R procedure are discussed including the error estimators, refinement procedure and automatic remeshing technique. An example involving an L-shaped plate is presented to demonstrate the performance of different error estimators and smoothing procedures. The penultimate section reviews the needs of an integrated shape optimization procedure and describes the different adaptive strategies which are designed to yield a specific accuracy in the FE computation. Finally, examples are presented to illustrate the effectiveness of the proposed schemes.

2.1. Plate with a circular hole

The material properties and loading conditions of this problem are the same as in Part I. The objective is to minimize the weight of the plate with an allowable maximum equivalent stress equal to 7.0 mm -2. For comparison we consider only two types of curve: the quadratic Bezier and cubic Bspline curves.t Figure 1 shows the optimum shapes obtained for the different types of curves. This figure also highlights the effect on the final optimum structural shape of using an increasing number of design variables. Figure 1 illustrates the fact that the use of the cubic B-spline curve for the boundary definition gives smoother shapes compared with the quadratic Bezier curve. Also, in spite of imposing several geometric constraints, we find cusps and kinks with the quadratic Bezier curve which do not usually occur with cubic B-splines. This is more evident in the three and four design variable case. One way to avoid the problem of cusps and kinks is to increase the number of design variables, but this may lead to additional problems as it becomes difficult for the Bezier curve to maintain continuity of the tangent at the junctions of the parabolic segments (i.e. the location of the design variables) and the curve breaks down. For example, the Bezier curve was not able to handle the six design variables case. Also, by increasing the Quadratic Bezier curve

Cubic B-spline curve

N D V A B ---- 3

N D V A B ----- 4

N D V A B ----- 5 2. S H A P E

DEFINITION

In this section the effect of using different degrees and types of curve for shape definition is demonstrated using two different examples: a plate with a circular hole and the tall beam considered earlier in Part I. In addition, the effect on the final optimum shape of the structure of the number of design variables used in the problem is also studied for the plate with a circular hole.

N D V A B ---- O

NDVAB

t The quadratic Bezier curve was earlier implemented in the mesh generator developed by Peraire et all* The authors have modified this mesh generator to include cubic B-spline curves.

(Number of

Delign VAriaBles)

Fig. 1. Comparison of the effect on the shape optimization of adopting a shape definition using quadratic Bezier and cubic B-spline curves with different numbers of design variables.

Integrated approach to structural shape optimization--II number of design variables, the cost of the sensitivity calculation by the implicit differentiation approach (IDA) increases. The use of B-spline curves requires the specification of the tangents at the two end points of the curve. With the help of careful specification of these tangents we can achieve better shape control. Alternatively, these tangents can be defined as design variables or the value of the tangents can be linked to the design variables, i.e. as the position of the design variables change the tangent values at the end segments also change as a factor of the arc length between the location of the two end design variables.

43

3. A M R P R O C E D U R E

3.1. Main concept It is well known that, in order to obtain a certain level of accuracy in a FE approximation, the mesh size should be locally refined in areas of stress concentrations and steep stress gradients. To obtain a desired accuracy with a uniform mesh refinement scheme it is often necessary to adopt a very fine mesh in the whole domain. This method is inefficient in most engineering problems; however recently developed AMR procedures produce solutions with minimum cost. To achieve a prescribed accuracy in an optimal way in a FE computation, we must design a mesh in which the error is nearly uniformly distributed.

2.2. Tall beam example 3.2. Basic algorithm

This problem was solved earlier in Part I, using quadratic Bezier curves. The optimum shape of the structure obtained using quadratic Bezier curves is shown in Fig. 2(a). Kinks can be observed on this curve. The final volume of the structure obtained is 3.60m 3 with a maximum stress violation of 5.5%. This problem was solved again using cubic Bsplines without imposing any geometric constraints to prevent cusps or kinks. Considerable differences are observed in the final optimum shape. The final volume is 3.53 m 3 with negligible stress violation. The optimum shape of this structure is shown in Fig. 2(b).

The AMR procedure is as follows: (1) Produce a starting mesh and carry out an initial FE analysis. (2) Based on the results of the FE analysis, evaluate the error estimate. (3) If the error is acceptable then the adaptive analysis is complete; otherwise continue. (4) Re-mesh the whole domain based on a new design mesh density evaluated using the error distribution. (5) Perform the FE analysis again based on the new mesh and go to step 2.

2.3. Discussion

Thus an adaptive refinement scheme is composed of three main ingredients:

Based on our experience we recommend that for practical applications it is better to use cubic Bsplines for boundary definition of the structure for the following reasons:

• a local and a global error estimator; • a mesh refinement (or derefinement) method; and • an automatic mesh generator.

• smoother shapes are always obtained with even fewer design variables; • there is better shape control (i.e. by specification of the end tangents); and • any complex shape can easily be represented.

In the following sections, the ingredients of the AMR procedure are presented in detail.

4.1. Governing equations of elasticity In order to describe the error estimators adopted in the current work we briefly review the sources of discretization error. The governing equilibrium equations of elasticity may be expressed as

Cubic B - s p l i n e c u r v e

Quadratic Bezier curve

4. ERROR ESTIMATORS

Lo - b = O

(i)

in the domain fl, subject to the conditions u = fi on the boundary F~, and

I (a)

I Cb)

Fig. 2. Final optimum shape of tall beam example: (a) using quadratic Bezier curve; (b) using cubic B-spline.

t = t on the boundary F~, where o is the "vector" of stresses, u is the vector of displacement fields, L is a matrix of differential operators, b is the vector of body forces, fi is the

44

E. HINTONet 01.

vector of prescribed displacements on boundary Fu and ~ is the vector of prescribed tractions on boundary F,. The strains ~ may than be calculated from the expression = Lru (2)

inter-element tractions and r are the residual forces and can be obtained from Eq. (7). A particular expression derived in Ref. 2 for two-dimensional problems defines an element contribution as

X~l/2

h2

,,ell,=(~-4--K~errrd.Q+2~f, JrJd,) ,

(9,

and the stresses a can be evaluated as a = De,

(3)

where h is the element size and K is dependent on the problem being solved: for plane stress

where D is the elastic modulus matrix.

K=

4.2. Finite e l e m e n t a p p r o x i m a t i o n In a FE representation, the displacement field [] may be approximated by [] -~ fi = N d ,

(4)

where fi is the FE approximation to u, N is the matrix of shape functions and d is the vector of nodal displacements. The FE strain representation is then given as = LrNd ~ lid

(5)

and the stresses are then obtained from the expression = Dlkl.

(6)

To check for satisfaction of the governing equilibrium equations we substitute # given by Eq. (6) into Eq. (1) and obtain L# - b =

r ~0,

(7)

where r is a vector of residual body forces. 4.3. R e s i d u a l b a s e d error e s t i m a t o r The FE solution gives an approximate displacement field which is continuous over the domain. However, the derivatives of this field are discontinuous across element boundaries. These discontinuities in the displacement derivatives imply strain and hence stress jumps across element boundaries. Various error estimators use the inter-element traction jump around the element boundary and the residual terms in the governing equation over the interior of the element to obtain an error estimate. Techniques using energy norms derived by several authors 2'3 have the general form

,Jell:(C~;rrrd~+C2JI" 3r Jdl)\1/2 , (8) where f~ is the total domain, C1 and C 2 a r e constants, I is the total interface between elements, J are the t It is assumed that appropriate stresses are discontinuous across boundaries between two different materials.

E

(10)

(1 --v)

and for plane strain problems K

E (l + v)(l -

2v)'

(1 1)

where E is the Young's modulus and v is Poisson's ratio. The interior residual term in Eq. (8) is the dominant error term when shape functions consist of piecewise biquadratic functions. However, the dominent error term in Eq. (8) is the boundary term when the shape functions consist of piecewise linear functions. Thus the boundary term can be neglected when estimating errors of quadratic FE approximations. 3 Furthermore, Babuska and Yu 8 have recently shown more generally that for odd-degree elements (elements with odd-degree polynomial shape functions), the inter-element traction jumps dominate the error, but for even degree elements, the interior residuals will be dominant. The error estimation scheme in this section uses these results to estimate the error energy Jl e II of piecewise quadratic approximations in meshes of six-noded triangular elements. The approach taken is to calculate only the interior residual portion of the error in the elements by neglecting the boundary term in Eq. (8). Another residual error estimator which makes use of r using a solution based on hierarchical bubble functions has been described by Baehmann et al. 9 and will be discussed later. 4.4. Z i e n k i e w i c z - Z h u error e s t i m a t o r The error ~ in the solution can be defined as the difference between the exact solution [] and the FE approximation tl. Even though the exact solution is generally not known it is possible to estimate it using the FE results. Recently, Zienkiewicz and Zhu I have developed an error estimator which makes use of projection techniques to define a continuous stress field over the domain.t The basic ideas in the derivation of this error estimator are: (1) The use of a continuous stress field which eliminates the stress jump around the element boundary.

Integrated approach to structural shape optimization--II (2) The use of specific "superconvergent" points in the element for the stress calculation where predicted FE results are more accurate than the other points. Consider the strain energy norm which may be expressed as

The local error for the approximate stresses O is given by ~

= a -

O.

(13)

An error energy norm can be written as lie tl = ( I ( o -

O ) r D - ~ ( o - #) df~) ~/2.

(t4)

As we do not have access to the exact solutions for stresses a, we substitute our best guesses a*, so that ~r ~ a * .

4.5

(t5)

Best-guessstresses

Several simple methods have been proposed in the past and are widely used in practice. These are nodal averaging, least-squares smoothing and Loubignac iteration, all of which have been considered in the present study. The component of the continuous stress field, which is to be found, is approximated ast

a* =~ N, 5,,

(16)

i

where N,, is a C o continuous shape function and #i is the smoothed nodal stress associated with node i. The continuous stress a* is determined by forcing the Galerkin-like condition

45

where typical components may be written as

so= faN~Njd~

and

qi= f N~6 d~.

(19)

For local least-squares smoothing, Eq. (18) is solved by considering individual element domains separately. To reduce the solution cost in global least-squares smoothing, the stress smoothing matrix S can be diagonalized. This reduces the effort required to solve Eq. (18) to one division operation for each stress component at each node. Implementation of this method is very easy--for example the terms of a row may be simply summed to the diagonal. However, if the S matrix for the six noded triangular element is lumped in this manner, the terms representing the corner nodes will be negative. Thus the smoothing matrix loses its positive definite character. Hinton et al.11have proposed a lumping scheme which yields a positive definite smoothing matrix. The total area of the element is divided among the nodes in proportion to the diagonal terms of the consistent S matrix and the shape functions are assumed to have unit value within this region and zero elsewhere. The details of this procedure, as well as other smoothing techniques and error estimators are discussed in Ref. 12. 4.6.

Loubignacalgorithm

Although the use of continuous stress fields eliminates the traction jumps between elements, it does not guarantee a better satisfaction of overall equilibrium. In fact, the continuous stress fields obtained using least-squares smoothing or simple nodal averaging do not satisfy the overall equilibrium conditions. The governing equations of static equilibrium for a domain f~ with a system of body forces and applied surface tractions can be simply established using the Principle of Virtual Work in which

fn&rad~-fnfurbd~-I

5urtdF =0,

(20)

J F~

faN.,

(a* - d) dr2 = O.

(17)

This is equivalent to the least-squares smoothing procedure of Hinton and Cambell. ~° One method which can be employed to obtain a continuous stress field using Eqs (16) and (17) is to assume that the stress tr* is interpolated by the same shape function N i as the individual displacement components. Note that other shape functions (e.g. functions which are one order less than Ni) may also be used. This yields Sd - q = O,

(18)

1"We deal with each stress component (e.g. a x, ay, zxy)in turn.

where virtual displacements 6 u and associated virtual strains & may be represented by the expressions 5u=~Nfd

and

&=~BSd.

(21)

The governing FE equations are obtained by substituting Eq. (21) into Eq. (20)

6d(foBradQ-fuNrbd~-Ijr, NrtdF~=O/

(22)

and as Eq. (22) is true for any 6d, we obtain rer d ~

1

f

O.

(23)

E. HINTON et al.

46

However, when the continuous stress field a* is substituted into Eq. (23), residual nodal forces q are obtained

resent the error in the displacements e using hierarchical (quartic) bubble functions so that 3

~

e = Z ~JgJ B r a * df~ - f = q .

These residual forces can be minimized by the application of Loubignac's iteration scheme. ~3.14The steps to be followed to implement the proposed algorithm are: Solve the problem by the classical method to obtain the displacements do and stresses 60 from Eqs (4) and (6) respectively and set counter i = 0. (2) Increase the count number by 1 (i = i + 1) and calculate smoothed nodal stresses, 0~ starting from d~_ l by global or local smoothing or by nodal averaging, e.g. compute mean nodal values of stresses a.

(31)

j=l

(24)

where l~lj= ~ I 2 and the shape functions ~. are written as g , = ¢~(1 - ¢ - ~)

g2=~2n(1-~-n)

(1)

I

N3 = Cr/:(l - ~ - r/);

gj are unknown nodeless coefficients associated with the distribution of the displacement error e within the element. To obtain these coefficients we solve the stiffness equations for the element which have the form g g = fr,

nec

#i = - - ~ dle-)~•

(32)

(33)

(25)

nse e = 1

where the stiffness matrix g has the form

(3) Interpolate the stresses using the nodal values to give a continuous field ¢* a?=Nd,.

(26)

(4) Evaluate the residual q~ from Eq. (24) corresponding to the current stress field a~ defined in (3). (5) Solve Adi = -- K - l q , .

(27)

(6) Update d and 0

g = fo, [LlqlrDLlq df~

(34)

and 1~1= [lql, 1~12,1KI3].The consistent nodal forces f, associated with the residual error body forces r may be expressed as ? fr = in, l~lrr d.Q

(35)

and the vector of the unknown coefficients g = [glr, gT, g3T]r. Remarks

di =

d~_ 1+

Adi

$i = DBd~.

(28) (29)

(7) Evaluate q, and check convergence II q, II

--<6;

Ilfll

(30)

go to step (2) if convergence is not obtained.

(1) Since the shape functions are quartic, a 13-point integration rule is needed to evaluate g. (2) For the quadratic six-node element the residual body forces are uniformly distributed over each element. (3) The strain energy error for element i associated with the quartic displacements is obtained for the expression 2 tl e I1~, = ½gTgg.

(36)

The convergence is generally fast and 2-6 iterations are normally used.

(4) The total error energy which is used to evaluate the global percentage error is obtained from

4.7. Error estimation based on hierarchical bubble functions

IIe I1~= ~ [1e 112,.

(37)

i=l

In this scheme we first evaluate the interior residual body forces r obtained when the FE stresses are substituted into the governing differential Eq. (7). For the six-noded triangular element we assume that there exists a better solution based on quartic variations of the displacements over the element. Thus we rep-

5. REFINEMENT P R O C E D U R E

5.1. Estimation o f global error The essential problem to be considered in this section is the identification of the adaptive principles

Integrated approach to structural shape optimization--II on which the decision for the refinement will be based. The main principle in h-refinement is to refine the mesh so as to equally distribute the error within each element and to reduce the total error to an acceptable user-specified level. The error energy norms presented in the previous sections are a measure of the absolute value of the predicted error over the domain. The predicted value of the percentage global error r/ is given by the expression

,~ =

II e II If w II

(38)

in which the strain energy of the exact solution is given in Eq. (12). In Ref. 1 the exact strain energy norm of the solution is approximated by

II w II ~ II w* II,

(39)

II w * 11 = (11 ~ II2 + IIe 112)~/2

(40)

where

47

where nej is the total number of elements in the domain under investigation. As the error is evaluated at the element level, we define the parameter ~ for each element as (~ -

lie II, ,

(46)

where ¢~ is called a "refinement indicator". Depending on the values of ~ we can identify three possible element states: • if cg = 1 we have the optimal element size; • if ~i > 1 refinement is necessary; and • if ~ < 1 derefinement is possible. 5.3.

Evaluation of the mesh density

The characteristic feature of the A M R procedure is the use of the current solution to predict the new mesh size/~i. For instance, if the current element size is h~ and the rate of convergence of the adopted element is O(h 1) then we can design the new element size to be h~ which is given by the expression

in which [[ff II is the FE approximation to the strain energy norm and is defined as

~ _ hi

I1~' 1[ =(~adrD-ld d.Q) 1/2.

The values of 1 depend on the smoothness of the solution and on the norm used to evaluate the error. Generally l is taken to be equal to the polynomial degree of the FE approximation. However, near a singularity it has been shown that 0.5 ~< l < 1 where the value of ! represents the strength of the singularity. Equation (47) can be used to evaluate the design mesh density for an automatic mesh generator.

(41)

Alternatively Sibai and Hinton 15approximated [[w* II using the 'smoothed' continuous stress field

II w* II =

(fo

[ a * l r D - l a * df~

(42)

since the strain energy of the FE solution IIif' II2 can be upper or lower bounded by the exact strain energy. The estimation given by Eq. (40) is not valid if II ~ II2 is lower bounded. Equation (38) allows an effective adaptive process to be developed with the principle objective of achieving a specified overall percentage accuracy q - - s a y 5% in many engineering applicationsl6~in the energy norm, i.e. r/~< r7

(43)

or noting Eq. (38), as well 14e II ~< ~ II w* II. 5.2.

(44)

Refinement indicators

As the aim is to obtain a uniform error distribution for all elements, the permissible error for each element is determined by ~6

q 11w* II

n~t2 ,

(45)

i - ~/l/r

5.4.

(47)

Effectivity index

The reliability of the error estimators and various smoothing methods is measured in terms of the effectivity index in the energy norm which is defined by 0 =

ble II

II ~/["

(48)

The error estimator is called asymptotically correct if 0 converges to unity when the errors converge to zero. 17 In practice, we require that IIe II must be close to the actual error II ~ II when the accuracy of the FE solution is in the range of the prescribed values, i.e. when q is close to q. The exact solution may be obtained using analytical methods or approximated by a fine mesh FE computation. 6. REMESHING TECHNIQUES

Another main ingredient of the A M R procedure is automatic mesh generation and remeshing. Remeshing may involve either refining or coarsening the

E. HINTONet al.

48

\

\

very coarse background mesh is always used for the initial mesh generation. The linear triangular background mesh must completely cover the domain to be triangulated. (b ) Specification of mesh parameters. The geometrical characteristics of the mesh are locally defined by means of mesh parameters which represent the size and shape of the element. Figure 4 shows the definition of the mesh parameters for a triangular element: 3 is the node spacing which defines the element size, s is the stretching parameter (if s = 1 we get a quasi-equilateral element), and ~q and ct2 define the local orientation of the mesh.t

Fig. 3. Initial background mesh over an L-shaped domain. 6.2. Mesh generation in A D O P T previous mesh. The following methods can be used to generate new meshes: • Remeshing of the whole domain (AMR). • Splitting/uniting elements (mesh enrichment, MER). • Utilization of a finer uniform mesh in the whole domain (uniform mesh refinement, UMR). 6.1. Mesh generator The A M R procedure requires a convenient means of prescribing the mesh density over the domain, to generate the so called optimal mesh; however, this type of feature is not commonly available in most mesh generators. Based on a recent study, 19in which many mesh generators were compared, the advancing front method of Peraire et all s appears to be one of the best approaches for mesh generation for problems involving adaptive analysis since it incorporates a remeshing facility to allow the possibility of refinement (or derefinement) coupled with directional refinement and allows a significant variation of mesh spacing throughout the region of interest. In the present work, the remeshing of the whole domain is performed using a mesh generator developed by Peraire et al. TM which is ideal for AMR procedures. The basic steps required for mesh generation are (a) the specification of a background grid and (b) the specification of mesh parameters. (a) Specification of a background mesh. The purpose of the background mesh is to control the spatial distribution of element sizes throughout the domain. The mesh parameters such as mesh density, mesh size, etc. are specified at the nodes of the background mesh. From these parameters, the element sizes inside the domain are obtained by linear interpolation. Figure 3 shows a domain to be triangulated and its associated background mesh. The generation of the FE mesh is constructed on the background mesh; a

When carrying out an initial FE analysis in the shape optimization process, the background mesh and the mesh parameters must be specified to generate FE meshes. Depending on whether adaptivity is requested or not we define the values of the mesh parameters to obtain a reasonable mesh. (a) Meshing without adaptivity. In Part I, where adaptivity was not used, the background mesh consisted of only a few triangles covering the domain. Also, since the background mesh remained constant throughout the optimization process, the size of the background mesh was selected anticipating possible shape changes, in such a way that the background mesh always covered the whole domain to be gridded. The mesh parameter 3 is usually selected at each node of the background grid (by user experience) in such a way that a reasonable mesh is obtained. The use of a constant 3 value at all nodes produces similar sized elements at each iteration stage, i.e. a nearly uniform mesh. The number of elements decreases when the volume of the domain is reduced. Figure 5 shows the background mesh and the corresponding generated mesh in the initial and final iterations when no adaptivity is used in the optimization process. (b ) Meshing with adaptivity. When using adaptivity in the shape optimization process the initial background mesh is the same as before but the mesh parameters are chosen in such a way that a very

R2 k., ....

~"s, cq and et2are only used when stretching is required and were included in the original advancing front program by Peraire et aLTMin computational fluid dynamics applications. We have taken s = 1, ctt = 1 and ~2= 0 in all cases.

l X2

I j.....,~.. \ ./ "'w {X ---XI Fig. 4. Definition of the mesh parameters.

Integrated approach to structural shape optimization--II Background mesh

Generated F E m e s h

49

v = 0

\

o_

\\ \

j

y

¢."

7x~

Fig. 7. An L-shaped domain problem definition. (and prescribed accuracy) at the nodes of the current mesh using the expression 6~ = -

\

1

~ /~,

(49)

rise i = 1

Fig. 5. Background mesh and generated FE mesh without adaptivity. coarse mesh is obtained. Generally, these mesh parameters are defined intuitively by experience. A judicious choice of mesh parameters in the initial background mesh is necessary for fast convergence of the solution during the A M R procedure. For example, the value of the size parameter 6 could be taken to be approximately b/5 (depending on the element used) where b is the largest dimension (e.g. side length) in the domain. After the initial FE analysis new values of ~ are calculated by the use of an error estimator Generated F E m e s h

B a c k g r o u n d mesh

where /~ is the new element size determined from adaptive analysis using Eq. (47) and nse is the total number of elements surrounding the node k. There are several types of strategy for incorporating adaptive meshing in structural shape optimiza t i o n - s e e Sec. 8. For example, in strategy (4), at a given stage of the analysis, the mesh used in the previous analysis becomes the background mesh. Thus, in this process the background mesh changes continually during the adaptive analysis at each optimization stage. The number of elements in the background mesh is always changing. In the optimization process without the use of adaptivity the background mesh remains constant. Figure 6 illustrates the changing background mesh in a typical iteration of the optimization process when adaptivity is used. 7. A D A P T I V E E X A M P L E A N D D I S C U S S I O N

7.1. L-shaped plate under edge pressure load

\ \ \

\

To check that the AMR procedure is working correctly for plane stress problems we consider an L-shaped plate under edge pressure load 2° as shown in Fig. 7. The following properties and dimensions are used: elastic modulus E = 100,000.0, Poisson's ratio v = 0.3, thickness t = 1.0, side length L = 100.0 and edge pressure intensity of p = 1.0.

/ /

7.2. Comparison of error estimators The initial mesh of an L-shaped domain shown in Fig. 8, has 20 elements, 103 degrees of freedom and /

Fig. 6. Background mesh and generated FE mesh in a typical adaptivity iteration.

Fig. 8. Initial mesh for L-shaped domain for use with different error estimators.

E. HINTONet al.

50

Table 1. L-shaped plate subjected to uniformly distributed edge loading Method ZZ (NA) ZZ (L) ZZ (G) ZZ (NA/Lg) R (H) R (I)

DOFf

IIff II}

O,

O/

~l,

1707 1449 1451 2265 1149 1217

0.31129 0.31128 0.31127 0.31130 0.31123 0.31124

0.81 0.71 0.67 0.99 0.46 0.45

0.96 1.05 1.14 1.10 1.07 1.15

11.2 9.81 9.19 13.7 6.38 6.25

a strain energy II ~ IIff equal to 0.30544. A solution error of r7 = 1% is assumed to compare the different error estimators. Figure 11 shows the effectivity index and consecutive Fig. 12 illustrates the corresponding rate of convergence obtained using the ZienkiewiczZhu error estimator with various smoothing procedures and residual method error estimators. Note that the effectivity index and rate of convergence curves are coincident for the residual methods. Also the rate of convergence curves of the local and global least-squares smoothing procedures are identical. Table 1 summarizes the results for this problem. In this table D O F means number of degrees of freedom

(a)

(c)

0.88 4.45 0.95 5.35 0.90 8.20 0.89 15.6 0.96 4.60 0.96 2.80

the effectivity index, ~/is the global percentage error and subscripts i a n d f r e f e r to the values at the initial and final adaptivity iteration. The values of the strain energy obtained are compared with the exact solution of 0.3112536. 20 We achieve the prescribed percentage error of ~ = 1% for all error estimators and smoothing procedures. • The Zienkiewicz-Zhu error estimator with simple nodal averaging, ZZ(NA), is considered in order to overcome the problem of high C P U time. Simple nodal averaging has a very high effectivity

DOFf

DOFf 1449

DOFf 1451

DOFf (d)

CPU

required, II ~ II2 is the unsmoothed strain energy, 0 is

1707

(b)

~If

2265

0~

Of

~h

rlf

0.81 0.96i 11.2 0.88

0~

Of

~

~f

0.71 1.05 9.81 0.95

8~

Of

Th

rlI

0.67 1.14 9.19 0.90

Oi

Of

~?i

tlf

0.99 1.10 13.7 0.89g

Fig. 9. Final meshes obtained using Zienkiewicz-Zhu error estimator for L-shaped domain: (a) ZZ(NA); (b) ZZ(L); (c) ZZ(G); and (d) ZZ(NA/Lg).

Integrated approach to structural shape optimization--ll

•

•

•

•

index value in the first iteration but has a slow convergence to the desired accuracy. The final mesh is shown in Fig. 9a. The Zienkiewicz-Zhu error estimator with local smoothing, ZZ(L), gives similar results but a much smaller computational cost when it is compared with global smoothing. The final mesh distribution is presented in Fig. 9b. The use of the Zienkiewicz-Zhu error estimator with global least-squares smoothing, ZZ(G), appears to be sufficiently accurate for shape optimization. The final mesh of the adaptive analysis is shown in Fig. 9c. The approximated strain energy is very close to the exact strain energy. Global stress smoothing has slightly better convergence characteristics but a higher computational cost compared with the other smoothing techniques. Simple nodal averaging (or other smoothing methods) with Loubignac iteration, ZZ(NA/Lg), provides a method for reducing all equilibrium errors. The addition of Loubignac iterations improves the values of the effectivity index. However, this method converges to the exact effectivity index value of 1.0 from above and is a relatively expensive procedure. To achieve the prescribed accuracy the number of degrees of freedom is 2265 as compared with only 1149 degrees of freedom required with the R(H) method (see Fig. 9d). Its use increases the consumed CPU time by a factor of 2or3. The residual method with the higher order bubble shape function, R(H), gives very impressive results giving the best error distribution and the optimal mesh with the smallest number of elements. The rate of convergence of this method is higher than those of other methods. We obtain the desired accuracy using a minimum number of degrees of freedom (see Fig. 10a).

(a)

~

~

51

• When we use the residual method, R(I) [see Eq. (9)], the results are almost identical to the other residual method, R(H). The effectivity index grows very fast with increasing degrees of freedom and is only slightly greater than 1 at the final iteration. This method is the cheapest error estimation method (see Fig. 10b). 7.3. Selection o f error estimators for use with A D O P T In the previous sections we described various types of error estimator that can be used and the different smoothing procedures that are available to calculate the best-guess stress values. The performance of these error estimators on the L-shaped domain was also presented. It is important to note that the different error estimators give different convergence characteristics and different degrees of refinement and have different computational costs. As our main aim is to obtain the optimal solutions accurately, we have selected the residual error estimator with the hierarchical shape functions for adaptive shape optimization because it requires fewer degrees of freedom and gives very accurate results at the final computation compared with other error estimators. The effectiveness and reliability of the residual error estimator with the hierarchical shape functions is also demonstrated by ()zakqa et al. ~z Another aspect of the adaptive process is the selection of ~/(the prescribed accuracy). We only need to obtain a near optimal mesh and hence it is not necessary to give extremely stringent values of q. When a very stringent accuracy is prescribed, the number of elements in the final mesh becomes large and also requires more adaptivity iterations to reach that accuracy. An accuracy between 3 and 8% is reasonable depending on the type of problem and application. Also a limit can be imposed oil the

DOFf

O~

Of

~

~f

1149 0.46i1.07 6.38 0.96

Cb)

DOFf

0~

8f

~

~f

1217 0.45 1.15 6.25 0.96

Fig. 10. Final meshes obtained using residual methods for L-shaped domain: (a) R(H); (b) R(I).

52

E. HINTONet al.

number of adaptivity iterations in each optimization iteration• It has been shown that the number of the elements in the initial mesh greatly influences the rate of convergence and number of elements in the final optimum mesh. The following strategies can be adopted regarding the selection of initial mesh: • Always use a reasonable mesh to start• • Use judgement and experience to grade the starting mesh. • Identify points of singularity in the domain.

8. A D A P T I V E

5.0O4.00-

.::..~?\ -...

70o~

::~::.~??-.."--,,

8.1/0-

%~. "\

\

L

........

..~• ZZ(NA)

I

• ,(. ZZ(G)

I ao [ \

~iL)"

.a.- ZZ'(N'A/Lg) I

[ - * - ra)

,,0o.

i #

STRATEGIES

]

i 5

"~L" ...

I ~

~

i 6

Io

\

" ~ : ~ ~. "%'.~

%~b i 7

'..

i 8

log N

At this stage we may pose the question: in the iterative shape optimization process where A M R is necessary, is it necessary to refine the mesh at every stage of the optimization or only once, or is there any other algorithm? Actually the assurance of the accuracy of a given FE analysis solution by the A M R procedure and the design of the structure by FE shape optimization schemes are two different problems which may be dealt with separately. However, shape change in the design can only be predicted after the required degree of accuracy in the FE solution is obtained. Computational efficiency and accuracy in the final design can be achieved by the proper integration of the two procedures into a single process. Thus, true shape changes are predicted as the A M R improvement is carried out. Wh.en considering the integration of shape optimization and adaptive analysis we have also tried to: • reduce the overall cost of the analysis to at least the cost without adaptivity; • ensure a FE solution of desired accuracy at every iteration of the shape optimization, and • allow for facilities to make the whole analysis fully automatic if necessary•

Fig. 12. Convergence curves for different error estimators. updated background mesh strategies. The various steps involved in implementing the different strategies are discussed below. 8.1. C o n s t a n t b a c k g r o u n d mesh strategies

Here, we consider the strategies in which the background mesh remains unchanged at the beginning of the each optimization iteration. S t r a t e g y 1. Use A M R to evaluate the ~ values in the first iteration--keep the 6 values constant throughout the optimization and check the final accuracy• (1) Estimate the J values using A M R in optimization iteration 1.

(a)

The adaptive strategies can be broadly classified into (a) constant background mesh strategies and (b)

feasible solution 120

38~0o ~

I00

(b) ~ .,,o~oo ~

0.8o

o.6o

(9. . . . . . . . . . . er""' z(" ."

×........

•

i

,

i

~

1..~ R(H)

//

•

t

200 400 600

,

i..~.. ZZCNA) t" .s. ZZ(L) I" ,(• ZZ(G) i-.~- ZZ(NAILg)

t ,

i , i

•

i

,

i

,

i

, i

,

i

,

80o i00o 1200 1400 1600 lS0o 2000 2200 2400

DEGREES

OF FREEDOM

Fig. 11. Effectivity index vs degrees of freedom for L-shaped domain.

ITERATION NO.

Fig. 13. Square plate with a square hole: optimization process with strategy 2: (a) boundary shape changes; (b) design history•

Integrated approach to structural shape optimization--ll (2) Keep the background mesh and the 6 values unchanged throughout the optimization iterations. (3) When the optimum solution is reached check the accuracy of the FE solution. (4) If the error estimate is acceptable in the final analysis then the optimization is complete; otherwise continue. (5) Restart the optimization from the final iteration and use adaptivity in subsequent iterations. In this strategy there is no guarantee on the accuracy of the FE computation after the first iteration. When some of the constraints are activated the global percentage error may increase. It may turn out that the optimization will converge to a solution which is not the real optimum. Another disadvantage of this method is that the cost of the computation may increase if the optimization has to restart from the final iteration. Strategy 2. Use A M R to update the 6 values when the constraints are activated. (1) Select the 6 values to give a coarse mesh. (2) Carry out the FE analysis and calculate the constraints. (3) If the constraints are activated then use AMR to improve the mesh; otherwise continue the optimization. In this algorithm the cost of computation is reduced and the accuracy of the FE computation is improved at the critical design points. Strategy 3. Use AMR to update the 6 values in every iteration. (1) Use a coarse mesh and constant background mesh. Mesh

Equivalent stresses

53

(~)

feasible solution solution 38000o0 300o000(b)

' ':" '~ ' ~ . ' , ' , ' o ' , ' , . ' . ' ~ ' , ' o ~ ITERATION NO.

Fig. 15. Square plate with a square hole: optimization process with strategy 3: (a) boundary shape changes; and (b) design history. (2) Carry out the AMR procedure during every iteration. This can yield considerable computational savings in the overall iterative solution. Here, we get high accuracy and fast solution convergence. 8.2. Updated background mesh strategy Here we consider a single strategy in which the background mesh is updated at every optimization and adaptive iteration. Strategy 4. Use AMR to evaluate ~ values in all of the optimization iterations and update the background mesh for the next iteration. (1) Estimate the 6 values using AMR and update background mesh for every optimization iteration. The accuracy of the solution is assured throughout the optimization. However, the cost of the algorithm dramatically increases. Also in this strategy if the structure expands locally beyond the background mesh then the available mesh generator cannot handle this specific case. 9. E X A M P L E S

9.1. Illustrative example

Fig. 14. Square plate with a square hole: initial and final shapes and associated equivalent stress distributions for strategy 2.

In order to assess the efficiency of the various adaptive strategies the problem of the square plate with a square hole is considered. The geometry, material properties and loading conditions of the

E. HZNTONet al.

54

problem are assumed to be similar to the problem described in Part I. To integrate the shape optimization and the A M R procedure we have used constant background mesh strategies with quadratic Bezier boundary curve definition.

.......................................... jSSs:::::::::,:...... . . . . . . . . . . . . . . .

(a)

,-,,-~,.,

.........................

~.. "~:}x

.....

""

'\

'iX

9.2. Results and discussion 40~0.00

Strategy 1. The A M R procedure is used only in the first iteration to evaluate the ~ values and later these 6 values are used for subsequent optimization iterations, i.e. the first FE analysis mesh is kept as a constant background mesh. This strategy fails after the first optimization iteration for the square plate with a square hole problem since one of the boundary points of the second iteration is outside the background mesh. As mentioned before the mesh generator cannot handle the situation when the boundary of the structure is outside the background mesh. Strategy 2. In this strategy the A M R procedure is triggered when one or more stress constraints are activated during optimization. The boundary changes and design history are shown in Fig. 13. Figure 14 presents the initial and final meshes and the corresponding equivalent stress distributions. In the initial iterations, the volume is only slightly reduced because of constraint violation at the singularity points. Local refinement occurs near the singularity points. In the first five optimization iterations, the global percentage error ~/is over 6% and less than 3% in the next six optimization iterations. In the last three iterations ~/again increases to more than 8%. However, for a given prescribed accuracy, the number of elements in each optimization iteration is always less than the number obtained in the solution without adaptivity which is presented in Part I. The volume of the plate is reduced from 40,250.0 to 28,736.0mm 3 in 14 Mesh

Equivalent

stresses

~asiblesolution

3 5 0 ~ aO

solution 3OOOO.O0

g 2OOO0.OO

I~000.00

IC000.~

50~0.00

~ _ ~ . ~ - - - ~

i

2

•

i

4

,

i

,

6 ITERATION

i

i

~

Io

i2

I

•

14

NO.

Fig. 17. Optimization process for connecting rod problem: (a) boundary shape changes; (b) design history. iterations with no violation of the stress constraints in the final iteration. Strategy 3. The same constant background mesh adaptivity strategy used earlier is now employed in each optimization iteration. Figure 15 displays the boundary changes and design history. Inital and final meshes and corresponding equivalent stress distributions are presented in Fig. 16. The final optimum shape is obtained in 17 optimization iterations and the volume has dropped from 40,250.0 to 28,500.0 mm 3. As expected, we get a smaller volume and a better final optimum shape using this method. On average, two adaptivity iterations are performed in each optimization iteration to reduce the global error below 3%. As with the other adaptive strategy, a considerable reduction in the CPU time is achieved. From the results of the various strategies it is seen that strategy 3, which involves adaptivity analysis at every optimization iteration, gives the most satisfying results. This strategy is adopted in the following example. 9.3. Connecting rod example

Fig. 16. Square plate with a square hole: initial and final shapes and associated equivalent stress distributions for strategy 3.

Using the third adaptivity strategy in conjunction with a quadratic Bezier curve boundary definition and a residual error estimator based on the hierarchical bubble shape functions, the connecting rod problem is studied. The geometry, material properties and loading conditions of the problem are assumed to be identical to those given for the problem described in Part I. Figure 17 shows the boundary changes and corresponding design history. The initial and final meshes and the corresponding equivalent stress distribution for a 3% prescribed accuracy are presented in Fig. 18; because of the use of an error estimator and adaptivity strategy, the final optimum weight and

Integrated approach to structural shape optimization--II Mesh

55

E q u i v a l e n t stresses

Fig. 18. Initial and final shapes and associated equivalent stress distributions for connecting rod problem. computation time are less than those given in the solution presented in Part I. The initial volume of 35,473.6 m m 3 reduces to final volume of 3815.4 mm 3 in 13 optimization iterations.

• A b o v e all, adaptive analysis (or at least error estimation) should be an important part of any optimization iteration where stress concentrations or dramatic shape changes are likely to occur.

10. DISCUSSION AND CONCLUSIONS

Acknowledgements--The authors wish to thank Dr J. Peraire for the use of an automatic triangular mesh generator based on the advancing front method.

The advantages of using cubic B-spline curves for the definition of the boundaries of the structure have been demonstrated using several examples. Various error estimators have been discussed and compared in the context of shape optimization. Several strategies have been presented to integrate the shape optimization and adaptive analysis. An A M R procedure using the residual method error estimator based on hierarchical bubble shape functions was applied at every optimization iteration to improve the F E solution. By comparing the results in Part I with those in the previous section we have shown that comparable or even better results can be obtained using the A M R procedure in integrated shape optimization. The method behaves well in the square plate and connecting rod examples and provides the best approximation to the true optimum solution. Furthermore, a degree of confidence in the reliability may be attached to the optimum solution. However, the accuracy of the A M R procedure depends on the error estimation, refinement procedure and mesh generator. The following general conclusions can be made on the various aspects of shape optimization and adaptive analysis in an integrated shape optimization process: • The use of' an automatic mesh generator with a cubic spline boundary definition is recommended. • The residual method with the higher order bubble shape function is the most suitable error estimator for adaptive shape optimization. It gives satisfactory reliable results with low computational cost. • The adaptive strategies 2 and 3 can be used in any integrated adaptive shape optimization problem. Both methods work well and assure the accuracy of the solution.

REFERENCES

1. O. C. Zienkiewicz and J. Z. Zhu, "A simple error estimator and adaptive procedure for practical engineering analysis," International Journal for Numerical Methods in Engineering 24, 337 357 (1987). 2. I. Babuska and W. C. Rheinboldt, "A posteriori error estimates for the finite element method," International Journal for Numerical Methods in Engineering 12, 1597-1615 (1978). 3. D. W. Kelly, J. P. de S. R. Gago and O. C. Zienkiewicz, "A posteriori error analysis and adaptive processes in the finite element method: part I--error analysis," International Journal for Numerical Methods in Engineering 19, 1593 1619 (1983). 4. M. H. Imam, "Three dimensional shape optimisation," International Journal for Numerical Methods in Engineering 18, 661~73 (1982). 5. M. E. Botkin and J. A. Bennett, "The application of the adaptive mesh refinement to shape optimisation of plate structures," Accuracy Estimates and Adaptive Refinements in Finite Element Computations (edited by I. Babuska, O. C. Zienkiewicz, J. Gago and E. R. de A. Otiveira), Ch. 3, John Wiley, Chichester, 1986. 6. M. S. Shephard and M. A. Yerry, "'Automatic finite element modelling for use with three dimensional shape optimisation," Discretisation Methods and Structural Optimisations--Procedure and Applications (edited by H. A. Eschenauer and G. Thierauf), Springer, Berlin, 1989. 7. N. Kikuchi, K. Y. Chung, T. Torigaki and J. E. l~aylor, "Adaptive finite element methods for shape optimisation of linearly elastic structures," Computer Methods Applied Mechanics and Engineering 57, 67 89 (1986). 8. I. Babuska and D. Yu, "Asymptotically exact a posteriori error estimator for biquadratic elements," Technical Report, Institute for Physical Science and Technology, Laboratory for Numerical Analysis, University of Maryland, Technical Note BN-1050, 1986. 9. P. L. Baehmann, M. S. Shephard and E. J. Flaherty, "A posteriori error estimation for triangular and tetrahedral quadratic elements using interior residuals," SCOREC Report, Scientific Computation Research Center, Rensselaer Polytechic Institute, 1990.

56

E. HINTONet al.

10. E. Hinton and J. S. Campbell, "Local and global smoothing of discontinuous finite element functions using a least squares method." International Journal for Numerical Methods in Engineering 8, 461-480 (1974). 11. E. Hinton, T. Rock and O. C. Zienkiewicz, "A note on mass lumping and related processes in the finite element method," International Journal of Earthquate Engineering and Structural Dynamics 4, 245-249 (1976). 12. M. Ozakqa, N. V. R. Rao and E. Hinton, "A comparison of error estimators for triangular elements," Int. Research Report, University College of Swansea, July 1991. 13. G. Loubignac, G. Cantin and G. Touzot, "Continuous stress field in finite element analysis," AIAA Journal 17, 1645-1647 (1977). 14. G. Cantin, G. Loubignac and G. Touzot, "An iterative algorithm to build continuous stress and displacement solutions," International Journal for Numerical Methods in Engineering 12, 1493-1506 (1978). 15. W. A. Sibai and E. Hinton, "Adaptive mesh refinement with the Morley plate element," Proceedings of the

16.

17.

18.

19. 20.

NUMETA 90 Conference, Swansea, Vol. 2, pp. 1044-1057, Elsevier, London, 1990. I. Babuska, O. C. Zienkiewicz, J. Gago and E. R. De A. Oliveira, Accuracy Estimates and Adaptive Refinements in Finite Element Computations, John Wiley, Chichester, 1986. M. Ainsworth, J. Z. Zhu, A. W. Craig and O. C. Zienkiewicz, "Analysis of the Zienkiewicz-Zhu a posteriori error estimator in the finite element method," International Journal for Numerical Methods in Engineering 28, 2161-2174 (1989). J. Peraire, M. Vahdati, K. Morgan and O. C. Zienkiewicz, "Adaptive remeshing for compressible flow computations." J. Computational Physics 72, 449-466 (1987). J. Sienz, "Finite element mesh design with adaptive procedures," M.Sc. thesis, Department of Civil Engineering, University College of Swansea, 1990. J. Z. Zhu, "Error estimation adaptivity and multigrid techniques in the finite element method," Ph.D. thesis, Department of Civil Engineering, University College of Swansea, 1987.