Evolutionary structural optimisation based on boundary representation of NURBS. Part II: 3D algorithms

Computers and Structures 83 (2005) 1917–1929 www.elsevier.com/locate/compstruc

Evolutionary structural optimisation based on boundary representation of NURBS. Part II: 3D algorithms E. Cervera, J. Trevelyan

*

School of Engineering, University of Durham, South Road, Durham DH1 3LE, UK Received 21 June 2004; accepted 16 February 2005 Available online 19 April 2005

Abstract This paper presents an optimisation algorithm applied to three-dimensional structural optimisation problems. This algorithm is based on the evolutionary structural optimisation method (ESO) and boundary elements (BE) are used to carry out the structural analysis. The geometries are represented using non-uniform rational B-spline (NURBS) curves and surfaces. Examples are presented to show the effectiveness of the algorithm for some preliminary results.
2005 Elsevier Ltd. All rights reserved. Keywords: Evolutionary structural optimisation; 3D structural optimisation; Boundary element method; NURBS

1. Introduction There are different techniques used in 3D structural optimisation problems, most of which are a straightforward extension from 2D optimisation approaches. This paper presents such an extension to 3D of an algorithm that has been successfully applied for planar problems [1]. Deterministic-gradient based methods applied to 3D can be divided into gradient-based and gradientless methods. In both, the finite element [2] or boundary element [3] equations are used to provide stress results on a frame or continuum representation of the structure to be optimised. In the gradient approaches, the variational equilibrium equations are differentiated and discretised, and this applies equally to finite element [4] and boundary element [3] implementations. On the other hand, applications of gradientless, heuristic or scholastic meth* Corresponding author. Tel.: +44 191 3342522; fax: +44 191 3342390. E-mail address: [email protected] (J. Trevelyan).

ods in 3D include, amongst others, the evolutionary structural optimisation method (ESO) [5] and genetic algorithms [6]. An important feature in three-dimensional optimisation problems is the model representation and how it is related to the structural analysis results. Initially, some authors defined the nodal coordinates of the discrete finite/boundary [3] element model as design variables. However, this requires a large number of design variables and also it is difficult to maintain an adequate finite/boundary element mesh during the optimisation process [7]. Alternative methods appeared to overcome these initial drawbacks such as the mesh parameterisation [4] and solid modelling, which comprises constructive solid geometry (CSG) [8] and boundary representation (B-Rep) [9] models. In this paper a structural optimisation algorithm based on ESO [10] is applied to 3D problems. The boundary element method (BEM) is used to carry out the elastostatic structural analysis. In addition, the boundary shape is represented by NURBS curves and

0045-7949/$ - see front matter
2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2005.02.017

1918

E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1917–1929

surfaces, defined by control points. The NURBS representation is attractive since it can be exchanged with CAD/CAM systems. Moving the basic parameters of the NURBS (control points) can modify the geometry of the NURBS curve or surface.

2. Algorithm Similarly to the successfully implemented 2D optimisation algorithm [1] the optimum topology evolves by slowly removing and adding material from/to low and high-stressed areas, respectively. The boundaries are represented using NURBS surfaces [11] defined by control point nets. Use is made of the locations of these control points as design variables of the problem. Fig. 1 shows a basic flow chart of the optimisation process. This process can be divided into several steps which are associated with either the structural analysis or the optimiser (boundary-ESO algorithm). Step 1 (Structural analysis): The geometry of the structure is defined. The initial design is subjected to a set of loads and constraints. Step 2 (Structural analysis): A boundary element analysis (BEA) follows the model description. Step 3 (Optimiser): The geometry is imported together with the structural analysis results. Step 4 (Optimiser): Removal of material. The nearest control points to the least stressed boundary nodes

are identified and moved in the direction of the normal vector to that area. Step 5 (Optimiser): Addition of material. If any boundary node is found with a stress higher than the yield stress or any other maximum stress criterion, then a similar process to removal is undertaken but, however, the opposite direction of movement of the control points results in a material addition. Step 6 (Optimiser): The output file is written for the modified geometry. Step 7: Such a procedure is repeated (from step 2) until the stopping criterion is reached. Step 8 (Structural analysis): Post-process final design.

2.1. Geometry definition NURBS surfaces are adopted to define the geometry in order to control the curvature and tangency of the changeable boundary. A control point net is specified for each surface. Adjacent surfaces are required to share control points over their coincident edge and therefore ensure continuity between surfaces throughout the process. 2.2. Boundary element model The boundary element method is used for the structural analysis of the 3D topologies. The boundary element software BEASY [12] is used though the

3D-OPTIMISER

START

Reads results file

STRUCTURAL ANALYSIS

Pre-Process

PROCESS NURBS DATA & ANALYSIS RESULTS REMOVE MATERIAL

MODEL DESCRIPTION

ADD MATERIAL

Structural Analysis

BEA

MODIFY GEOMETRY Writes data file Reads results file

Post-Process

PLOT RESULTS

IS STOPPING CRITERION SATISFIED ? Yes

STOP

Fig. 1. Flow chart of the basic optimisation process in 3D.

No

E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1917–1929

algorithm is of course transferable to other BEM codes. The use of commercial software requires some specific features, further described, in the problem definition and also while processing and managing the analysis data. If an open source BEM system is available, including automatic meshing and 3D analysis capability, the optimisation algorithm described herein can be integrated into the code, and greatly enhanced efficiency and easy of use can be achieved. This work is under development. 2.3. Removal and addition of material 2.3.1. Identifying inefficient areas Material would be removed from the structure if there is any node p that satisfies the equation ð1Þ

rp 6 RRrmax and added to the structure if any node p satisfies ðrp P ry Þ OR

ðrp P ARrmax Þ

ð2Þ

where rp is the node von Mises stress or any other selected criterion, rmax is the maximum von Mises stress or any other selected criterion, which varies as the optimisation progresses, ry is the yield stress or any other maximum stress criterion, RR is the removal ratio and AR is the addition ratio (0 6 RR, AR 6 1). Following the classical ESO procedure [10], if a steady state is reached RR is increased by the evolutionary rate for removal, ERR, as follows RRi ¼ RRi1 þ ERR

expresses the influence of the node on the associated control point. The closer the node to the control point the higher its influence and therefore the weight. Consequently the weights are calculated as follows aj x0j ¼ Pm1 ð4Þ l¼0 bl with aj the distance between the node and the control point. m is the number of nodes associated to the control point. The distances, b, are equal to the distances a, of these m nodes from the control point. Notice that for a unique node related to the control point (m = 1) this weight is 1 (x0j ¼ 1). In order to afford nodes close to the control point a higher weight, a reciprocal is introduced x00j ¼

1 x0j

ð5Þ

To avoid numerical instabilities there is a minimum value allowed for aj. This value, w, is a factor related to the maximum dimension of the structure. After numerical tests this value is defined in this work as w ¼ max dimension=100

ð6Þ

Finally, the weight xj x00j xj ¼ Pm1 l¼0

x00l

ð7Þ

and therefore as a result of defining xj in such a way, P m1 j¼0 xj ¼ 1.

ð3Þ

where i is the current iteration. Typical values for suitable ratios are RR0 = 0.01, ERR = 0.01, AR = 0.99. These values are determined from numerical experiments. For each selected node pj, the algorithm searches its nearest control point Pj. Therefore, each node has a control point related as shown in Fig. 2, where p0 is related to P0 situated a distance a0 away. Since nodes are likely to be more closely spaced than control points, the same control point will be associated to different nodes as shown also in Fig. 2 for P1. To consider such situations each node has associated with it a weight xj. This weight

S P0

p4

a0 p3

a3

p1

a5

a4

a1 a2 P1

2.3.2. Direction of movement The direction of movement of each control point is perpendicular to the NURBS surface defining the associated area, i.e. where the nodes (pj) are placed. The perpendicular movement (see Fig. 3 (a)) follows the normal vector nj. This vector is calculated using the analytical definition of NURBS surfaces [11], its derivatives, i.e. vj and uj, and vector relationships (cross-product uj · vj = nj). The movement is inwards when removing material and outwards when adding it. For situations in which a single control point is associated to a group of nodes, as is illustrated in Fig. 4, the direction of movement n is calculated using the weight factors (Eq. (7)) and mathematically stated n¼

p5

p0

1919

p2

Fig. 2. Selection of control points associated to nodes.

m1 X

xj nj

ð8Þ

j¼0

where nj is the normal vector to the node j and xj is its associated weight. In the case of nodes placed on edges between surfaces (Fig. 3(b)) an averaged normal is calculated. Therefore, knowing the vectors nSj and nTj which are the normal vectors at pj to the surfaces S and T, respectively, the direction of the new normal vector at pj is calculated as that

1920

E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1917–1929

S P

nSj

uj, vj ∈ S nj ⊥ S

pj

nj

T

S

vj

⊥ S;

nj =

uj

(a)

nTj

⊥

pj

T nSj

n Sj + n Tj

nTj P nj

n Sj + n Tj (b)

Fig. 3. Direction of movement for control point P associated to node pj: (a) pj on surface S, (b) pj on surfaces S and T.

di ¼ p0 p5

S p4

a4

a5

1 AF ai P w ai

ð10Þ

T a0 a1

a2

P a3

n

For cases in which several nodes influence a control point (Fig. 4), a more general expression of the distance of movement for the control point Pi is given as follows " # m1 X 1 RF ð11Þ di ¼ aj j¼0

p1 p2

p3

2.4. Geometry control Fig. 4. Direction of movement for control point P associated to a group of nodes.

of a vector addition of vectors nSj and nTj . This criterion is not only applicable to two surfaces but to any number of surfaces sharing a node p. 2.3.3. Distance to move In any iteration the material is either removed or added to the structure by migrating control points. Each control point is moved a distance related to the following parameters [13]: • Distance of the control point from its associated node. This distance is denoted as a. • A factor related to the stress situation within the structure at the current iteration. This factor, associated to geometric dimensions of the problem, is called the removal factor (RF) if removing material, and the addition factor (AF) if adding it. For example a control point Pi situated at a distance ai from the node is moved a quantity di calculated as follows di ¼

1 RF ai

ai P w

ð9Þ

where w is defined in Eq. (6). Likewise, for addition of material

2.4.1. Smoothing algorithm It is found that as the optimisation process evolves, corners in the surfaces may become sharper, increasing the stresses in that region. This is an artefact of the method of averaging the direction of control point movement on edges (Eq. (8)). To overcome this drawback an algorithm to identify and to smooth sharp areas is implemented. Sharp edges are defined here as those for which, in the notation of Fig. 5, h < 45 or h > 180. The reason for choosing these limits is because any values out of this range would produce numerical errors in the meshing and Jacobian calculations. To smooth a corner the control point in corner (P) is moved to the midpoint of the line joining the previous (P1) and next (P2) control points as shown in Fig. 5. To calculate the angle in the corner the shape of the corner is approximated using two vectors u and v. The vector u points from P to its previous control point P1 whereas the vector v points from P to its following control point P2. 2.4.2. Corners effect Control points situated in corners are treated differently according to the subtended angle h as defined in Fig. 5(a). The motivation for this is because moving control points a uniform distance gives rise to undesirable distortion that tends to sharpen corners as shown in Fig. 6(a).

E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1917–1929

A relationship is established between the movement of the outer control points and the inner control points by using the concept of serendipity finite element shape functions [14]. Thus, for a surface defined by a given number of boundary control points (numBCP), the coordinates (x, y, z) at any inner control point can be found from the values of the coordinates at the control points on the boundary of the control point net

θ ≥180 θ

v

P

CP2

CP2

s u

d/2

v d

P u CP1

CP1

(a)

1921

x¼

(b)

numBCP1 X

LSi xi

i¼0

Fig. 5. Smoothing edges and corners: (a) initial model, (b) final model.

y¼

numBCP1 X

LSi y i

ð13Þ

i¼0

The effect of moving material in corners which are very convex/concave or sharp is larger than in areas on an even surface. A factor c, related to the angle in the corner h, is introduced and thus, Eq. (9) is modified as follows d i ¼ cd i

z¼

numBCP1 X

LSi zi

i¼0

where LSi terms are the serendipity element shape functions and (xi, yi, zi) locates the ith boundary control point. In processes with no restrictions over a plane of movement the current position of the inner control points is also affected by the effect of moving adjacent control points on the boundaries. Therefore for these control points, the global movement can be split into two movements; i.e. the main movement and the propagation movement. The main movement follows the general Eqs. (9) and (10) and the propagation movement depends on the movement of the edges. Hence, for a surface defined by a given number of boundary control points (numBCP), the coordinates (x, y, z) at any inner control point can be defined as follows

ð12Þ

The factor c is calculated following a linear interpolation. Fig. 6(b) plots the values of c according to the angle h. 2.4.3. Mapping the control points In certain problems, the control point movement can be restricted to be performed exclusively in one plane, i.e. movement in XY, or XZ or YZ. This might be desirable, for example, to test a 3D algorithm by reproducing conditions found in 2D modelling [13], and therefore seems a good starting point for the evaluation and development of the 3D optimiser. For these situations control points on the boundaries of the control point net are moved in the plane of movement. Control points situated inside the net, although also moving in the plane, are moved in a manner influenced by the movement of the boundaries. This is necessary to ensure that the essential (u, v ) grid based layout of the control point is maintained even following significant distortion of the edges of the surface.

x ¼ xm þ kDxp ;

numBCP1 X

Dxp ¼

Li Dxi

i¼0

y ¼ y m þ kDy p ;

Dy p ¼

numBCP1 X

ð14Þ

Li Dy i

i¼0

z ¼ zm þ kDzp ;

Dzp ¼

numBCP1 X

Li Dzi

i¼0

2.5

P2,j

γ

P1,j d1

d2

2.0

P1,j+1 P2,j+1 d0 P0,j+1

1.5

P0,j 1.0

0.5 0

(a)

45

90

135

180

225

(b)

Fig. 6. (a) Illustration of sharpening a corner when d0 = d1 = d2, (b) corner effect factor c.

θ°

1922

E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1917–1929

where xm, ym, zm are the coordinates resulting from the main movement; k is a parameter that determines the influence of the propagation movement over the global movement; Dxp, Dyp, Dzp are the quantities related to the propagation movement; and the Li(n, n) terms are the shape functions calculated over the local coordinates n, g

where V refers to the volume of the structure and U to the strain energy which is defined in BEA notation using the equivalence of strain energy to the external work done over boundary C, i.e. Z 1 U¼ T u dC ð17Þ C 2

Li ðn; gÞ ¼ ðnmax  jn  ni jÞðgmax  jg  gi jÞ

where T are the tractions over the boundary and u the displacements over the part of the boundary C. Accordingly, the stopping criterion e related to an arbitrary quantity fU can be quantified as follows  iþ1  f  f i  ð18Þ e ¼  U i U  6 104 fU

ð15Þ

It should be noted that the global movement of the inner control points is mostly determined by the main movement. The propagation movement effect is smaller and simply avoids overlapping due to the faster movement of the edges.

3. Restart procedures Restart procedures are implemented updating the problem definition at specific intermediate steps along the process. Once the problem is redefined a BEA is carried out and the optimisation process continues with the updated topology. For example, to take control of localised geometry distortion, in 2D problems an algorithm for adding and removing control points [1] was implemented. However, in 3D, the possibilities for geometric distortion are more numerous and varied, and restarts are applied to insert and remove control points. Another application of the restart procedures is to avoid distorted meshes. In these cases, even a finer mesh would not improve the convergence of the solution. Finally, restart procedures are used to perform topology changes. Inner low-stressed areas are removed from the structure by creating cavities. These areas are identified when low stressed nodes are present in the same region in two opposite close surfaces. Since elastostatic problems are assumed and no body forces are present, the material between the two facing surfaces can be considered to be at the same low stress levels. In such cases, the changes to the higher level NURBS surface description become significant and the new topology is accommodated by a totally new NURBS based boundary representation. The complicated issues of automated surface management remain to be addressed and are the subject of further research.

where fUi is the value of the quantity at iteration i and fUiþ1 is the value at iteration i + 1. Quantities used in this way in the current work include strain energy, stiffness and uniformity of stresses over portions of the boundary.

5. Examples 5.1. Beam under vertical load subject to restricted movement This simple example shows the optimisation of a beam with movement restricted to the plane YZ. Since the movement is constrained to one direction control points are mapped as explained in Section 2.4.3. The following isotropic material properties are assumed: YoungÕs modulus E = 210,000 N/mm2, PoissonÕs ratio m = 0 (fictitious to avoid undesirable stress concentrations at the fixed supports). Fig. 7 depicts the problem definition. The beam (100 mm · 50 mm · 20 mm) is described by surfaces S1 to S15, the surfaces S1 to S3 having prescribed boundary conditions and forming the non-design domain. It is fixed at the top and bottom surfaces (S1, S2) on the plane Y = 0. A vertical load of 100 N is applied on the surface (S3) situated in the middle of the plane Y = 50 mm. The objective is to minimise f = UV (Eq. (16)). The von Mises stress contour plots every 10 iterations are displayed in Fig. 8. The final design is found at iteration 43 for a volume ratio of V/V0 = 56%. 5.2. Beam under vertical load without movement restriction

4. Stopping criteria The stopping criterion adopted in 3D optimisation problems is related to the monitoring of the specific strain energy fU in the structure, i.e. fU ¼ UV

ð16Þ

The same example as that in Section 5.1, of a beam under vertical load, is now optimised to minimise fU (Eq. (16)). However, in this case, there are no restrictions on the 3D movement. The material properties are YoungÕs modulus E = 210,000 N/mm2 and PoissonÕs ratio m = 0.3.

E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1917–1929

Fig. 7. Beam problem: (a) geometry definition, (b) surface definition.

The evolution of the function (strain energy · volume) is plotted in Fig. 9. After 6 iterations the geometry is redefined in the restart step, the new geometry at this stage is shown in Fig. 10(b). Furthermore, a new restart is performed for reasons of geometric distortion at iter-

1923

ation 21, defined as restart2 and geometrically depicted in Fig. 10(c). Finally, the process evolves until the stopping criterion (Eq. (18)) is satisfied at iteration 44 (Fig. 10(d)). As shown in Fig. 9, further iterations do not improve the strain energy function. Control points situated inside the control point net are influenced by the movement of the control points on the boundaries. Hence, their movement follows Eq. (14) using the parameter k = 0.5. The RR is set up initially to 0.05 and increases steadily to the value of 0.35 when the optimum is reached (iteration 44). To control the convergence the factor RF is also considered. Initially, this factor is 1.0 but through the optimisation process this value is reduced to RF = 0.5, as a result of the stress levels that develop. Fig. 11 depicts the von Mises stress contour plots at different stages of the process. The optimum design is obtained after 45 iterations for a volume ratio V/V0 = 0.40, i.e. 60% reduction from the initial design. At iteration 44 the value of the von Mises stress for most of the material is between 40% and 50% of the maximum von Mises stress within the structure, which occurs in a localised fashion at boundary condition discontinuities.

Fig. 8. von Mises stress contour plots: isometric view.

1924

E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1917–1929

5.0E+07

f = U *V (Nmm4)

4.8E+07 4.6E+07 4.4E+07 Restart

4.2E+07 4.0E+07

Restart 2

3.8E+07

optimum

3.6E+07 3.4E+07 3.2E+07 0

5

10

15

20

25

30

35

40

45

50

Iteration

Fig. 9. Strain energy function evolution.

5.3. Beam under bending This final example shows the optimisation of a beam under bending. Isotropic material properties are assumed (E = 210,000 N/mm2; m = 0.3). The dimensions of the beam in this case are 100 mm · 50 mm · 30 mm. The problem is defined in Fig. 12(a), the initial geometry is divided into 10 surfaces as dictated by the requirement to isolate design and non-design domains in different surfaces. Surface S1 is fixed and surface S2 has a horizontal load applied of 1000 N. The optimisation follows a stress driven criterion which aims to minimise the function f = UV (Eq. (16)). The initial number of quadratic boundary elements is 237 (see Fig. 12(b)). However, the boundary element mesh is updated at each iteration, and the final element mesh consists of 504 elements. Fig. 13 shows the evolution of the objective function. Initially a RR of 0.05 is set; this ratio is increased gradually by the ERR(ERR = 0.01) (Eq. (3)). At the end of the process RR has increased to 0.11. The AR is set constant and equal to 0.90. Moreover from Fig. 13 it can be

seen that restarts are performed at iteration 5 (Restart 1) and iteration 21 (Restart 2), to avoid meshing problems caused by excessive distortion of a surface. Fig. 14 shows a wireframe view of the geometry for the different critical steps. Following the evolution of the objective function, the feasible minimum (for shape optimisation) is found at iteration 42, since further iterations result in a very distorted shape. Fig. 15 displays the optimum feasible design obtained at iteration 42. At this iteration there has been a 38% reduction from the initial volume (V0 = 150,000 mm3). The display of the elements (Fig. 15(b)) makes it evident that the section is undergoing a significant amount of thinning in two locations, as might be expected from the 2D study of the related plane stress case. This is the optimum for shape optimisation; topology optimisation would be accomplished doing a new redefinition of the surfaces. Looking at the lateral view in Fig. 15(b) there are some areas from the bottom and top surfaces that become closer, therefore surfaces between these areas would be potential sites for new cavities to develop. 5.3.1. Topology optimisation step In order to identify the areas to be removed the shape contours of both surfaces S3 and S4 are computed and displayed in figures Fig. 16(a) and (b), respectively. It happens that for both surfaces the convex/concave areas occur at similar coordinates (x, y). Thus, joining the isolines of maximum curvature in Fig. 17 two cylindrical cavities are created in the solid. It should be noted that the asymmetric solution of this symmetric problem is due to the automesher, since it produces non-symmetric meshes which later influence the movement of the control points. Fig. 18 shows the evolution of the whole process including topology optimisation. Topology changes are

Fig. 10. Boundary element mesh: (a) initial design, (b) restart, (c) restart 2, (d) optimum design.

E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1917–1929

1925

Fig. 11. von Mises stress contour plots: (a) isometric view, (b) lateral view (XY).

performed manually at Restart 3 (iteration 43), Restart 4 (iteration 56) and Restart 5 (iteration 59). These restart procedures create new cavities in the structure at inner low stressed areas in facing surfaces separated by a small distance. Inner low stressed areas must satisfy Eq. (1).

The minimum distance dmin is taken as a percentage of the minimum dimension of the problem, in this case dmin = 10 mm; 33% of 30 mm. It is observed in Fig. 18 that the effect of performing such a dramatic topology change, from one step to another, considerably

1926

E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1917–1929

Fig. 11 (continued)

Fig. 12. (a) Problem definition, (b) initial boundary element mesh.

4.7E+07

f = U*V (Nmm4)

4.5E+07

Restart 1

4.3E+07 Restart 2

4.1E+07

Feasible optimum Shape optimisation

3.9E+07 3.7E+07 RR=0.050

0.06R

R=0.08

RR=0.09R

R=0.10

RR=0.11

0.07

3.5E+07 0

5

10

15

20

25

30

35

40

45

Iteration Fig. 13. Evolution of the strain energy function.

increments the objective function. If care is not taken in restarting with a new topology with moderately small changes, the objective function can degrade severely to an unrecoverable state, as can be seen in Fig. 18.

Fig. 19 shows the von Mises stress at the critical steps of the process. Restart 4 performs changes enlarging the cavities and at Restart 5 two new cavities are created. From the related 2D plane stress case [1,13] it should

E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1917–1929

1927

Fig. 14. Critical steps of the optimisation process.

Fig. 17. 2D map of shape contours: (a) surface S3, (b) surface S4. Fig. 15. Final design.

Fig. 16. Shape contours: (a) surface S3, (b) surface S4.

1928

E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1917–1929

5.5E+07 Restart 4

5.3E+07

f = U*V (Nmm4)

5.1E+07 4.9E+07

Restart 5

4.7E+07

Restart 1

4.5E+07 Restart 3

Restart 2

4.3E+07 4.1E+07 3.9E+07 3.7E+07

Feasible optimumTopology optimisation

ther topology changes do not contribute to reach a minimum. Looking at the iterations of the process, we can observe the presence of very low stressed material near the fixed region (S1) which is not removed due to design constraints. This situation shows the clear dependence of the final design on the constraints imposed. Further investigation should consider the possibility of redefining such areas which do not contribute to increase the performance of the structure in spite of being constrained.

3.5E+07 0

10

20

30

40

50

60

Iteration

6. Conclusions

Fig. 18. Evolution of the strain energy function.

be expected that these changes would increase the performance of the structure. However, the evolution of the objective function shows that from iteration 55 fur-

This paper has presented an initial algorithm for 3D structural optimisation in which the structural analysis is boundary element-based and the geometry is parameterised with NURBS. Control points, which define the NURBS, migrate evolving to the optimum design.

Fig. 19. von Mises stress contour plots.

E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1917–1929

Geometry control features are inserted, such as treatment of corners, control point mapping and restart procedures. Corners and shape edges are specially handled by using smoothing and corner effect algorithms, both implemented to ensure smoothness of the surfaces. The complexity of dealing with surfaces has restricted the automatic geometry control and in this initial study some parts of the algorithm have been implemented as a manual procedure. A defined algorithm has been strictly adhered to in all manual modifications to the model. Restarts are inserted to perform addition and removal of control points as well as topology changes. These restarts imply redefinition of the surfaces. Computational applications of the algorithm have been presented proving the algorithm for shape optimisation and showing preliminary results for topology optimisation.

[3]

[4]

[5]

[6]

[7] [8]

Acknowledgments [9]

The authors thank Prof. G.P. Steven for valuable discussions. The first author wishes to acknowledge the financial support for this research work provided by Durham Research studentship Isabel Fleck awards, Doreen Bretherton studentship.

[10] [11] [12]

References [1] Cervera E, Trevelyan J. Evolutionary structural optimisation based on boundary representation of NURBS. Part I: 2D algorithms. Comput Struct, in press, doi:10.1016/ j.compstruc.2005.02.016. [2] Jacobsen JB, Olhoff N, Rønholt E. Generalized shape optimization of three-dimensional structures using materi-

[13]

[14]

1929

als with optimum microstructures. Mech Mater 1998;28: 207–25. Kane JH, Zhao G, Wang H, Guru Prasad K. Boundary formulations for three-dimensional continuum structural shape sensitivity analysis. J Appl Mech 1992;59:827– 34. Haug EJ, Choi KK, Komkov V. Design sensitivity analysis of structural systems. New York, NY: Academic Press; 1986. Young VO, Querin OM, Steven GP. 3D and multiple load case bi-directional evolutionary structural optimization (BESO). Struct Optim 1999;18:183–92. Annicchiarico W, Cerrolaza M. Structural shape optimization 3D finite-element models based on genetic algorithms and geometric modeling. Finite Elem Anal Des 2001;37:403–15. Haftka RT, Grandhi RV. Structural shape optimization, a survey. Comput Meth Appl Mech Eng 1986;57(1):91–106. Kodiyalam S, Vanderplaats GN. Constructive solid geometry approach to three-dimensional structural shape optimization. AIAA J 1992;30(5):1408. Schramm U, Pilkey WD. The application of rational Bsplines for shape optimization. In: Proceedings of the First World Congress of Structural Multidisciplinary Optimization (WCSMO-1), Germany, 1995. Xie YM, Steven GP. A simple evolutionary procedure for structural optimization. Comput Struct 1993;49(5):885–96. Piegl L, Tiller W. The NURBS book. 2nd ed. SpringerVerlag; 1997. BEASY Users Manual. Computational Mechanics BEASY Ltd Ashurst Lodge, Ashurst, Southampton, Hampshire, SO40 7AA, England, 2002. Cervera E. Evolutionary structural optimisation based on boundary representation of B-spline geometry. PhD thesis, University of Durham, Durham, UK, 2003. Bathe K. Finite element procedures in engineering analysis. Englewood Cliffs, New Jersey: Prentice-Hall, Inc.; 1982.