Mesh generation with adaptive finite element analysis

Mesh generation with adaptive finite element analysis E. Hinton, N.V.R. ILao and M. Ozakga Department of Civil Engineering, University College of Swansea

This paper deals with key aspects of the coupling of a u t o m a t i c mesh generation with adaptive analysis and shape optimisation of planar and surface structures. Adaptive analysis requires the design of nearly optimal meshes with varying element sizes to achieve a prescribed accuracy. In any shape optimisation process as the region to be discretised changes continuously and sometimes dreustically, mesh generation should ideally take place without intervention by the analyst and without any excessive distortion of elements. To carry out these processes in a fully automatic and efficient manner a highly robust, versatile and flexible mesh generator is required. A brief introduction to the features of one such automatic mesh generator, based on the advancing front method, is described. Several examples are presented exploiting the potential of this mesh generator for adaptive analysis and shape optimisation. The integration of the rnesh generator with these procedures is discussed as well as some issues which required special attention.

Key Words: finite elements, adaptive analysis, error estimation, automatic mesh generation, advancing front method, structural shape optimisation.

1. I N T R O D U C T I O N

I. I Preamble The finite element (FE) method is now firmly established as an engineering tool of wide applicabilil, y. When carefully used by experienced and diligent engineers, FE based models can provide much valuable insight into structural behaviour. However, designers who do not possess extensive expertise in numerical analysis, may unquestioningly believe that all FE results are highly accurate or even exact. This has led to recent emphasis on adaptive analysis or adaptive mesh refinement (AMR) procedures to improve the reliability of the FE method and to ensure that the results produced are of appropriate accuracy.

1.2 Automatic mesh generation in AMR procedures A*IR procedures involve the design of nearly optimal meshes with varying element sizes to achieve a prescribed accuracy. The importance of AMR procedures in industrial applications has led to increased research on fully automatic mesh generators which require only the specification of the boundary and mesh size distribution over the domain. The success of AMR proce-

dures depends to a large extent on the efficient coupling between the adaptive FE analysis and automatic mesh generation. Tile 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. A method for estimating the new element sizes is required in order to guide the generation of an entirely new mesh to achieve a prescribed accuracy. The new element sizes can be evaluated using an error estimate from the results of a previous FE analysis and there are different ways of estimating this error based on residual or projection methods. Once information concerning the new element sizes has been estimated it has to be linked to the mesh generator to perform the desired discretisation by transferring the new element sizes to tile nodes of the current mesh or to the boundary of the domain to be discretised. The extension of this procedure for surface structures further complicates the issue, as mesh generation and error estimation on curved surfaces are required coupled with accurate estimation of mesh parameters such as element size. Another aspect which requires attention in this case is the mapping of elements and

@1991 Elsevier Science Publishers Ltd 238

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined

Mesh generation with adaptive finite element analysis: E. Hinton et al.

mesh parameters onto the curved surface. The mesh generation algorithm should be efficient and fully automatic.

1.3 Automatic mesh generation in shape optimisation Shape optimisation is a further area where proper coupling between the mesh generation and mesh optimisation procedures is an essential requirement. In any shape optimisation procedure as the region to be discretised changes continuously and sometimes drastically, mesh generation should ideally take place without intervention by the analyst and without excessive distortion of the elements during the whole process. Recent efforts to integrate shape optimisation and AMR procedures have highlighted the need to achieve proper coupling with mesh generators if the whole process is to be carried out in a fully automatic and efficient manner. Also the way in which the design variables are selected in shape optimisation procedures depends on the input data required by the mesh generators.

It should incorporate a convenient geometric representation of boundaries and be able to represent complex shapes easily. It should be able to automatically generate meshes of different sizes over general domains as specified by the user. It should possess a convenient means of prescribing the element size variation over tile domain. • The input data to the mesh generator should be minimised. o The mesh generator should be portable, so that it can be attached as a part of the analysis module.

I. 4 Objectives and layout

• It should be flexible, in such a way that its potential can be exploited in other applications (e.g. mesh generation on curved surfaces and with shape optimisation).

In this paper the main emphasis will be on the efficient coupling of a mesh generator with AMR and shape optimisation procedures. The following section deals with the description of key features of a mesh generator based on the advancing front method (AFM) which is used in this study and had originally been developed for AMR processes in fluid dynamics applications'. Its extension to mesh generation on curved surfaces is then described. The next section presents an overview of current AMR procedures, error estimation and the calculation of new element sizes for the automatic mesh generator. Several examples illustrate the effectiveness of proper coupling between AMR and AF,\I procedures. Some issues are raised concerning adaptive FE analysis of surface structures (such as element size estimation and degeneracies in the shape of elements) and strategies to address these issues are discussed. The penultimate section discusses the coupling of automatic mesh generation with shape optimisation procedures and some examples are then presented. Finally, some conclusions are listed.

2. A U T O M A T I C

2.2 The advancing front method (AFI~D The above requirements of the mesh generator do not seem to be fully satisfied by generation using structured grids. In recent years a wide variety of algorithms have been devised for the generation of unstructured grids around complex geometrical shapes. A survey of different approaches for unstructured grid generation has been presented by Shepherd 2. Based on a recent study 3 the AFM ofPeraire et all 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. Another important feature is that elements and points are generated simultaneously as compared to other algorithms which require the generation of all the interior nodes in the domain before their connectivities can be defined. A brief description of the algorithmic steps involved in this procedure is now presented. For more details interested readers can refer elsewhere 1,4.

MESH GENERATION

2.1 General requirements

2.2.1 Selling up of background mesh

Due to the many challenging applications of the FE method with the ever increasing emphasis on improving the reliability of tile method, the requirements of the mesh generator are quite substantial. In general, for the automatic mesh generator to be robust, versatile and efficient, it should have the following features:

In the AFM it is necessary to specify a background mesh the main purpose of which is to control ttle spatial distribution of element sizes or mesh density through-

1The original AFM mesh generator x adopted here used

quadratic Bezier curves. It was later modified to include cubic B-spline curves.

out the domain. Tile background mesh consists of triangles with corner nodes at which the values of various mesh parameters (such as mesh size, stretching, and orientation of the mesh) are specified. The background mesh completely covers the whole domain to be discre-

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined 239

.'+f~'-h 2c,~eration u'ith adaptit, e fiTzite ele~ent ,It~as,~js~:,.- E. ;qis~c,y..:~:-cal.

background \\

mesh

generated

FE

mesh

-\

\

\

\

\\ \

\ \,

\

\

\

\

Fig. l Specification of background mesh for generating uniform and graded mesh.

j j

....\, ....

tised. Generally, one or two elements will be sufficient if only a uniform mesh density is required. Typical background meshes for uniform and graded refinement are illustrated in Fig. 1. Definitions of the mesh parameters specified on the background mesh are presented in Fig. 2: ~ is the nodal spacing which defines the element size, s is the strelching parameter and oil and a.~ define the local orientation of tile mesh 1.

2.2.~ Representation of the domain boundaries X

CCZ t~

I

\\\.\L...... ,j" ~X

1

Fig. 2 Definilion of mesh paramelers.

240

Adv. Eng. Software, 199i, Vol. 13, No. 5/6 combined

Tile shapes of structural domains encountered in practice are so arbitrary and complex that it is essential tt~at they should be represented in a convenient way using CAD tools such ms Bezier or B-spline curves. In most structural applications there is a need for smoothness (i.e.. continuity of curvature) of the structural shape. In such applications the use of cubic Bspline curves is recommended t because it allows a flex-

Mesh ,],3n.eralion with adaptive finite element analysis: E. Hinton et al.

9

® 13

-,

12 6

7

©

18

®

5

19 ¢

([~)

3

1¢

15

0 2

(b Orientation

o

Key points (D

Segments

Fig. 3 Representation of the boundaries of the structural domain showing segments, key points and orientation. ible system of both interpolation and interactive manipulation to co-exist. To avoid confusion we define certain terms which will be used frequently in this paper. The boundary of a typical structure is shown in Fig. 3 and is formed by an assembly of segments. Each segment is a cubic Bspline curve passing through certain key points. Some key points are common to different segments at their points of intersection. To represent a straight edge the user has to provide a minimum of two key points lying on the boundary of the segment and three points to represent a curved edger. The orientation of the curve is shown by arrows. In structural domains with openings, the segments of the interior boundary are oriented in a clockwise manner whereas the segments of the outer boundary are oriented in an anticlockwise manner. By orienting the segments in this way the region to be triangulated always lies to the left of the segment 1.

2.2.3 Discretisation of boundary segments Based on the information stored in tile background grid, tile boundary of the structural domain is discretised into a series of lines which will form the bases of triangles adjacent to the boundary. The nodes on the boundary are spaced according to the values of t To represent a cubic spline uniquely two end points and the tangent values at those points are sufficient, but to ensure that the cubic spline passes through certain desired points we use a minimum of three points.

6(a) representing the variation of the local values of the spacing along the length of the curve, a being the arc length p a r a m e t e r of the curve. The 6(a) values on the boundary are interpolated from the background mesh. The number of sides N,a into which the boundary is divided is computed (to the nearest integer value) by the expression

A,d = /oL ~-~a)da

(1)

where L is the length of the boundary edge and the positions of the node along the curve are determined from 4 ¢(a) = ~

da

(2)

2.2.4 Generation of initial front and triangles As the AFM proceeds the initial front is formed by a list of oriented elements of the discretised boundary of the region to be meshed. Tile orientation of an element in the front is defined by the order in which its nodal points are listed. From this front an item is selected which is used as a base for a new element in the mesh. During the generation of tile triangles, any item from the front which is available to form an element side is termed active. The new element is generated by creating a new point or selecting one of the nodes of the front. After the

Adv. Eng. Software, 1991, Vol. I3, No. 5/6 combined 241

~.~,z j_.,~. . . . . a¢in,z~.~, w,tt. a~L~pli,_,e j, ntte element ~na[,.jsts: E. ]]~nton c't a[.

12 ;>

11 o

10 © ~

9 ©

8

13 1/,

initial front

15

--'-----0

0

0

;--------0

\

0"

0"'

\ \ / J

18 o \

-h 12

11

10

12

9

11

9

10

element generation

8

o

13

13

,

)

1/-,

18

< 6

15

t5

5

17x~o' 1

2

3

1

2

3

Fig. # Generation of triangular elements in the advancin 9 front technique.

242

Adv. Eng. Software, I991, Vol. 13, No. 5 / 6 combined

updated front

3le.sf~ jeneration with adaptive finite element analysis: E. fIinton et al.

generation of the new element, the generation front is subsequently updated by removing those items from the front which are no longer active. This procedure is illustrated in Fig. 4. Thus in the AFM, the generation front changes continuously and has to be updated whenever a new element is formed. This process is continued until the generation front is empty 4.

In this paper we will concentrate on the use of triangular elements to demonstrate the efficient coupling between adaptive analysis and automatic mesh generation.

2.3 Quadrilateral element mesh generation using AFM

To carry out adaptive FE analysis of surface structures (e.g. shell structures), it is necessary to represent surfaces in such a way that 2D mesh generators can be used to generate meshes on these surfaces. In this paper we have used a parametric representation of surfaces with Coons patches for mesh generation on curved surfaces. The algorithmic steps for mesh generation on curved surfaces are described below.

3. M E S H FACES

Recently Zhu el al 5 have extended the AFM for quadrilateral element mesh generation. The methodology is based on the facts that (i) a region can always be subdivided entirely into quadrilaterals if the polygon which forms the boundary of the region has an even number of sides and (ii) a quadrilateral can be formed from two triangles which share a common side. This procedure has been successfully used by Zhu el al for adaptive FE analysis of planar structures. Figure 5 illustrates the mesh of a machine part after adaptive

P=l

~ . - _ _ _ ~ - - -

.

• "--;

-

~

,.

,

-

>

GENERATION

ON CURVED

3. i Parametric representation of surface structures Most surfaces can be divided into a network of patches. These patches can be constructed solely in terms of the information given on the surface boundary and certain auxiliary functions of u and v, where u and v are parametric coordinates which vary from 0 to 1 along the relevant boundaries of each patch as shown in Fig. 6. The patches so defined are known as Coons patches if interpolation within each patch is carried using the expression z 3

3

i=0

j.=O

=

•

~:'~--~

~

degreesoffreedom: 565 .~_

~

percentage

~

accuracy

(0) = 9 . 7 5 %

SUR-

(3)

o,j,,,¢

where r is the position vector of a point on the surface and alj are the coefficients which contain linear combinations of the position and derivative (tangents and twists) vectors at patch corners.

finalmesh

r(u,7)

y--

....'

u=O

- - " " - - ; A ~ - ' ~ I ~ " degreesof freedom=3155 percentage

accuracy

• r~(u.v )

u=1

I r(1,v)

(r/) = 4 . 8 5

u=O Fig. 5 Quadrilateral mesh obtained from adaptive analysis to achieve prescribed accuracy of 5 %.

analysis has been carried out to achieve a prescribed accuracy of 5% 6.

v=O

r'(l,l,O}

u=1

v=O

Fig. 6 A parametric surface patch.

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined 243

Mesh

y~neration with adaptive ,finite e~emen~ ' ~ analysis:

~P.

Hinton. et~ a!.

3.2 Mesh parameters in the parametric plane

patches the stretching factor s is also used where the value of the stretching factor s is taken as equal to tile aspect ratio of the patches, llaving determined these parameters the element sizes 6 are scaled accordingly in tile parametric plane. In this way the resulting surface discretisation is compatible with the local mesh parameters in the 2D parametric plane.

To maintain the desired element sizes on tile 3D surface, the element sizes in the parametric plane should be determined as accurately as possible. To achieve this for approximately square patches on the surface, we determine the parametric ratio of the patches, which is equal to the ratio between the arc length da oll the surface patch to the paranaetric length dl on the parameter plane in the appropriate direction. For rectangular

"---~

/

/

'

7

y

surface

1

-~x

/ /

r(u,v)

v

parametric plane

J

i

! 1

! i! 3 l

ff

I E

_T

I

1!

1

i i

i !

Fig. 7 Parametric representation of surfaces.

244

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined

U

M e s h ~jer~eration w i t h a d a p t i v e f i n i t e e l e m e n t analysis: E. H i n t o n et al

sponding parametric plane. Knowing the positions of the nodal points on the parametric plane, the corresponding positions of the nodal points on the surface can be determined using (3). The element connectivity on the surface remains identical to the one on the parametric plane as there is a one-to-one correspondence between the surface and its corresponding parameter plane. Figure 8 illustrates the mapping procedure for generating meshes on curved surfaces.

3 . 3 M e s h g e n e r a t i o n on the 2D p a r a m e t r i c p l a n e

Once the mesh parameters in the parametric plane have been determined, mesh generation using the 2D AFM mesh generator is carried out. The parametric plane consists of squares of unit dimensions corresponding to each patch on the surface. Therefore, the number of squares on the parametric plane is always equal to the number of patches on the surface. Figure 7 shows a general surface divided into patches and its corre-

,r-,

,....... ,,--:.:'-v<. . . . 7...... Y ...... ":. . . . ~" i ....~ ...... ~....... '~....... ~-"

~.. :: ......,, " """i -'"

/......r

/. ..... ; ....

//

i

i

!

/ -I- - /

surface

!

,.,

.........!.....

/.......... j...1-". . . . . . . . I

,,,

,'

I

:~Z~ ;

.*°°~

I

l ..........

"--.4 . . . . . . . . . . . . . . . . . . . . . .

',,........

,

L . . . . . . . . . . . . . . . . .

--.1. ............

×

y/ plane

< < < < > >>.,~

k------+

.

.

.

.

.

.

.

.

.

.

.

-I

U L

t

x

Fig. 8 M a p p i n g curved surfaces.

~

~

~

a

procedure f o r g e n e r a t i n g

meshes

on

A d v . E n g . S o f t w a r e , 1991, Vol. 13, N o . 5 / 6 c o m b i n e d

245

}.[esh 9enerutzo~'~ ~rith adapti,'e .finzte eLz.,v~e;z~ (zn,~lgszs: Z. ,:~:ortr~ '

et aL

V

1

C

~..-X

Y

3-O

w

v

A

B

"-U v

v

A

2-13 I0, lxI0,71

tJ'

B 2-0

Fig. 9 Stretching and shearing of unit sq~,are patch to approacimale the surface segment. 3.4 Mesh generation for non-square patches on the surfs ce Usually it is necessary to use non-square or even degenerate patches to represent complicated geometries. In such cases to achieve very accurate discretisation we may use the rever~e mapping technique. In this technique instead of directly generating meshes on the paramctric plane, the meshes are generated in another

plane in which tile squares of the paranaetric plane are stretched and sheared so as to approximate the surface segment in a 2D domain. Then, using bilinear isoparametric mapping the generated points are mapped onto the parametric plane from which they are subsequently mapped onto tile 3D surface. This procedure is illustrated in Fig. 9. This approach h ~ also been used by Lohner et al S.

,./

Fig. I0 Ezamples of mesh generation on curved surfaces

246

Adv. Eng. oCoftware, 1991, Vol. 13, No. 5/6 combined

3.resh ,je~eration with adaptive .finite element analysis: E. Hinton et al.

II

3.5 Examples To illustrate the mesh generation on curved surfaces we consider a few examples as shown in Fig. 10. It is important to mention here that in this paper we have restricted ourselves to rectangular geometries. IIowever, extension to non-rectangular geometries can be achieved easily by using degenerate triangular patches. This technique has also been used to represent the complex geometry of an aeroplane: see Peiro a.

uate this energy norm II e two approaches have been used: one based on the residual forces and the other based on besl guess stress values. In this paper only a brief discussion of these two approaches is given. For more details readers can refer to Ozakca et al 9

5.1 Residual body forces The governing equilibrium equations of elasticity may be expressed as LTer--b = 0

4. INTEGRATION DURES

WITtI

AMR

PROCE-

(4)

in a domain ~, subject to tile conditions u = fi on the boundary Fu, and

Ilaving described the important features of the AFM mesh generator and demonstrated its extension to mesh generation on curved surfaces we now concentrate on the coupling of this mesh generator with AMR procedures. A typical AMR procedure is usually based on the following algorithm: 1. Produce a starting mesh and carry out an initial FE analysis. 2. Based on the results of the FE analysis produce an error estimate. 3. If the error is acceptable stop; otherwise continue. 4. Usi.ng an automatic mesh generator re-mesh based on the information from the error estimator. 5. Carry out a further FE analysis based on the new mesh and go to step 2. As we can see in step 4 we have to give some information to the mesh generator based on an error estimate. This information is given in the form of a distribution of new element sizes or spacing, 6. . . . throughout the domain. The values of the 5new are specified at the nodes of the current mesh, which then becomes the background mesh for the generation of an entirely new mesh. The 5ne~ values can be determined from an error estimate based on the current FE analysis. Thus the type of error estimator used will determine tile degree of refinement of the mesh. In the following sections we briefly deal with different types of error estimators and demonstrate how new element sizes for the mesh generator are calculated.

t = t

on the boundary I',

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, ii is the vector of prescribed displacements on boundary Fu and t is the vector of prescribed tractions on boundary Ft. In a FE representation, the stresses o" m a y be approximated by 5" obtained from the expression 5- = D B d

where D is the elastic modulus matrix, B is the FE strain-displacement matrix and d is the vector of nodal displacements. To check for satisfaction of the governing equilibrium equations at any point, we substitute the FE stresses & given by (5) into (4) and obtain L T& - b = r ¢ 0

5.2 Residual based error estimator The FE solution gives an approximate displacement field which is continuous over the domain. IIowever, the derivatives of this field are discontinuous across element boundaries. These discontinuities in the displacemcnt derivatives imply strain and hence stress j u m p s across element boundaries. The error estimates based on residual methods 1°,11 use the inter-element traction j u m p around the element boundary and the residual terms in the governing equation over the interior of the element to obtain an error estimate. The error energy norm llcll is then given by

=

--

(6)

where r is a vector of residual body forces which is not usually a null vector.

l,ol, @l£r rdn + c2fJTJ ')l/' 5. ERROR ESTIMATION BASED METHODS

(.5)

(r)

RESIDUAL

To give a quantitative measure to the magnitude of error in the FE analysis, error norms are used. The most frequently used norm is the energy norm II~ll. To eval-

where ~ is the total domain, C1 and C2 are constants, I is the total interface between elements, J are the interelement tractions and r are the residual forces and can be obtained from (6). A particular expression derived in Ref. 10 for two

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined 247

.Uc.:}~ 5~;rzeT,zt,,:,r~ ~zth adaphce finite elemeT~t ,",n,~lysis: E. iiinton eta[.

dimensional problems defines an element contribution

5"1 = 5.~i( 1 - < -

~l)

O-S

iie;n= [ ~

rrrdf!.+ ~

j J ' j d [ ]1/2

e

e

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

E (1 -

'~'3 = ( q 2 ( l

(8)

- 4 - rT)

(10)

gj are unknown nodeless coefficients associated with the distribution of the displacement error p within the element. To obtain these coefficients we solve the stiffness equations for each element which have the form

u)

7.

K g = fr

(11)

and for plane strain problems where the stiffness matrix I< h ~ the form K=

E (: + u ) ( : - 2 u )

I< = / ~

where E is the "t%ung's modulus and u is the Poisson's ratio. ]'he interior residual term in (7) is the dominant error term when shape functions consist of piecewise biquadratic functions and the inter-element traction j u m p is the dominant error term in (7) when the shape fimctions consist of piecewise bilinear functions. Babuska and Yu t~" have recently shown that for odddegree 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.

[L1Q]rD[LtN']d~

(12)

1

and 1q = IN,, N2,1qa]. The consistent nodal forces fr may be expressed as fr = / a - [1N]rrdf~

(13)

1

and the vector of the unknown

coefficients g

[g T ,g2, T g ]T.

=

The strain energy error for element i ~Lssociated with the quartie displacements is obtained for the expression Uell~,= ~grI
(14)

5.3 Residual and hierarchical based error estimators Another residual error estimator which makes use of r with a solution based on hierarchical bubble functions has been described by Baehmann el al la and is now presented. In this scheme we first evaluate the interior residual body forces r obtained when the FE stresses are substituted into the governing differential equation (6). A beller solulion is assumed involving hierarchical shape functions of a higher order than the ones used in the FB solution; they also have zero values on the boundaries of the element. The interior residual body forces r are then applied to each element and the hierarchical displacements may then be evaluated for each element independently. This leads to a better estimate of the strain energy which in turn leads to an estimate of the strain energy error of the original FE solution. We demonstrate this error estimator using the sixnoded triangle and assume that there exists a better solution based on quartic variations of the displacements over the element. Thus we represent the error in the displacements p using hierarchical (quartic) l)ubble functions so that 3

v =

r:bg;

(9)

j=l

where 1Qj = iv'jI2 and the shape functions Nj are written asl3

248

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined

For a mesh o f n elements, the total error energy which is used to evaluate the global percentage error is obtained from I ' l

[[el[~,

Ilell =

(15)

i=1

6.

ERROR

STRESS

ESTIMATION

-

BEST

GUESS

METHODS

6.1 Zienkiewicz - Zhu error estimator Recently, Zienkiewicz and Zhu 14 have developed an error estimator which, in its standard form, makes use of projection techniques to define a continuous stress field over the domain *+. Tile basic ideas in the derivation of this error estinaator are: 1. The use of a continuous stress field which eliminates the stress j u m p around the element boundary. 2. The use of specific 'superconvergent' points in the element for tile stress calculation where predicted FE results are more accurate than the other points. They evaluate the error energy norm 14 from the expression, ; t i t is a s s u m e d t h a t a p p r o p r i a t e s t r e s s e s a r e d i s c o n t i n u o u s across boundaries between two different materials.

3resh ,ler,:ration with adaptive finite element analy.sis: E. Hmton et al.

llelt

=

(~ - 6")TD-X(o" -

&) df~

(16)

In this expression as we do not have access to the exact solutions for the stresses o-, we substitute our best guesses ~r', so that o" ~ ~"

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. We can derive, from the governing equations of equilibrium, the expression for nodal equilibrium n B T o " dr2 - f = 0

(17)

6.2 Best guess stresses

(22)

ttowever, when the continuous stress field ~r" is substituted into (22), residual nodal forces q are obtained

The best guess values for (16) for standard isoparametric elements can be determined in several ways: nodal averaging, least squares smoothing and Loubignac iteration. The smoothing procedures are now briefly described.

f Br~r'df2 - f = q

(23)

6.2. I Smoothing procedures

These residual forces can be minimised by the application of Loubignac's iteration scheme ls,lz. This procedure leads to continuous stress fields for which the residual forces q are minimised.

The component of the continuous stress field, which is to be found, is approximated as~

6.4 Ad hoc best guess stress for a Mindlin-Reissner eL ement

(is)

For other applications different procedures for obtaining best guess stresses (or more precisely best guess stress resultants) have evolved. For example, Selman et al is have demonstrated the use of tile AFM mesh generator and AMR procedures in accurately reproducing boundary layers in plates with certain types of boundary conditions. They use the six-noded triangular Mindlin-Reissner (MR) element introduced by Zienkiewicz el al 19 The error energy norm II e II for ,~rt plates is computed in terms of bending moments and shear forces and m a y be written as

i

where No, is a C ° continuous shape function and 6"i is the smoothed nodal stress associated with node /. The continuous stress a ° is determined by forcing the Galerkin-like condition ~ N~,(~r" - &) dr2 = 0

(19)

This is equivalent to the least squares smoothing procedure of Ilinton and Campbell is. One method which can be employed to obtain a continuous stress field using (18) and (19) is to assume that the stress a" is interpolated by the same shape function Ni as the individual displacement components. This yields S# - q = 0 (20) where typical components may be written as

Sij = j[NN?NjdQ and qi = ~Nif'df'2

(21)

For local least squares smoothing, (19) is solved by considering individual element domains separately. For elements other than the standard isoparametric element various ad hoc or empirical methods have been suggested. Two examples of such approaches will be presented later.

6.3 Loubignac algorithm Although the use of continuous stress fields eliminates the traction j u m p s between elements, it does not 4We dcM with eada stress componcnt (e.g. az, cry, r.u ) in

turn.

Ilell

= [ffl[M - 1QIITD~-'[M -- 1VIJdfl + ffl[S - S]TD21[S - S]df2 ]1/2

(24)

where M = [M,, My, M,y] T and 1(¢I= [~[~, ~[~, lQ, y]T are the exact and FE bcnding m o m e n t s respectively and S = [S~, Sy] r and f; = [5"~, oOy]T are the exact and FE shear forces respectively. Db and D , are matrices of bending and shear rigidities respectively. As we do not have access to the exact solutions for the bending moments M and shear forces S, we substitute our best guess values M" and S ' , so that M .~ M" and S .~ S*

(25)

Procedures for obtaining these best guess stress resultants have been found by experimentation. As the bending moments resulting from tile FE analysis are linear, the improved moments are interpolated with quadratic shape functions from the averaged nodal values. Thus the terms of the first integral in (24) are quartic polynomials and a seven-point integration rule is used in the evaluation of the bending contribution to II ell. The FE shear forces are also linear over the element

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined 249

Mesh generation a'zth adaptive finite element analysis: E. Hir~lorz et ai.

and are evaluated at the centroid of the element and a constant distribution over the element is assumed. Averaged nodal values are subsequently obtained at the vertices and a linear interpolation is then assumed. Therefore, the terms in the second integral in (24) are quadratic and a three-point integration rule is used in the evaluation of the transverse shear contribution to

z >

I1~11. The strain energy of the exact solution is estimated as:

Ilw" I1~ [fa[M']rDg~M'dfl (26)

+ fa[s']rDT~S'dn 1~/~-

nA @

6.5 Ad hoe best guess stress for.facet shell element For a triangular facet thin shell element 2°, we use an error energy norm expressed in terms of the bending moments and membrane forces

Iloll

fl B

= [ffl[M - 1Q]TD~'t[M - l~Ildfi + f n [ N - IQ]TDL~[N -- lqldfi]'/-"

(27)

where M = [Mz, My, M~:y]T and ~I = [f[j:, flu, f[~:y]T are the exact and FE bending moments respectively and N = [N,, Ny, N,~] T and N = [/~'~, ]gy, .~,ylT are the exact and FE membrane forces respectively. Db and D,n are the matrices of flexural and membrane rigidities respectively. As we do not have access to the exact solutions for the bending moments M and membrane forces N, we substitute our best guesses M" and N ' , so that M~M"

and N ~ N "

Fig. 1I Definition of tangent and normal coordinate directions at midside node of element sides of facet element.

4. For each element mid-side, transform the bending moments and membrane forces to the local n,t coordinate system, so that l~I = [f&, 2Q~, and lq = [N,~, Ne, i~rnt]T.

5&,]T

5. Obtain average values M" = [Mg, M ( , M2,] T and N" = [b~, N;, N2,] T at each element midside.

(28)

To evaluate the best guess values we use the following procedure: 1. Evaluate FE bending moments 1VI and membrane forces N for each element. Note that these values are constant over each element. 2. At each element mid-side determine the local tangent direction ' t ' the positive direction of which is given by a vector coinciding with the element side and pointing to the corner node with the higher number. (Note that this provides a unique way of defining the positive tangent direction for each element side.) 3. Set up the local normal direction 'n' at the element midside. Note that for adjacent elements the two normal directions will not coincide exactly but the projection of one normal onto the extended plane of the adjacent element will coincide with the other n o r m a l - see Figure 11.

250

r

The error energy norm Ilell may then be evaluated using a three-point quadrature rule. As I[el[ is a scalar quantity and as the quantities I ~ I , M ' , N and N" are already known at the mid-side sampling points, very little extra effort is required. The strain energy of the exact solution []w']l" is estimated as w

')

IIw'll ~ = IIw ltx~- + IIw

~

IIN9

(29)

where

IIw'll~r ~ /a[M'ITDbl[M'] dfl and

IIw'll~ ~ /r [N']TD~I[N'] dfl Note that this is different to the way in which Zienkiewicz and Zhu evaluate llw'll.

Adv. Eng. Software, 1991, t/ol. 13, No. 5/6 combined

Mesh generation with adaptive finite element analysis: E. Hinton et al

7. R E F I N E M E N T

PROCEDURE -

The error energy norms presented iu the previous sections are measures of the absolute value of the predicted error over the domain. The predicted values of the percentage global error r/is given by the expression

llell The quality expression

IIw It = (/i"l ¢rrD-Xo " dgt ) 1 , 2

~/2

(38)

Tie/

7.1 Estimation of global error

0Ilwll can be computed Ilwll

IIw-II

where n~t is the total number of elements in the domain under investigation. As the error is evaluated at the element level, we define the parameter ~i for each element as ¢, = Ilell, (39)

p

(30) from

the

where ~i is called a 'refinement indicator'. Depending on the values of ~i we can identify three possible element states: o if (i = 1 we have the optimal element size

(31)

* if ~i > 1 refinement is necessary, and In Reference 14 the exact strain energy norm of the solution is approximated by

o if (i < 1 de-refinement is possible.

7.3 Evaluation of the mesh density

Ilwll ~ IIw* II

(32)

The characteristic feature of the AMR procedure is the use of the current solution to predict the new mesh size hl- For instance, if the current element size is hi and the rate of convergence of the adopted element is O(h l) then we can design the new element size to be hi which is given by the expression

where

IIw'lt=

(11,~11-~ + II~llb w-

(33)

in which II ~" 11 is the FE approximation to the strain energy norm and is defined as

=

h i - ~/~ hl

(34)

Alternatively Sibai and Hinton 2~ approximate ]1 w* ][ using the smoothed continuous stress field

w

=

(35)

since the strain energy of the FE solution [[ ,~" [[2 can be upper or lower bounded by the exact strain energy. The estimation given by (33) is not valid if II '~ I12 is lower bounded. Equation (30) allows an effective adaptive process to be developed with the principal objective of achieving a specified overall percentage accuracy f / - - say 5% in many engineering applications - - in the energy norm, i.e.

< f]

(36)

or noting (30), as well

I1oll _< ~ IIw" II

(37)

7.2 Refinement indicators As the ultimate aim in AMR procedures is to obtain a uniform error distribution for all elements, the permissible error for each element is determined by ¶ Note that for convenience here we provide expressions for Iiwll for 2D and 3D elasticity situation ~ for plates and shells appropriate expressions should be inserted.

(40)

The values of g depends on the smoothness of the solution and on the norm used to evaluate the error. Generally £ is taken to be equal to the polynomial degree of the FE approximation, tlowever, near a singularity it has been shown that 0.5 _< g < 1 where the value of g represents the strength of the singularity. Equation (40) can be used to evaluate the design mesh density for an automatic mesh generator. For the facet shell element 2°, which is a combination of the linear plane and the quadratic Morley plate element, the value of g should represent the characterstics of the dominant element. The value of £ in (40) required to design the new element size can be estimated using

e = p - I f " " I1~ /II w" 112

(41)

where p is the polynomial degree for the lateral displacement of the Morley plate element - - in other words p=2

7.4 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 as

0_ Ilell ]]E][

(42)

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined 251

Mesh generation with adaptive finite element analysis. £. Hinton et al.

8. I AMR for plane strain problem

estimators we consider a plate with a crack under uniform tensile load"" as shown in Fig. 12. The following properties and dimensions are used: elastic modulus E = 1.0, Poisson's ratio L, = 0.3, side length L=2.0, length of the crack ~=0.5 and edge pressure intensity of p = l . 0 II. Tile exact strain energy for this problem is known to be 1.468762 ~2. In these comparisons a 6noded triangular element is used exclusively. The initial mesh on the cracked plate domain is also shown in Fig. 12 and has 22 elements and 108 degrees of freedom. A prescribed error of ~ = 1% is assumed to compare the different error estimators. Figure 13 shows the effec tivity index versus the number of degrees of freedom and Fig. 14 illustrates tile corresponding rates of convergence obtained for the best guess stress and residual type of error estinaators. 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 almost identical. Table 1 summarises the results for this problem. In this table DOF means number of degrees of freedom required, tl ~ II" is the unsmoothed strain energy, 0 is the effectivity index, q is the global percentage error and subscripts i and f refer to the values at the initial and final adaptivity iterations respectively.

To demonstrate the efficient use of the AFM mesh generator for adaptive FE analysis of plane strain problems and to compare the performance of different error

used.

where II°II is the estimated error and "iiei] - ' is the exact error. As the adaptive solui.ion proceeds 0 should converge to 1.

7.5 Determination of 6 for mesh generator Once the values of hi a r e determined the new values of 8 are calculated at the nodes of the current mesh using the expression -rlel .=

(43)

i

where n,i is the total number of elements surrounding node k.

8. E X A M P L E S

OF AMR

WITH

AFM

MESH

GENERATOR

In the previous sections we described the methodology which may be used to estimate the error energy norm 11e 1[ and evaluate the mesh paranaeters 5 for the mesh generator. We now present some applications, where these error estimators have been used.

[]In these examples, a

p=l.0

l

p=l.0

tt tit

0 . 5 _,_

C

tttttt

1.0

vF

-i o

c,,i

i

self-consistent

2.0

I;

ga

i

L v I

Fig. 12 Plate with a crack." problem definition and initial mesh.

252

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined

set

of units

have been

Mesh generation with adaptive finite element analysis: E. Hinton et al.

1.I0

Table 1 -- Plate with a crack under uniform tensile load

i

I !

1.00

method

0.90

/e / 9""

i

L

0.70

io, Lof!

n. J n! I C P U (sec) i

z z (,~x)

2350 1.46863 0.75 0.98 16.9 i0.84 i

zz (L)

2308 1.46863 0.69 0,89 15.7110.76 I

/

°...y

[

14.44

1540 1.46852 0.46 I1 0,83 9.83 i 0.89

4.32

1252 P 1.46841 0.3810.74 i82410.84 l

3.41

:'

r"

0.40

'

I

500

'

1000

R (H)

/500

[

rt (i)

LEGEND • :. ZZ(NA) •.~,.. Z Z ( L ) ZZ(Lg/NA) • ~<. Z Z ( G ) . ~ , - R(I) -~° R(H)

2000

2500

5.14 6.82

2182 1.46863 0.64 0.92 14.18 i 0.77 i I ZZ (Lg/NA) 3000 1.46814 0,9310.96 20.16 i 1,93

1

........... ~_i~ .-'-~'' 0.50

W!i~

zz (G)

08O

.~

DOFf

j

i

17,5,t

Figure 15 shows the final meshes obtained with the best guess stress methods with different smoothing procedures. Figure 16 illustrates the final meshes obtained using the residual methods.

3000

3500

degrees of freedom Fig. 13 Plate with a crack: effeclivily inde:c versus degrees of freedom for Zienkiewicz-Zhu with nodal averaging (ZZ(NA)), Zienkiewicz-Zhu with local smoothing (ZZ(L)), Zienkiewicz-Zhu with global least squares smoothing (ZZ(G)), simple nodal averaging with Loubignac iteration (ZZ(Lg/NA), residual based error estimator (R(D) and residual and hierarchical based error estimator (R(H))

(a) ~

DOF I

DOF!

(b)

~ .........

l':'::

?

'--...~ •..'i-77: ,

(e) ~

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

DOF/

0!

rh

7/I

LEGEND [] Z.Z(NA) ..4~.. Z Z ( L ) • . a . Z,Z(Lg/NA)

":t,,L N4%

• ~<. ZZ(G) - ~ - R(I) -¢~- R(H)

200

ti

".20

240

0/

rh

rtt

) J

260

0i

2182 0.64 0.92 14.2 0.77

q.60

-220

0i

"'"-.A ".:..,

;~ "%

.2.00

nI

~ .....I

-I.40

.1.80

ni

E *-..

-1.20

~

01

2308 0.69 0.89 15.7 0.76

-0.80

~

i

2350 0.75 i 0.98 16.9 0.84

-060

-1.00

0~

2.80

I

3.00

w 3 20

-340t

(d) 360

log N

) ) )\

DOF!

0~

Of

r/i

rH

3000 0.93 0.96 20.2 1.93

5

Fig. 14 Plate with a crack: convergence curve for different error estimators. Fig. 15 Final meshes obtained using Zienkiewicz-Zhu

error estimator for cracked plate problem: (a) ZZ(NA), (b) ZZ(L), (c) ZZ(G) and (d) ZZ(Lg/NA).

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined 2 5 3

Mesh generation with adaptive finite element analysis: E. Hinton et al.

(a) i

1540 o46 o83 9.83t0.891

(b)~"

• When we use the residual method, R(T) (equation (7)), the results are almost identical to the other residual method, R(tI). 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 but the effectivity index is far from 1.0. The final mesh is shown in Figure 16(b).

8.2 AMR for plate bending DOFI- Oi

0s

rh

'1!

1252 0.38 0.7,1 8.24 0.84

Fi 9. i6 Final meshes obtained using residual methods for cracked plate problem: (a) R(tl) and (b) R([).

The following points are noted: • The Zienkiewicz-Zhu method with simple nodal averaging, ZZ(NA), has a very good effectivity index value of 0.98 in the last iteration. It is the easiest method to implement. The final mesh for ZZ(NA) is shown in Fig. 15(a). o The Zienkiewicz-Zhu error estimator with local smoothing, ZZ(L), gives similar results but at a much smaller computational cost when it is compared with global smoothing. The final mesh distribution is presented in Fig. 15(b).

MR plate bending is another area where the design of suitable FE meshes is essential to highlight certain phenomena such as edge effects. To demonstrate the occurence of boundary layers we consider the adaptive FE analysis of a uniformly loaded square plate with soft simple supports (i.e. the lateral displacement w is constrained to zero around the plate boundary). The plate has a thickness to side length ratio of 0.05 and a Poisson's ratio of 0.3. From considerations of symmetry only the quadrant shown in Figure 17 is analysed. Figure 18 shows the mesh progression and the associated distribution of twisting moments M,~t and shear forces St obtained using the AMR procedure (n and t coordinates are the normal and tangential axes to the supported edge). The presence of the boundary layers for both M,~t and St is clearly shown and the ability of the AMR and in particular the AFM mesh generator to automatically grade the mesh appropriately is highlighted. If uniform mesh refinement (UMR) had been adopted many more elements would have been required to reproduce the edge effect.

Y

• The use of the Zienkiewicz-Zhu error estimator with global least squares smoothing, ZZ(G), appears to be more accurate than the other smoothing methods. The final mesh of the adaptivity analysis is shown in Fig. 15(c). o Simple nodal averaging (or other smoothing method) with Loubignac iteration, ZZ(Lg/NA), provies a method for reducing all equilibrium errors. In the Loubignac method the values of the effectivity index are very close to 1.0 in all adaptivity iterations. The final mesh obtained for this case is shown in Fig. 15(d). • The residual method with the higher order bubble shape function, R(tI), gives very impressive results. The rate of convergence of this method is excellent. Ilowever the effectivity index is lower than the ones obtained with the best guess stress methods. The effectivity index grows very fast with increasing degrees of freedom. The final mesh obtained is shown in Fig. 16(a).

254

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined

X

Fig. 17 Symmetric quadrant of a uniformly loaded square plate with 'soft' simple supports.

Mesh #eneralion with adaptive finite element analysis: E. Hinton et ah

mesh

M,~t

f

St

i

Fig. 18 Sequence of adaptive analyses showing meshes, lwisling moments, lt[nt and shear forces St. 8.3 AMR for thin shells To demonstrate AMR for surface structures we choose the triangular facet thin shell element ~°. In this problem, as we are dealing with curved surfaces, we use the procedure described earlier for surface structures to carry out the A M R procedure. We present two examples to demonstrate the coupling of the mesh generation and A M R procedures.

(a) Clamped quadratic shell under central point load The clamped quadratic shell shown in Fig. 19 is analysed. The span of the shell is taken as L = 6 and the thickness t=0.2. The following material properties have been assumed: elastic modulus E = 30 × 10 z and Poisson's ratio u = 0.3. The shell is subjected to a central concentrated load of intensity P = 1.0. Using s y m m e t r y only a quarter of the segment

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined 255

Mesh generation with adaptive finite element analysis: E. Hintoa et el.

Table

Y

2 -- C l a m p e d quadratic poi~:t load ~ith

x;-- : • ' 0

IS5

07214

n57

O +}7~#.

061

U6473

2 8~,57

462

71 :,8213

@T3

i) ~740

i074

0 I0gO

2 , 3 2 5 ~,

Q 74 {

0.0:;6':~

2.2'J73

-'3 iJ

00232

22910

6

0 57 t3

,

w" :

0 7~

. ',0 "

0 57!r3

~"

e

z , !3'

cer~trat

eieczents

I !22

~'

~h,~i[ u n d e r ANIR

.~.

• 10 ~

.";~

!31i1 13 1

1 1807

}

05735

074

{).57[7

,c~ f,' = a ' * : .

Table

e I e m e n t s I ';',~ i

3 -- C l a m p e d point

A =

074{ aw

quadratic load with

.... i7,.,, w ~,~

shell under UNIR

central

:': 1-0-~

128

0.71 i 0

450

0.6263

1152

0.5987

1800

0.5912

:~06"5 I

:

f'l 5 s g s

j

J . 10.72 I

i

0.3657

2.5052

i 24.9

1

0=5%, ), = ~ : i W H~': N / "ItW!L'4

Fig. 19 Clamped quadratic shell under central point load." symmelric quadrant used in analysis. is analysed. Figure 20 shows the progression of meshes o b t a i n e d a n d the associated error d i s t r i b u t i o n using the A M R procedure with 7) = 5%. T h e initial mesh c o n t a i n s l a 5 elements giving a percentage error of 31.0 %. T h e final

e l e m e n t s = 135 r/% = 31.0

e l e m e n t s = 462 0% = 13.4

mesh contains 1807 elements and tile error reduces to 6.4 %. T h e results o b t a i n e d using A M R are s u m m a r i s e d in Table 2. Similar studies were carried out for this p r o b l e m using U M R and Table 3 lists the results for this problem. Figure 21 shows the sequence of meshes and associated error distrib u t i o n using the U M R procedure. T h e rates of convergence are compared in Fig. 22. T h e convergence curves indicate the a d v a n t a g e of using A M R procedures.

e l e m e n t s = 1122 0% = 8.0

e l e m e n t s = 1807 r/% = 6.4

meshes

error distribw

I

,on

m [] Fig.

20 Clamped quadratic shell under central point load: sequence of meshes and associated error distributions obtained using AMR.

256

Adv. Eng. Software, 199i, Vol. 13, No. 5/6 combined

Mesh generation with adaptive finite element ana'ysis: E. Einton et al.

elements

elements = 450 ,7% = 24.9

elements = 128 ,7% = a7.9

=

elements

1152

--- 1 8 0 0

,7% = t4.6

~7% = 1 7 . 4

meshes

error

distribt

m

m m m~

Fig. 21 Clamped quadratic shell under central point load." sequence of meshes and associated error distributions obtained using UMR.

-I.90 -

Z

.2.00

-

-2.10 -

-220

.2JO

%

d

-f

-240

-2.50

.2.60

-270

-

= 25m, ~4= .5, d/t = 12.5

-

"m

-2.80

220

240

260

280

3.00

320

3.40

log N Fig. 22 Clamped quadratic shell under central point load: convergence plot AMR vs UMR. (b) Clamped hypar shell under uniformly distributed normal load T h e clamped hypar shell shown in Fig. 23 is analysed. T h e hypar shell is subjected

Fig. 23 Clamped hypar shell under uniformly dislributed normal load: symmetric segment used in analysis. to a uniformly distributed load p of intensity 10.0. Tile following material properties are assumed: Y o u n g ' s modulus E = 1.0 x 105 , and Poisson's ratio v = 0.4. O t h e r i m p o r t a n t dimensions are indicated in Figure 23. Figure 24 shows the sequence of meshes which are obtained when the A M R procedure is used with 0 = 5%. Table 4 lists the results

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined 2 5 7

Mesh qeneratton with adaptive finite element analysis: E. Hinton et al.

elements = 141 r/% = 18.5

elements -r7% = 10.5

390

elements = 793 r7% = 7.7

elements ~% --

6.2

7- 1185

meshes

error distribt

Fig.

-

--

7

2~/ Clamped hgpar shell under uniformlg distributed normal load: sequence of meshes and associated e r r o r distributions using A MR.

e l e m e n t s -" 144

elements = 400

e l e m e n t s = 784

elements = 1156

r/°-/o = 2 0 . 5

'7% = 1 4 . 3

~% = I 0 . 8

r/% = 9 . 1

meshes

de ir sr torri b u t i o n s #

Fig.

-- i.-~ :..f~.. :~::- . / . . : : j~ " ~. .: .:~.j .:- : : : : -j

:

j;:::~s::::

~5 Clamped hypar shell under uniformly distributed normal load: sequence of meshes and associated e r r o r distributions using UMR.

of this p r o b l e m . S t u d i e s with U M R were also carried out for this p r o b l e m a n d Fig. 25 shows the sequence of meshes and associated error

258

.:~z:::~;:i~:z: .....

Adv. Eng. Software, 1991, Vol. i3, No. 5/6 combined

d i s t r i b u t i o n s with U M R . T a b l e 5 s u m m a r i s e s tile results for this p r o b l e m using U M R . T h e r a t e s of convergence are c o m p a r e d in Fig. 26.

Mesh generation with adaptive finite element analysis: E. Hinton et al.

Table 4 - - C l a m p e d h y p a r shell u n d e r d i s t r i b u t e d normal load with AMR

It !i:

i

i

'

0.1367

390

0.1320

7.01

0.1304

793

0.1313

7.07

1185

0.1312

7.07~

i e x

L

0132

141

i929i

045 1 !18.5

7.58

0.1443

10.5

0.1305

7.34

0.0777

7.7

0.1306

]7.25!

0.0499

6.2

o = 5%, A = iwli~./;w'-.,, Table 5 - - C l a m p e d h y p a r shell u n d e r d i s t r i b u t e d n o r m a l load w i t h U M R

lO3 .%

'ieii: x i 10.27

0.5575

20.5

0.1312

8.29

0.2698

14.3

0.1310

I 7.71

0.1522

10.8

0.1309

! 7.52

0.1073

9.1

144

0.1380 17.11 I

0.1324

400

0.1339

17 0 ~

784

0.1325

7 04

0.1320

17.05

1156

1.40

The shape of the structure is modified by changing certain parameters which are known as design variables. It is desirable to have an evident connection between the values of the design variables and the actual geometry. Often the main aim is to minimise the weight or volume of the structure although other objective functions may be used (e.g. stress levelfing). The optimization is carried out subject to certain constraints which may, for example, be designed to limit the equivalent stress. Also certain geometric constraints are often employed to avoid buckling, vibration, or manufacturing problems. In shape optimisation the original mesh gradually changes as tile optimisation proceeds and tile structural geometry changes. Unfortunately, these changes in the structural geometry can lead to gross element distortions and also to substantial discretisation errors. The problem of distortion and discretisation errors can be solved by integrating the AFM mesh generator, AMR and shape optimisation procedures. To achieve a fully efficient and integrated approach the potentials of the AFM mesh generator are exploited in the following ways: • By using tile key points defined earlier as design variables.

130

By using, for shape optimisation, the same shape definition of boundaries adopted ill the mesh generator.

120 "4

~

II0

~,,~

_q

By moOifying the output data from the optimiser to match with the input data required by the mesh generator so that the process becomes wholly automatic.

........"'" ~ 1"""''''"'~)

1.00 LEGEND [] A M R •. ~ . . U M R

0.90

I \

• By allowing for remeshing to avoid element distortion if necessary.

f '

I 2.20

'

2 30

240

2-~0

I

'

2.60

l

2.70

'

I

2.80

'

i

2.90

'

I

'

300

'

• By integrating adaptive analysis and shape optimisation procedures.

310

log N Fig. 26 Clamped hypar shell under uniformly distributed normal load: convergence plot A M R vs UMR.

The coupling of the mesh generator, AMR and shape optimisation procedures is illustrated in the flow chart shown in Fig. 27. 9.1 Connectin 9 rod ezample

9. SHAPE OPTIMISATION MESH GENERATOR

WITH

AFM

VCe now present an example from structural shape optimisation, to demonstrate the versatility of the AFM mesh generator. Shape optimisation involves the determination of an appropriate geometric shape for a structure so that it can carry the imposed loads safely and economically.

To demonstrate the efficiency of the above integration we consider the optimisation of a connecting rod. The objective in this problem is to mininaise the weight subject to the equivalent stress remaining below a certain critical value and to use FE meshes which lead to results of prescribed accuracy. An engine connecting rod is subjected to a distributed load varying linearly as its direction of application changes from 90 o to 180 o as shown in Fig. 28.

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined 259

:~,[esh generation with adaptive finite element analysis: E. tfinton

: L.

OpI'im~zahon Parameters Design Parameters

et al.

i

I T

HESHGENERATO IN

~ 1 .

Update boundary condihons I I Update loading I I

.

.

.

.

.

v

.

ANALYSIS

ANALYSIS{FEN) !i '

Eshmafe

Flew

element

size i

t

i /

.

A D A P T IV E

STRUCTURAL \

adaphv,ty ~

<,I-

.

J

i

__

.

:

f YES

Error esHmafe

! I

/Is ...... ---'~ error

~ , ,

~..accepIab~

'"-"---

NO

["'STRUCTURALSENSITIVITY t! II I

Obje(:tive derivatives Constraintsd e r i v a t i v e s

•

I i t

r

OPTIMISATION

L i

Update design variables

L

Z50 12

-- i *,~.

12 S 11

,L2"~

T L

I4

Z2 S

~11

5

ANALYSIS

i

Fig. 27 Flow chart illustrating the integration of AFM mesh generator, AMR and shape optimisation procedures.

75

~0

~D

J

-3--'--

9

~"

iI

Jl

4t

h3

,lz

~

ho ~9

i

r

A

~

/#

k

,~.: z2 ~oE6 .... 2 , "~'',//'ll///).;'.;'/~//W':';-C.22>C';'.'.-y/C;;:.s->;;:)'2 ~'Y//~,7/J'>;'-//II/'/~;'~7 ,-

.

.

.

.

.

.

"

6o

,

(.,)

~5

3o

!

(h)

All dlmens~ons in mm

Fig. 28 Connecting rod example: (a) problem defthilton, (b) location of design variables and constraint points.

260

•

• I

Adv. Eng. Software, 1991, Vol. lg, No. 5//6 combined

Design variables C o n s t r a i n t points

Mesh generation with adaptive finite element analysis: E. Hinton

etal

(a)

40000.00

'~' feasiblesolution

35000.00

30000.00

25000.00

(b)

~

20000.00

nrr1

I0000.00

5000.00

I

I

I

2

4

6

'

[

'

8

I

I

10

12

14

ITERATION NO. Fig. 29 Optimisation process of connecting rod problem: (a) boundary shape changes and (b) design history.

Equivalentstresses

Mesh

V Fig. 30 Initial and final mcshes and associated equivalent stress distributions for connecting rod problem.

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined

261

Mesh generation with adaptive finite element analysis: E. IIinton et al.

The following properties and dimensions are adopted: Young's modulus of elasticity g=210000, Poisson's ratio u = 0.3, thickness t=10.0, and linearly distributed loading intensity of Po = 68. By' taking advantage of symmetry only one half of the connecting rod is analysed. The allowable equivalent stress constraint used for the analysis is chosen as 200. Figure 29 shows the boundary' changes and corresponding design history. For a 3% prescribed accuracy and using a residual based error estinaator, the initial and final meshes and the corresponding equivalent stress distribution are presented in Fig. 30. The initial volume of 35473.6 mm 3 reduces to final volume of 3815.4 sam 3 in 13 optimisation iterations.

10. C O N C L U S I O N S The main features of the AFM mesh generator and its extension to mesh generation on curved surfaces have been described. Algorithms for adaptive FE analysis have been presented indicating the importance of mesh generation in the overall process. Various methodologies for evaluating the error energy norm II e II (required to compute the new elenlent sizes for the mesh generator) have been presented for different applications. Several examples have been used to demonstrate the potential of this mesh generator for adaptive analysis of planar and surface structures. All the examples indicate the efficient coupling achieved between the mesh generation process and AMR procedure. The application to structural shape optimisation has also been shown illustrating the versatility of the AFM mesh generator. Other possible areas of application of this mesh generator include the simulation of solidification, contact-impact, localisation and crash mechanics where remeshing has to be performed fre-

of Swansea, Sial', 1991 4. Pelro, J., A finite element procedure for the solution of the euler equations on unstructured meshes, Ph.D. Thesis. Dept. of Cir. Eng., University College of Swan.sea, 1989 5. Zhu, J.Z., Zienkiewicz, O.C., Hinton, E. and Wu, J., A new approach to the development of automatic quadrilateral mesh generation. Int. J. Num. Moth. Eng., (to be published) 6. Zhu, J.Z., Hinton, E. and Zienkiewicz, O.C., Adaptive finite element anal)sis with quadrilaterals, Computers 0 Structures, 40, 1097-1104, 1991 7. Fau×, I,D. and Pratt, M.J., Computational geometry for design and manufacture, Ellis tlorwood, Chlchester 1979 8. Lohner, R. and Parikh, P., Generation of three-dimensioned unstructured grids by the advancing front method, Int. J. Num. M-eth. Fluids., 8, 1135-1149, 1988 9. Ozakca, M., Rao, N.V.R., and Hinton, E., A comparison of error estimators for triangular elements, Int. report, Univ. Col. of Swansea, April 1991 10. Babuska, I. and Rheinboldt, W.C., A posterlori error estimators for the finite element method, Int. J. Num. Moth. Eng., 12, 1597-1615, 1978 11. Kelly, D.W., Gago, J.P.de.S.R., and Zienkiewicz, O.C., A posteriorl error analysis and adaptive processes in the finite element method, Part I: error analysis Int. J. Num. Moth. Eng., 19, 1593-1619, 1983 12. Babusk~, I. and Yu, D., Asymptotically exact a ponteriori error estimator for biqaadratic elements, Tech. Report, Institute for Physical Science and Technology, Laboratory for Numerical Analysis, University of Maryland, Technical Note BN-1050, 1986 13. Baehmarm, P.L., Shephard, M.S. and Flaherty, E.J., A posteriori error estimator for triangular and tetrahederal quadratic elements using interior residuals, SCOREC, Report, Scientific Computation Research Center, Rensselacr Polytechnic Institute, 1990 14. Zienkiewicz, O.C. and Zhu, J.Z., A simple error estimator and adaptive procedure for practical engineering analysis, 1at. J. Num. Moth. Eng., 24, 337-357, 1987

quently.

15. tlinton, E. and Campbell, J.S., Local and Global smoothing of discontinuous finite element functions using a least squares method, Int. J. Num. Moth. Eng., 8, 461-480, 197,t

ACKNOWLEDGEMENTS

16. Loubignac, G., Cantin, G. and Touzot, G., Continuous stress field in finite element analysis, AIAA journal, 17, 16,t5-16.t7, 1977

The authors wish to thank Dr. J. Peraire for the use of an automatic mesh generation program for triangular elements and Dr. J.Z. Zhu and A. Selman for useful discussions.

REFERENCES

1. PerMre, J., Vahdati, M., Morgan, K. and Zienkiewicz, O.C., Adaptive remeshing for compressible flow computations, J. Comp. Phys., 72, 449-466, 1987 2. Shephard, M.S., Approaches to the automatic generation and control of finite element meshes, Appl. Mech. Roy., 41, 169-185, 1988 3. Sienz, J., Finite element mesh design with adaptive proce. dures, Msc. Thesis, Dept. of Civil Eng. University college

262

Adv. Eng. Software, 1991, Vol. 13, No. 5/6 combined

17. Cantln, G., Loubignac, G. and Touzot, G., An iterative algorithm to build continuous stress and displacement solutions, Int. J. Num. Math. Eng., 12, 1493-1506, 1978 18. Selman, A., Hinton, E., and Atamaz-Sibai, W., Edge effects in Mindlln-Relssner plates using adaptive mesh refinement, Eng. Comput., 7, 217-226, 1990 19. Zienkiewicz, O.C., Taylor, R.L., Papadopoulos, P. and Onate, E., Plate bending elements with discrete constraints: new triangular elements, Comp. Struct., 35, 505522, 1990 20. Dawe, D.J., Shell analysis using a simple facet element, J. Strain Anal., 7, 266-270, 1972 21. Sibai, W.A. and Hinton, E., Adaptive mesh refinement with the Morley plate element, Proc. NUMETA 90 Conf., Swanson, Jan. 1990 22. Shephard, NI.S., Niu, Q.X. and Baehmann, P.L., Some results using stress projectors for error indication and estimation. SCOREC Report No. 4-1989, Scientific Computation Research Center, Renssela, er Polytechnic Institute, 1989