- Email: [email protected]
Evolutionary structural optimisation based on boundary representation of NURBS. Part I: 2D algorithms
Computers and Structures 83 (2005) 1902–1916 www.elsevier.com/locate/compstruc
Evolutionary structural optimisation based on boundary representation of NURBS. Part I: 2D 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 evolutionary structural optimisation (ESO) approach based on the boundary element method. Non-uniform rational B-splines (NURBS) are used to define the geometry of the component and, since the shape of these B-splines is governed by a set of control points, use can be made of the locations of these control points as design variables. The developed algorithm creates NURBS-based internal cavities to accomplish topology changes. The optimum topologies evolve allowing cavities to merge between each other and to their closest outer boundary. Two-dimensional structural optimisation is investigated in detail exploring single and multiple load case elastic problems. 2005 Elsevier Ltd. All rights reserved. Keywords: Evolutionary structural optimisation; Boundary element method; NURBS
1. Introduction In shape optimisation, the topology of the structure is fixed and only the shape of the boundary can change. Extensive work has been done in shape optimisation [1] using numerical methods, such as the finite element method (FEM) [2,3] and the boundary element method (BEM) [4], for structural analysis. Of special interest to this work has been the application of BEM when design sensitivity analysis is employed in the optimisation algorithm [4]. This is due to the high accuracy of the BEA results on the boundary. Indeed, since the method deals with integrals over the boundary, only the boundary needs to be discretised which is a clear advantage for * Corresponding author. Tel.: +44 191 334 2522; fax: +44 191 334 2390. E-mail address: [email protected] (J. Trevelyan).
remeshing purposes. Heuristic methods have been also applied to shape optimisation. These methods are known as gradientless methods since they require no sensitivity calculations. Schnack and Spo¨rl [5] introduced a method for the reduction of stress concentration on boundaries through the gradual removal of low stress material. In this fashion, the biological growth method [6] is based on the hypothesis that in nature structures such as trees have a constant stress distribution over the boundaries. In the field of structural optimisation, topology optimisation refers to optimal design problems in which the topology of the structure is allowed to change in order to improve the performance of the structure. For that reason, this optimisation class is regarded as one of the most challenging optimisation problems. Much research has been devoted to topology optimisation over the last 10 years; this has included methods based on microstruc-
0045-7949/$ - see front matter 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2005.02.016
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
tures such as the homogenisation method [7] and the solid isotropic microstructure with penalty (SIMP) method [8], and methods that deal with elements on a macro basis approach such as the soft-kill option (SKO) [9], the evolutionary structural optimisation (ESO) method [10] and the bubble method [11]. Hassani and Hinton [12] present a useful review of topology optimisation methods. More recently, methods such as genetic algorithms (GAs) [13] have been applied to topology optimisation since they are particularly robust in finding global optima. However, this can come at a high computational cost. Also, the concepts of robust and reliable design [14] are of growing interest in topology optimisation since they consider input variations and uncertainties in the design and manufacturing process. From the different topology optimisation techniques presented, the evolutionary structural optimisation method ESO [10] is chosen as the basis for the current research. Most of the work carried out on ESO is based on the FEM. The classical ESO [15] approach is based on the idea of removing inefficient material from an initially oversized domain. The removal process is carried out by deleting regions occupied by elements with low stresses. By repeating this process and removing small amounts of material at each stage, the topology for the structure gradually evolves to a more efficient structure. Following this basic approach, there have since been a number of modifications and refinements such as BESO [16] where not only are elements removed but are also added in high stressed areas. Recent applications of ESO [17] cover a wide range of physical situations by considering element modification sensitivity terms. These sensitivities are used to drive the removal and addition process in order to achieve a minimum (or maximum) of the objective function. Nevertheless, ESO has some drawbacks and weaknesses, notably those related to the rejection criteria, as has been reported by Zhou and Rozvany [18], those due to mesh dependency and those due to the tendency to converge to meshes exhibiting an undesirable single corner contact [19]. Moreover, the final solution, due to the nature of the mesh, can result in jagged edges. It is clear that this is an undesirable situation, since the accuracy of the FE results is in doubt. In addition, post-processing, such as spline construction, must be carried out to smooth the boundary for manufacturing reasons [20]. Similar ideas are presented by Maute and Ramm [21] and Hammer and Olhoff [22] in the material topology optimisation approach, in which the geometry is smoothed using splines based on density distributions. Another disadvantage, common in fixed grid FE-based structural optimisation methods, is the presence of checkerboard patterns [23]. These refer to the phenomenon of alternating presence of void and solid elements ordered in a checkerboard fashion. Such patterns result in structural analysis inaccuracies which complicate the
1903
interpretation of optimal material distribution and geometry extraction. Checkerboard formation is related to the use of low-order elements in the finite element approximations. The use of higher order elements can reduce this effect, but at the cost of increasing the computational time required. Alternatively, to improve the estimation quality of elemental sensitivity in low-order elements, a smoothing algorithm [23] can be implemented. On the other hand, in the ESO-based intelligent cavity creation (ICC) algorithm [24] checkerboard patterns (with numerous cavities) can be eliminated through controlling the number and scale of structural cavities in the final topology. Alternatively, a perimeter control technique [25] can be incorporated into BESO to overcome checkerboard patterns. This technique has also been shown to reduce the dependency of the converged optimum on the initial finite element discretisation. Finally, an important task pointed out by Querin et al. [16] is related to the solution time, that is, the necessity to make the ESO process faster so that the designer is able to obtain the optimum design within a few seconds of describing the environment. Research in this topic has been carried out combining ESO and the fixed grid (FG) method [24], instead of classical FE, to improve the performance. In order to overcome these drawbacks relating largely to domain-based discretisation, the boundary element method (BEM) [26] may be implemented as the tool to carry out structural analysis. In addition, in the current work the geometry is described using the nonuniform rational B-spline curves (NURBS) [27]. The optimisation process is fully integrated within in-house boundary element software [28] allowing a straightforward and rapid communication between the analysis software and optimisation code since they share the same database. The advantages of the proposed algorithm over FEM based ESO procedures are, therefore, the easy remeshing and smoothness of geometry that provides for a smoothness in the stress solution. The use of boundary elements enables further computational efficiency gains beyond a smoothed FE approach (e.g. [21,22]) since many computations may be saved by reusing influence coefficients relating to the boundaries that remain unchanged in successive iterations. The use of the BEM, however, restricts the class of problems to those to which are suitable for BEM treatment. The present work is limited to linear elastostatics, but this might be extended through, for example, dual reciprocity methods, to non-linear problems.
2. Algorithm The developed algorithm considers the BEM to carry out the structural analysis of the component under study. The optimisation approach is stress-based selecting low
1904
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
and high-stressed regions as potential areas to be modified. The geometry of the component is described using NURBS. These parametric curves allow free-form representation with total geometry control over the curve. This feature is very interesting since, a priori, the final design is not known. NURBS are defined by control points which become the design variables of the problem. Another important property of these curves is the locality, that is, changes in one area do not propagate throughout the whole B-spline. The complete optimisation algorithm [29] is presented in detail in the following sections of this paper. First, the main steps of the basic algorithm are summarised as follows: Step 1: The geometry of the structure is defined (quite possibly as a simple rectangle bounding a given design space). Loads and constraints are applied. Step 2: A boundary element analysis (BEA) is carried out. Step 3: The material removal process is performed by selecting the least stressed nodes within the boundary mesh and effectively moving the control points nearest to those nodes. Following a similar process, material addition is carried out if a node is found with a stress higher than the yield stress or a certain maximum stress criterion. The control points nearest to the high-stressed nodes would be moved. However, movement takes place in the opposite direction to that in the removal process. Step 4: To complete the topology optimisation, the inside regions are also considered and internal holes are created in the inner low-stressed areas. As a natural evolution of the process, holes can merge between each other or into the outer boundary. Step 5: Finally, such a procedure is repeated, from Step 2, until the evolution of the objective function shows no improvements. Notice that since, during the optimisation process, the control points are moved and in consequence the shape changes, the boundary element mesh is updated every iteration. The algorithms that are presented in this work have been developed heuristically and, although use is made of parameters that have been used successfully by authors in finite element implementations of ESO, the values of these parameters have been varied according to the results of numerous numerical experiments. These experiments are not described in detail here. Further work is in progress to determine the sensitivity of the process to changes in the parametric values.
2.1. Geometry definition According to the applied constraints, loads or any other design requirements, the boundary domain can be divided into three different types of curves; i.e. design domain, non-design domain and symmetry lines. Lines that can change freely along the process are identified as design domain, whereas those lines that cannot change due to constraints or design restrictions are identified as non-design domain. Symmetry lines can be regarded as an intermediate step between design and non-design domain. They can be modified by changing their length but only if the adjoining line is a design domain line; furthermore, their changes are as a result of changes in the design domain line. 2.2. Boundary element model In the boundary element method (BEM) the boundary integral equations are approximated by a set of discretised integral equations. As a result, the boundary surface is divided into elements, thus the response is given at the nodal points associated with the elements. The main advantage of this method relates to the mesh, since only the surface of the structure needs to be discretised, though their accuracy of boundary stress solution is of course advantageous to this work. The results inside the structure are calculated at an arbitrary number of internal points. In the software used [28], these points are randomly distributed throughout the interior domain. 2.3. Removal and addition of material 2.3.1. Identifying inefficient areas Following the elastostatic boundary element analysis, for any iterative step, material is removed from areas of low stress and added in areas of high stress. Material can be removed from the structure if any node p satisfies ð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). These ratios are conceptually similar to the ones used in FE-ESO [10]. In the customary ESO fashion, if a steady state is reached in which no nodes, or only a few nodes, can satisfy Eq. (1) then RR is increased by the evolutionary rate for removal, ERR, as follows:
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
RRi ¼ RRi1 þ ERR
ð3Þ
where i is the current iteration. Similarly, if only a few nodes can satisfy Eq. (2) then AR is decreased by the evolutionary rate for addition, ERA, ARi ¼ ARi1 ERA
ð4Þ
Typical values for suitable ratios are RR0 = 0.01, ERR = 0.01, AR0 = 0.99, ERA = 0. These values are determined from numerical experience. It should be noted that the von Mises stress is the criterion generally used in this work. Nevertheless, other stress criteria can be considered. For example, the principal stresses r1 and r2 take into account the differences between tension and compression, which may be important in cases having two materials; for example, steel and concrete. 2.3.2. Multi-load case problems In most engineering disciplines multiple load case problems are frequent. For such situations, the final solution depends on the scheme adopted to compute the effect of each load case. The criterion explored in this work is the logical AND/OR algorithm [30]. Here, a node is selected for the removal process only if is lowstressed in all load cases, and it is selected for the addition process if is high-stressed in any of the load cases. Eqs. (1) and (2) refer to a single load case problem. In the case of multiple load cases these equations are slightly changed, thus, for each node p and each load case j, a ratio between stresses is calculated as follows: rp;j Rp;j ¼ ð5Þ rmax;j Computing the above equation for all the load cases would give a set of ratios associated to each node. For example the node p has associated N ratios, where N is the number of load cases. Rp ¼ fRp;0 ; Rp;1 ; . . . ; Rp;N g
ð6Þ
Thus, applying the logical AND/OR, material can be removed from the structure if all the ratios {Rp} are lowstressed ðRp;0 6 RRÞ AND ðRp;1 6 RRÞ . . . AND ð7Þ
ðRp;N 6 RRÞ In the same way, material can be added if ðRp;0 P ARÞ OR
ðRp;1 P ARÞ . . . OR
ðRp;N P ARÞ ð8Þ
where RR and AR are the same optimisation parameters as the ones related to Eqs. (1) and (2).
1905
control points is associated to each specific area to be moved. Every set consists of the three nearest control points to this area. The justification for this choice of three points comes from numerical tests with different numbers of control points. For sets of less than three control points it has been found that the distortion of the curve is very localised causing extreme curvature, while for sets of more than three control points the distortion is propagated to a wider region of the curve and therefore, losing control of the areas to be removed/ added. Each set of control points is moved a distance related to the following parameters: • Length of the least/most stressed element, Le. • Distances of the three control points from the least/ most stressed boundary element. These distances are denoted a, b, c, respectively. • A factor related to the stress situation within the structure at the current iteration, which is called the removal factor (RF) if removing material, and the addition factor (AF) if adding it. This factor allows larger geometric changes in early iterations and smaller changes when fine-tuning is required in the later iterations. The factors RF and AF allow us to gain in computational efficiency without losing the accuracy of the final solution. For example, in Fig. 1(a) the control point (P2) situated at a distance a from the node of lowest stress is moved a distance calculated as follows: 1 0 B Movement of P 2 ¼ Le @
1 C1 A RF 1 1 1 a þ þ a b c
ð9Þ
2.3.4. Direction of movement The direction of movement for each set is perpendicular to its nearest least/most stressed boundary element, i.e. follows the normal vector, n, to the element as shown in Fig. 1(a). The only difference between removal and addition is the inwards movement when removing material and outwards movement when adding it. In the case of the control point at the intersection of two B-splines, as depicted in Fig. 1(b), the direction of movement follows an averaged normal vector. This vector, n, is calculated as the resultant of the two normal vectors, nA and nB , to the B-spline curves involved. 2.4. Creation of holes
2.3.3. Distance to move In any iteration, material is either removed or added to the structure by changing its boundary definition. Since control points define the NURBS curves, a set of
Topology changes are carried out by creating internal cavities or holes in low-stressed areas. These holes are shaped according to the internal von Mises (or any other
1906
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
0
e
1
P2
n
Least stressed node
a
Least stressed element
P4
P Pl
nA
b
P0
P0
B-spline B
P3
c
nB
n = n A + nB
B-spline A
P1 (a)
(b)
Fig. 1. Direction of movement when removing material: (a) one B-spline; (b) two B-splines.
selected criterion) stress contours. The free-shapes are created using NURBS curves to define their geo-metry. In a BEA the results at points contained within the boundary of the object are calculated at internal points. Recalling Section 2.2, in the software used in the current study, these internal points are scattered at random (though concentrated near stress concentrations and their distribution smoothed using a Laplacian algorithm). It is the results at these points, and not at the boundary nodes, that form the starting point of the algorithm for the creation of holes in the domain. Nevertheless, the basis to identify the inner low-stressed areas is similar to the one used for the outside boundary (Eq. (1)). The process of creation of holes is summarised in the following steps:
Step 2:
Step 3:
Step 4:
Step 5:
Step 1: Check for any internal point that satisfies Eq. (10). If so, add this internal point to a stack denoted stackIP ð10Þ
rIP 6 RRrmax
where rIP is the internal point von Mises stress or any other selected criterion, rmax is the maximum von Mises stress or any other se-
Step 6:
IP8
IP14
IP14 IP7
IP8
IP1 IP2
P 10
P20 IP0
IP4
IP9
IP6 IP3
IP13
lected criterion and RR is the removal ratio used in Eq. (1). The stackIP containing the selected internal points is sorted in ascending von Mises stress order. The first point of the stackIP is identified as the minimum and removed from the stackIP to another stack called temporarystackIP. This minimum is taken as the initial central point to be enclosed by the new contour. Internal points surrounding this minimum and having a von Mises stress lower than a threshold (generally related to Eq. (10)) are also moved from the stackIP to the temporary-stackIP. A convex polygon is created from the internal points in the temporarystackIP, as displayed in Fig. 2(a). This polygon is split up into two open polygons and the internal points are taken as the control points. Therefore, these polygons together with the control points define the two NURBS curves that create the new cavity (Fig. 2(b)). If the remaining stackIP is not empty, the process is repeated from point 2 until there are no more internal points in the stackIP.
IP13
NURBS 0
P01 P30
IP10 IP11
(a)
Least stressed internal point IP0
P00
NUR NURBS1
IP5 IP12
P41
P 11
P31
P21 IP11
(b)
Fig. 2. (a) Creation of holes from internal points; (b) new hole consists of two NURBS.
IP9
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
Holes, which are described by two NURBS curves per hole, have the same form of boundary representation as the outside boundary. Consequently, they have similar (but not identical) properties and behaviour. Finally, note that if the von Mises stress in most of the internal points surrounding the minimum is higher than the threshold, then the temporarystackIP consists of only a few internal points and in cases such as this no holes are inserted. 2.4.1. Distance to move control points defining holes To calculate the distance to move the control points in internal holes, a similar procedure to the outer boundary is followed. However, the set of control points associated to each specific area to be moved consists of only one control point, the closest point to this area. The reason for moving exclusively the closest control point is that it is common for stresses around cavities to oscillate between high and low stress regions. A good example is the classical case of a circular hole in a plate under uniaxial tension, in which the tangential stress around the hole
1907
varies in a well known cosine form, containing both tensile and compressive stresses and passing through zero in four locations around the circumference. The use of a three control point set would cause geometric changes to be too widespread and would greatly slow down, or probably entirely prevent, convergence to an optimum. Thus, the control point is moved a distance related to the following parameters: • length of the least/most stressed element, Le; • distance, a, of the control point from the low/high stressed node; • the removal factor (RF) if removing material, or the addition factor (AF) if adding it; • thickness of the component, t. For example, in Fig. 3(a) the closest control point to the node of lowest stress is moved a quantity calculated as follows: 1 ð11Þ Movement of P 2 ¼ Le tRF a
Fig. 3. Direction of movement: (a) isolated hole; (b) hole with a close boundary.
P4 A
P3A
P7C
P6C
Hole A and Hole B merge
P4B P5B
P2 A
P5A P3B
P6A
P2B P1B
P8C Hole C
P0A P8 A
dist |P3B - P6A| < min low stress area
P4C
P9C
P7 A
Hole B P0B
Hole A P1A
P3C P0C
P1C
Fig. 4. Merging two holes (A, B) into a new one (C).
P2C
1908
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
Adjacent holes merge only if they satisfy these two conditions. As shown in Fig. 4, merging the holes implies the creation of a new one. Consequently, this follows a similar procedure to the hole creation process detailed in Section 2.4. In this case, stackIP consists of all the control points related to the holes to be merged.
2.4.2. Direction of movement in holes Generally, the direction of movement for a control point in a hole is perpendicular to its nearest least/most stressed boundary element (Fig. 3(a)). Nevertheless, in the case of a close outer boundary (Fig. 3(b)) the direction of movement is determined by this borderline, i.e. parallel to the boundary that is closest to the selected node.
2.4.3.2. Merging hole to outside boundary. The criterion for merging a hole with the outside boundary is similar to the one for merging two holes (Section 2.4.3.1), that is, minimum distance and under stressed area in between. Fig. 5 depicts how a hole is merged to its closest outer boundary. The process implies a change of the boundary definition and also the deletion of the hole. The extreme
2.4.3. Merging holes It appears that the algorithm works towards an optimum by producing multiple holes which later merge to form larger cavities. Moreover, it is found that, without a hole merging algorithm, the holes approach each other to produce multiple thin filaments of material which are not efficient as a use of material and present unnecessary difficulties in manufacture. 2.4.3.1. Merging two holes. This situation occurs when the cavities evolve approaching too close to each other. The criterion for merging holes is based on the following two premises:
P0B
P 0B
Hole P1 B
P 4H P2B
• Minimum distance. If the distance between two holes is less than a minimum, dmin, those holes are likely to be merged. The distances are calculated between control points from each hole. The largest element length, from the elements in the immediate vicinity, is taken as dmin. • Under stressed area in between. Internal points situated in the region between the holes satisfy Eq. (10).
P1 C(u)
P1H
P3B
P0 H
P4B
New Boundary
P6B
low stress area Fig. 5. Merging hole (hole) into boundary (boundary).
P1
Inserted control point P’
C’(u)
P2
P0
C(u) |P1-P2|i > k|P1-P2|0
(a)
(b) Fig. 6. Insertion of control point P 0 .
P2
P2
P4
C’(u)
C(u)
P3
P4
P1 C(u)
P’ P0
P 4B
dist |P3B – P5H| < min
P2
P1
P3B
Hole and Boundary merge P5B
P5 B
d0 = |P1-P2|0
P0
P2B
P2H
P5 H
Boundary
P1 B
P 3H
P’ P0
P3
|P’-P3| < k’ |P2-P’|
Removed control point (a)
(b) 0
Fig. 7. Removal of control point P .
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
The number of control points defining each NURBS curve is not fixed. On the contrary, control points are automatically either inserted or deleted as required throughout the process. This capability has been implemented in order to control geometry changes and keep smoothness. Another important feature related to the geometry control process is the mesh subdivision. This subdivision is performed in areas in which the NURBS curve distortion has caused excessive curvature of an element. 2.5.1. Insertion of control points The insertion of a new control point between two existing control points is done when these control points are moved away from each other more than a constant factor, k, times the initial distance between them jP j P j1 ji > kd 0
ð12Þ
where i is the iteration number, d0 = jPj Pj1j0 and generally k P 1.5. The control point insertion algorithm is implemented to avoid areas whose defining control points are placed too far away, causing loss of localised control over the B-spline curve. This effect is illustrated in Fig. 6. Fig. 6(a) depicts the original curve C(u) with d0 = P1 P2. As the optimisation evolves, the original curve C(u) changes according to the movement of its control points. In this example, P1 is moving away from the curve compared to the other control points. Fig. 6(b) shows the curve C(u) at iteration i for which Eq. (12) is satisfied. Thus, C(u) is updated to C 0 (u) inserting a control point P 0 to the control point set. The new control point is placed at the mid-point of the line joining both control points P1 and P2.
where k 0 < 1, jPj Pj1j is the length of the actual line checked, jPj1 Pj2j and jPj+1 Pjj are the lengths of the previous and next control polygon lines, respectively. Thus, if in Eq. (13), the inequality on the left is satisfied, then Pj1 would be removed. On the other hand, if it is the inequality on the right that is satisfied then Pj is the control point to be removed. For example, Fig. 7(a) displays an original curve C(u) and the distribution of control points defining the curve. Fig. 7(b) shows the curve C(u) some iterations after, when P 0 and P2 have been moved. At this stage P 0 has approached too close to P3 and the distance between the two control points satisfies Eq. (13) since jP3 P 0 j < k 0 jP2 P 0 j. Finally, C 0 (u) is created by removing control point P 0 and so modifying the control polygon definition.
σ y = 10 N/mm 2 120 mm
2.5. Geometry control
2.5.2. Removal of control points The removal of a control point is performed when a set of control points are moved to locations too close together. This gives rise to regions of excessive and unpredictable curvature of the boundary. The criterion to remove the point from this set is based on comparison of the different lengths of the lines of the control polygon P j P j1 < k 0 P j1 P j2 OR P j P j1 < k 0 P jþ1 P j ð13Þ
σ x = 20N/mm 2 B
30 mm
(maximum and minimum x y coordinates) control points defining the hole are identified and added to the set of control points describing the boundary. In addition, control points from the boundary which are between these extreme points are removed.
1909
y x
A
30 mm
120 mm
Fig. 8. Problem definition for a quarter of the plate.
Fig. 9. Initial and final (a) element mesh; (b) control point distribution.
1910
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
2.5.3. Mesh subdivision An essential element of a fast optimisation algorithm is a reliable automatic mesher to discretise each model. Such a mesher is available in the adopted software, but it has been found important to enhance its conventional automesh algorithm to handle the unusual level of distortion of B-splines in the later iterations. In areas of high curvature, an individual element may need to be subdivided in order to maintain suitable solution accuracy. An element is subdivided if its overall length exceeds certain factor of the distance between the ends of this element as stated in Eq. (14) Le > mpj pj1 ð14Þ
criterion. This stopping criterion can be quantified in the following normalised form: iþ1 f f i 6 104 e¼ ð15Þ fi where fi is the value of the monitored quantity at iteration i and fi+1 is the value of the monitored quantity one iteration after. 2.6.1. Stress levelling The failure of structures under service frequently takes place in areas of locally high stresses. Therefore, it is crucial to avoid stress peaks in order not only to
where Le is the length of the element, e, under study. pj and pj1 are the first and last node of the element and jpj pj1j is the distance between both nodes and generally, m P 1.05.
B Numerical optimum b/a = 0.5521 b
2.6. Stopping criteria Different stopping criteria can be adopted depending on the target to achieve. In this paper two different criteria have been considered, one being based on specific stiffness (which has an equivalence with strain energy) and one on achieving a uniformity of stress over the boundary. Monitoring of the required criterion towards, and ultimately beyond, an optimum provides a stopping
Analytical optimum b/a = 0.5 y A
x a
Fig. 10. Comparison between the theoretical and analytical solution.
Fig. 11. (a) Initial and (b) final von Mises stress contour plot.
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
prevent localised yielding but also to increase the fatigue life of the structure and avoid crack initiation. Stress levelling aims to make uniform the stress along a detailed part of the boundary (C). Therefore, the stress levelling objective deals with minimising the following function: Z fC ¼ ðrp rref Þ2 dC ð16Þ C
structure. However, simply maximising stiffness is not very useful in most engineering contexts since this will result in an optimum in which the design space is completely filled with material. Thus, dividing the stiffness (K) by volume (V) a function termed specific stiffness [17] can be defined as fK ¼
where rp is the node stress and rref is a reference stress. 2.6.2. Strain energy criterion A common objective in engineering is to obtain a structure stiff enough to produce a good performance under certain loads and constraints. Such an optimisation problem consists of maximising the stiffness in the
K V
Von Mises Stress (N/mm2)
Initial design Optimum numerical Optimum analytical
B
y
x
A
10
20
30
40 50 60 Angle (degrees)
ð18Þ
The simplest method of computing strain energy is to use the boundary integral Z 1 U¼ Tu dC ð19Þ C 2
θ
0
ð17Þ
Therefore, the overall specific stiffness is maximised which implies that the stiffness is maximised while the volume is reduced. An equivalent expression can be written in terms of strain energy (U). In this case the objective function adopted is related to the minimisation of the specific strain energy fU as follows: fU ¼ UV
90 80 70 60 50 40 30 20 10 0
1911
70
80
90
Fig. 12. Distribution of von Mises stress along the hole. Fig. 14. Initial and final boundary element mesh.
Fig. 13. Problem definition for a short cantilever beam.
Fig. 15. Initial and final control point distribution.
Fig. 16. Optimum designs: (a) Rozvany et al. [31]; (b) Chu et al. [32]; (c) boundary-ESO.
1912
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
Fig. 17. Von Mises stress contour plot.
where T are the tractions over the boundary and u the displacements over the part of the boundary C where the tractions are applied.
1.0E+04
In this section, several examples are presented in order to prove the developed algorithm. The following isotropic material properties are assumed: YoungÕs modulus E = 210,000 N/mm2, PoissonÕs ratio = 0.3 and an arbitrary thickness t = 1 mm. The yield stress for this material is 280 N/mm2. Plane stress conditions are also assumed.
f = U*V(Nmm4)
9.0E+03
3. Examples
8.0E+03 7.0E+03 6.0E+03 5.0E+03
0
5
10
15 20 25 Iteration number
30
35
Fig. 18. Evolution of the strain energy function.
40
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
3.1. Hole in a biaxial stress field This classical example addresses the shape optimisation of a hole in a large plate under a biaxial stress field. Due to the symmetry, only a quarter of the plate is modelled. Fig. 8 shows the initial shape and boundary conditions. The plate is subjected to stress rx = 20 N/mm2 and ry = 10 N/mm2. The hole is modelled by a quadratic NURBS curve (line B–A). The hole is monitored for stress levelling (Eq. (16)) and in this case, the maximum von Mises stress rmax is taken as rref. The theoretical solution [3] is an ellipse with an aspect ratio b/a equivalent to the loading ratio ry/rx, i.e. 1/2, with the stresses uniformly distributed around the hole. This represents the optimum shape condition, in which all the material at the edge of the hole is at the same maximum permissible stress. The problem was also studied by [4] using the BEM and design sensitivity analysis. The optimisation parameters are set to RR0 = 0.01, ERR = 0.01, AR = 0.99, ERA = 0, RF = 0.05 and AF = 0.025. RF and AF are taken to be constant and of small value in order to keep smoothness and thereby do not give rise to artificial stress peaks over the boundary. Initially, the element model contains 44 quadratic elements, a number that increases to only 45 in the optimised design. Fig. 9(a) shows the initial and final boundary mesh in the area of the hole. There is one quadratic NURBS curve defining the central hole. This curve consists initially of six control points. At the end of the process this number has increased to 9. See Fig. 9(b) for the initial and final distribution of the control points. The final design for the hole is reached after 60 iterations in a total CPU time of 14.162 s on a Pentium 4 (2 GHz). Fig. 10 compares the numerical optimum shape with the theory. In order to compare both results, the analytical optimum is calculated taking the major axis a = 47.8 mm. Since a = 2b, for this problem the minor axis is b = 23.9 mm. Thus, the ellipse (analytical optimum) is calculated following the equation for an ellipse being centred at (0, 0): x2 y 2 x2 y2 þ 2¼1! þ ¼1 2 2 a b ð47.8Þ ð23.9Þ2
1913
compared to the analytical 1.5 [3] and other numerical results 1.59 [4]. A reason for this can be that the theoretical results are calculated assuming an infinite plane whereas in our case the geometry is described by finite segments and the effects of the loading and displacement constraints cannot be avoided. It may also be related to the choice of objective function where rref is the maximum von Mises stress in preference to an average value. Fig. 12 shows the distribution of the von Mises stress along the hole for the initial and final designs. In the arbitrary initial design there is a stress peak of 80 N/mm2 that occurs in a localised area at h = 0. Comparing this result to the (numerical) optimum design; in the latter, the stress concentration has been effectively removed and the stress levels are clearly uniform around the hole. The optimum analytical solution [3] rh = (1 + ry/rx)rx, has a constant value rh = 30 N/mm2. Thus, there is a very good general agreement between the analytical and the numerical optimum. 3.2. Short cantilever beam A common benchmark problem is the cantilever beam with a central load and an aspect ratio of 1.6. The definition of the initial problem is shown in Fig. 13. This cantilever of dimensions 160 mm · 100 mm is fixed at the top and bottom corners of the left hand side; whereas it is subjected to a vertical load at the centre of the free end. The load has been chosen arbitrarily to have value 10 N. This value can be regarded to be very small compared to the failure stress; however it produces an appropriate stress distribution to study removal of
ð20Þ
The numerical results show a slight deviation from the theoretical solution. The ratio of the minor axis of the elliptic hole to the major is b/a = 25.29 mm/ 45.81 mm 0.5521 (analytical optimum b/a = 0.5). Comparing the initial (Fig. 11(a)) and final (Fig. 11(b)) stress distribution, as expected, the stresses around the curve become more uniform when approaching the final solution. The ratio between the maximum von Mises stress at the edge of the hole and the larger of the two biaxial components is 1.691 (rmax/rx = 33.828 N/mm2/20 N/mm2 = 1.691), which is high if
Fig. 19. Initial and final design for a bridge under moving load.
1914
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
material as well as the creation and evolution of holes. The optimisation parameters are AF = 1.0, constant AR = 0.98, ERR = 0.01, RR varies from 0.01 to 0.2, and RF changes from 0.8 to 0.5 according to the stress levels within the structure. The boundary domain is initially discretised into 37 quadratic elements (Fig. 14(a)). At the end of the process this number has increased to 165 elements (Fig. 14(b)).
The number and positions of the control points for the initial and final topology design are depicted in the Fig. 15. The initial five non-constrained straight lines are converted into five quadratic NURBS curves. In the final design this number has increased, with the creation of internal holes, to 11 NURBS. The analytical optimum for this case consists of an infinite number of infinitesimally thin members. Fig.
Fig. 20. Von Mises stress contour plots for LC0, LC1 and LC2.
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
16 compares the optimum design obtained with the approach (boundary-ESO) to those approximations using the optimal layout theory by Rozvany et al. [31] and to a solution obtained by Chu et al. [32] using FEESO. Fig. 16(b) shows the ESO solution obtained under stiffness constraints. The optimum design obtained with the current approach (Fig. 16(c)), shows a good agreement with the FE-ESO. Fig. 17 shows, in terms of von Mises stress contours, how the structure evolves towards the final solution. As the optimisation progresses, small cavities are merged into larger cavities. Also cavities that are very close together or to the outer boundary may merge. Specific stiffness is monitored to form a stopping criterion, and the optimum topology design is reached after 38 iterations. The volume reduction is approximately 63% (V/V0 0.37) from the initial design. Although the overall process shows a clear reduction of volume (removal of material) in the structure, it should be noted that material is also added if any stress concentration appears as a result of the movement of the control points. The final design is obtained in a CPU time of 5 min and 27 s (Pentium 4 (2 GHz)). The evolution of the objective function (Eq. (18)) is displayed in Fig. 18. As can be observed, the function decreases smoothly until iteration 16. At this iteration there is an increment in the function due to the insertion of the first hole. A few iterations later, at iteration 21, two new holes are created. The effect of the creation of these latter holes is not as dramatic as the first one because the size of the holes is much smaller. As the process evolves another two holes are created at iteration 27. These last cavities are very small compared to the existing ones and eventually are merged to the larger holes. The final topology is achieved when no further improvements (Eq. (15)) are shown in the objective function which shows a 28% reduction compared to the initial situation. 3.3. Bridge under a moving load This final example tackles the problem of a bridge structure under multiple loading conditions. The initial design for the bridge is shown in Fig. 19(a). The dimensions of the bridge are 160 mm · 50 mm. The bridge is supported at four points in the bottom surface. A point load of 100 N is moving along the top surface and this effect is approximated by five load cases. The logical AND/OR scheme represented by Eqs. (7) and (8) is adopted. Point boundary loads and constraints are difficult to apply in the BEM, and consequently the loads and constraints in this solution are applied over sections of boundary of finite length. The area over which the constraints are applied are evident in the converged optimal geometry.
1915
The top surface is considered to be in the non-design domain as well as the four supports. Moreover, there is a limit for the maximum vertical displacement of 0.02 mm. The rest of the domain is free to change and defined by 5 NURBS curves. The optimum design, shown in Fig. 19(b) that satisfies this problem is found after 42 iterations with RR varying from 0.01 to 0.27. Fig. 20 shows the von Mises stress contours, at the initial and final iteration, for load cases LC0, LC1, LC2. LC3, LC4 are symmetric to LC1, LC0, respectively. As expected, the final solution under multiple load cases is not necessarily a fully stressed design [20]. The ratio of the final volume to the initial volume is 0.54. The maximum vertical displacement (umax) for any load case is below the limit constraint of 0.02 mm. Certainly, these results would have varied if the constraints over the maximum vertical displacement for each load case had been considered differently.
4. Conclusions This paper presents a new algorithm based on the evolutionary structural optimisation (ESO) method. The algorithm exploits the concept that the optimum topology evolves by slow removal and addition of material. The boundary element method is used for the analysis and NURBS curves for describing the changeable lines. The control points defining these curves are selected to be the design variables of the problem. Topology optimisation is accomplished allowing the insertion of internal holes in the inner low-stressed areas of the structure. Holes are also described by NURBS curves and so they have similar behaviour to the outside boundary. To prove the algorithm, several examples have been reproduced showing the effectiveness of the method generating the final topologies as well as rapid solution times. The approach presents the advantage of decreasing the number of design variables compared to FEbased ESO. Also, it has been proved that working directly with the boundary provides a large flexibility in the design. It has been shown how at each iteration the boundary remains smooth without artificial stress concentrations.
Acknowledgments The authors would like to thank Prof. G.P. Steven for inspiring advice and useful discussion. The first author gratefully acknowledges the Durham Research studentship Isabel Fleck awards, Doreen Bretherton studentship and Lindsay scholarship provided during this study.
1916
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
References [1] Haftka RT, Grandhi RV. Structural shape optimization, a survey. Comput Methods Appl Mech Eng 1986;57(1): 91–106. [2] Braibant V, Fleury C. Shape optimal design using Bsplines. Comput Meth Appl Mech Eng 1984;44:247–67. [3] Kristensen ES, Madsen NF. On the optimum shape of fillets in plates subject to multiple in-plane loading cases. Int J Numer Methods Eng 1976;10:1007–19. [4] Tafreshi A, Fenner RT. Design optimization using the boundary element method. J Strain Anal 1991;26(4): 231–40. [5] Schnack E, Spo¨rl U. A mechanical dynamic programming algorithm for structure optimization. Int J Numer Methods Eng 1986;23:1985–2004. [6] Mattheck C, Burkhardt S. A new method of structural shape optimisation based on biological growth. Int J Fatigue 1990;12(3):185–90. [7] Bendsøe MP, Kikuchi N. Generating optimal topologies in structural design using a homogenisation method. Comput Methods Appl Mech Eng 1988;71:197–224. [8] Rozvany GIN, Zhou M, Birker T. Generalized shape optimization without homogenization. Struct Optim 1992;4:250–2. [9] Baumgartner A, Harzheim L, Mattheck C. SKO (soft kill option): the biological way to find an optimum structure topology. Int J Fatigue 1992;14(6):387–93. [10] Xie YM, Steven GP. Evolutionary structural optimization (ESO). Springer; 1997. [11] Eschenauer HA, Kobelev VV, Schumacher A. Bubble method for topology and shape optimisation of structures. Struct Optim 1994;8:42–51. [12] Hassani B, Hinton E. A review of homogenization and topology optimization III—topology optimization using optimality criteria. Comput Struct 1998;69:739–56. [13] Cerrolaza M, Annicchiarico W, Martinez M. Optimization of 2D boundary element models using b-splines and genetic algorithms. Eng Anal Boundary Elem 2000;24:427–40. [14] Frangopol DM, Maute K. Life-cycle reliability-based optimization of civil and aerospace structures. Comput Struct 2003;81:397–410. [15] Xie YM, Steven GP. A simple evolutionary procedure for structural optimization. Comput Struct 1993;49(5):885–96. [16] Querin OM, Xie YM, Steven GP. Evolutionary structural optimization (ESO) using a bidirectional algorithm. Eng Computat 1998;15:1031–48. [17] Steven GP, Qing L, Proos K, Xie YM. The role of physical sensitivity in evolutionary topology design optimisation
[18]
[19]
[20]
[21] [22]
[23]
[24]
[25]
[26] [27] [28]
[29]
[30]
[31] [32]
with multi-criteria and multi-physics. In: Proceedings of the Fifth World Congress on Computational Mechanics (WCCM V). Vienna, Austria, 2002. Zhou M, Rozvany GIN. On the validity of ESO type methods in topology optimization. Struct Optim 2001;21: 80–3. Sigmund O, Petersson J. Numerical instabilities in topology optimisation: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Struct Optim 1998;16(1):68–75. Hinton E, Sienz J. Fully stressed topological design of structures using an evolutionary procedure. Eng Computat 1995;12:229–44. Maute K, Ramm E. Adaptive topology optimisation of shell structures. AIAA J 1997;35:1767–73. Hammer VB, Olhoff N. Topology optimization of continuum structures subjected to pressure loading. Struct Multidiscip Optim 2000;19:85–92. Li Q, Steven GP, Xie YM. A simple checkerboard suppression algorithm for evolutionary structural optimization. Struct Multidiscip Optim 2001;22:230–9. Kim H, Querin OM, Steven GP, Xie YM. A method for varying the number of cavities in an optimized topology using evolutionary structural optimization. Struct Multidiscip Optim 2000;19:140–7. Yang XY, Xie YM, Liu JS, Parks GT, Clarkson PJ. Perimeter control in the bidirectional evolutionary optimization method. Struct Multidiscip Optim 2003;24:430– 40. Trevelyan J. Boundary elements for engineers. Computational Mechanics Publications; 1994. Piegl L, Tiller W. The NURBS book. second ed. SpringerVerlag; 1997. Trevelyan J, Wang P. Interactive re-analysis in mechanical design evolution. Part I. Background and implementation. Comput Struct 2001;79(9):929–38. Cervera E. Evolutionary structural optimisation based on boundary representation of B-spline geometry. PhD Thesis. University of Durham, Durham, UK, 2003. Li Q, Steven GP, Xie YM. On equivalence between stress criterion and stiffness criterion in evolutionary structural optimization. Struct Multidiscip Optim 1999;18:67– 73. Rozvany GIN, Bendsøe MP, Kirsch U. Layout optimization of structures. Appl Mech Rev 1995;48(2):41–119. Chu DN, Xie YM, Hira A, Steven GP. On various aspects of evolutionary structural optimization for problems with stiffness constraints. Finite Elem Anal Design 1997;24(4): 197–212.
Evolutionary structural optimisation based on boundary representation of NURBS. Part I: 2D 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 evolutionary structural optimisation (ESO) approach based on the boundary element method. Non-uniform rational B-splines (NURBS) are used to define the geometry of the component and, since the shape of these B-splines is governed by a set of control points, use can be made of the locations of these control points as design variables. The developed algorithm creates NURBS-based internal cavities to accomplish topology changes. The optimum topologies evolve allowing cavities to merge between each other and to their closest outer boundary. Two-dimensional structural optimisation is investigated in detail exploring single and multiple load case elastic problems. 2005 Elsevier Ltd. All rights reserved. Keywords: Evolutionary structural optimisation; Boundary element method; NURBS
1. Introduction In shape optimisation, the topology of the structure is fixed and only the shape of the boundary can change. Extensive work has been done in shape optimisation [1] using numerical methods, such as the finite element method (FEM) [2,3] and the boundary element method (BEM) [4], for structural analysis. Of special interest to this work has been the application of BEM when design sensitivity analysis is employed in the optimisation algorithm [4]. This is due to the high accuracy of the BEA results on the boundary. Indeed, since the method deals with integrals over the boundary, only the boundary needs to be discretised which is a clear advantage for * Corresponding author. Tel.: +44 191 334 2522; fax: +44 191 334 2390. E-mail address: [email protected] (J. Trevelyan).
remeshing purposes. Heuristic methods have been also applied to shape optimisation. These methods are known as gradientless methods since they require no sensitivity calculations. Schnack and Spo¨rl [5] introduced a method for the reduction of stress concentration on boundaries through the gradual removal of low stress material. In this fashion, the biological growth method [6] is based on the hypothesis that in nature structures such as trees have a constant stress distribution over the boundaries. In the field of structural optimisation, topology optimisation refers to optimal design problems in which the topology of the structure is allowed to change in order to improve the performance of the structure. For that reason, this optimisation class is regarded as one of the most challenging optimisation problems. Much research has been devoted to topology optimisation over the last 10 years; this has included methods based on microstruc-
0045-7949/$ - see front matter 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2005.02.016
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
tures such as the homogenisation method [7] and the solid isotropic microstructure with penalty (SIMP) method [8], and methods that deal with elements on a macro basis approach such as the soft-kill option (SKO) [9], the evolutionary structural optimisation (ESO) method [10] and the bubble method [11]. Hassani and Hinton [12] present a useful review of topology optimisation methods. More recently, methods such as genetic algorithms (GAs) [13] have been applied to topology optimisation since they are particularly robust in finding global optima. However, this can come at a high computational cost. Also, the concepts of robust and reliable design [14] are of growing interest in topology optimisation since they consider input variations and uncertainties in the design and manufacturing process. From the different topology optimisation techniques presented, the evolutionary structural optimisation method ESO [10] is chosen as the basis for the current research. Most of the work carried out on ESO is based on the FEM. The classical ESO [15] approach is based on the idea of removing inefficient material from an initially oversized domain. The removal process is carried out by deleting regions occupied by elements with low stresses. By repeating this process and removing small amounts of material at each stage, the topology for the structure gradually evolves to a more efficient structure. Following this basic approach, there have since been a number of modifications and refinements such as BESO [16] where not only are elements removed but are also added in high stressed areas. Recent applications of ESO [17] cover a wide range of physical situations by considering element modification sensitivity terms. These sensitivities are used to drive the removal and addition process in order to achieve a minimum (or maximum) of the objective function. Nevertheless, ESO has some drawbacks and weaknesses, notably those related to the rejection criteria, as has been reported by Zhou and Rozvany [18], those due to mesh dependency and those due to the tendency to converge to meshes exhibiting an undesirable single corner contact [19]. Moreover, the final solution, due to the nature of the mesh, can result in jagged edges. It is clear that this is an undesirable situation, since the accuracy of the FE results is in doubt. In addition, post-processing, such as spline construction, must be carried out to smooth the boundary for manufacturing reasons [20]. Similar ideas are presented by Maute and Ramm [21] and Hammer and Olhoff [22] in the material topology optimisation approach, in which the geometry is smoothed using splines based on density distributions. Another disadvantage, common in fixed grid FE-based structural optimisation methods, is the presence of checkerboard patterns [23]. These refer to the phenomenon of alternating presence of void and solid elements ordered in a checkerboard fashion. Such patterns result in structural analysis inaccuracies which complicate the
1903
interpretation of optimal material distribution and geometry extraction. Checkerboard formation is related to the use of low-order elements in the finite element approximations. The use of higher order elements can reduce this effect, but at the cost of increasing the computational time required. Alternatively, to improve the estimation quality of elemental sensitivity in low-order elements, a smoothing algorithm [23] can be implemented. On the other hand, in the ESO-based intelligent cavity creation (ICC) algorithm [24] checkerboard patterns (with numerous cavities) can be eliminated through controlling the number and scale of structural cavities in the final topology. Alternatively, a perimeter control technique [25] can be incorporated into BESO to overcome checkerboard patterns. This technique has also been shown to reduce the dependency of the converged optimum on the initial finite element discretisation. Finally, an important task pointed out by Querin et al. [16] is related to the solution time, that is, the necessity to make the ESO process faster so that the designer is able to obtain the optimum design within a few seconds of describing the environment. Research in this topic has been carried out combining ESO and the fixed grid (FG) method [24], instead of classical FE, to improve the performance. In order to overcome these drawbacks relating largely to domain-based discretisation, the boundary element method (BEM) [26] may be implemented as the tool to carry out structural analysis. In addition, in the current work the geometry is described using the nonuniform rational B-spline curves (NURBS) [27]. The optimisation process is fully integrated within in-house boundary element software [28] allowing a straightforward and rapid communication between the analysis software and optimisation code since they share the same database. The advantages of the proposed algorithm over FEM based ESO procedures are, therefore, the easy remeshing and smoothness of geometry that provides for a smoothness in the stress solution. The use of boundary elements enables further computational efficiency gains beyond a smoothed FE approach (e.g. [21,22]) since many computations may be saved by reusing influence coefficients relating to the boundaries that remain unchanged in successive iterations. The use of the BEM, however, restricts the class of problems to those to which are suitable for BEM treatment. The present work is limited to linear elastostatics, but this might be extended through, for example, dual reciprocity methods, to non-linear problems.
2. Algorithm The developed algorithm considers the BEM to carry out the structural analysis of the component under study. The optimisation approach is stress-based selecting low
1904
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
and high-stressed regions as potential areas to be modified. The geometry of the component is described using NURBS. These parametric curves allow free-form representation with total geometry control over the curve. This feature is very interesting since, a priori, the final design is not known. NURBS are defined by control points which become the design variables of the problem. Another important property of these curves is the locality, that is, changes in one area do not propagate throughout the whole B-spline. The complete optimisation algorithm [29] is presented in detail in the following sections of this paper. First, the main steps of the basic algorithm are summarised as follows: Step 1: The geometry of the structure is defined (quite possibly as a simple rectangle bounding a given design space). Loads and constraints are applied. Step 2: A boundary element analysis (BEA) is carried out. Step 3: The material removal process is performed by selecting the least stressed nodes within the boundary mesh and effectively moving the control points nearest to those nodes. Following a similar process, material addition is carried out if a node is found with a stress higher than the yield stress or a certain maximum stress criterion. The control points nearest to the high-stressed nodes would be moved. However, movement takes place in the opposite direction to that in the removal process. Step 4: To complete the topology optimisation, the inside regions are also considered and internal holes are created in the inner low-stressed areas. As a natural evolution of the process, holes can merge between each other or into the outer boundary. Step 5: Finally, such a procedure is repeated, from Step 2, until the evolution of the objective function shows no improvements. Notice that since, during the optimisation process, the control points are moved and in consequence the shape changes, the boundary element mesh is updated every iteration. The algorithms that are presented in this work have been developed heuristically and, although use is made of parameters that have been used successfully by authors in finite element implementations of ESO, the values of these parameters have been varied according to the results of numerous numerical experiments. These experiments are not described in detail here. Further work is in progress to determine the sensitivity of the process to changes in the parametric values.
2.1. Geometry definition According to the applied constraints, loads or any other design requirements, the boundary domain can be divided into three different types of curves; i.e. design domain, non-design domain and symmetry lines. Lines that can change freely along the process are identified as design domain, whereas those lines that cannot change due to constraints or design restrictions are identified as non-design domain. Symmetry lines can be regarded as an intermediate step between design and non-design domain. They can be modified by changing their length but only if the adjoining line is a design domain line; furthermore, their changes are as a result of changes in the design domain line. 2.2. Boundary element model In the boundary element method (BEM) the boundary integral equations are approximated by a set of discretised integral equations. As a result, the boundary surface is divided into elements, thus the response is given at the nodal points associated with the elements. The main advantage of this method relates to the mesh, since only the surface of the structure needs to be discretised, though their accuracy of boundary stress solution is of course advantageous to this work. The results inside the structure are calculated at an arbitrary number of internal points. In the software used [28], these points are randomly distributed throughout the interior domain. 2.3. Removal and addition of material 2.3.1. Identifying inefficient areas Following the elastostatic boundary element analysis, for any iterative step, material is removed from areas of low stress and added in areas of high stress. Material can be removed from the structure if any node p satisfies ð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). These ratios are conceptually similar to the ones used in FE-ESO [10]. In the customary ESO fashion, if a steady state is reached in which no nodes, or only a few nodes, can satisfy Eq. (1) then RR is increased by the evolutionary rate for removal, ERR, as follows:
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
RRi ¼ RRi1 þ ERR
ð3Þ
where i is the current iteration. Similarly, if only a few nodes can satisfy Eq. (2) then AR is decreased by the evolutionary rate for addition, ERA, ARi ¼ ARi1 ERA
ð4Þ
Typical values for suitable ratios are RR0 = 0.01, ERR = 0.01, AR0 = 0.99, ERA = 0. These values are determined from numerical experience. It should be noted that the von Mises stress is the criterion generally used in this work. Nevertheless, other stress criteria can be considered. For example, the principal stresses r1 and r2 take into account the differences between tension and compression, which may be important in cases having two materials; for example, steel and concrete. 2.3.2. Multi-load case problems In most engineering disciplines multiple load case problems are frequent. For such situations, the final solution depends on the scheme adopted to compute the effect of each load case. The criterion explored in this work is the logical AND/OR algorithm [30]. Here, a node is selected for the removal process only if is lowstressed in all load cases, and it is selected for the addition process if is high-stressed in any of the load cases. Eqs. (1) and (2) refer to a single load case problem. In the case of multiple load cases these equations are slightly changed, thus, for each node p and each load case j, a ratio between stresses is calculated as follows: rp;j Rp;j ¼ ð5Þ rmax;j Computing the above equation for all the load cases would give a set of ratios associated to each node. For example the node p has associated N ratios, where N is the number of load cases. Rp ¼ fRp;0 ; Rp;1 ; . . . ; Rp;N g
ð6Þ
Thus, applying the logical AND/OR, material can be removed from the structure if all the ratios {Rp} are lowstressed ðRp;0 6 RRÞ AND ðRp;1 6 RRÞ . . . AND ð7Þ
ðRp;N 6 RRÞ In the same way, material can be added if ðRp;0 P ARÞ OR
ðRp;1 P ARÞ . . . OR
ðRp;N P ARÞ ð8Þ
where RR and AR are the same optimisation parameters as the ones related to Eqs. (1) and (2).
1905
control points is associated to each specific area to be moved. Every set consists of the three nearest control points to this area. The justification for this choice of three points comes from numerical tests with different numbers of control points. For sets of less than three control points it has been found that the distortion of the curve is very localised causing extreme curvature, while for sets of more than three control points the distortion is propagated to a wider region of the curve and therefore, losing control of the areas to be removed/ added. Each set of control points is moved a distance related to the following parameters: • Length of the least/most stressed element, Le. • Distances of the three control points from the least/ most stressed boundary element. These distances are denoted a, b, c, respectively. • A factor related to the stress situation within the structure at the current iteration, which is called the removal factor (RF) if removing material, and the addition factor (AF) if adding it. This factor allows larger geometric changes in early iterations and smaller changes when fine-tuning is required in the later iterations. The factors RF and AF allow us to gain in computational efficiency without losing the accuracy of the final solution. For example, in Fig. 1(a) the control point (P2) situated at a distance a from the node of lowest stress is moved a distance calculated as follows: 1 0 B Movement of P 2 ¼ Le @
1 C1 A RF 1 1 1 a þ þ a b c
ð9Þ
2.3.4. Direction of movement The direction of movement for each set is perpendicular to its nearest least/most stressed boundary element, i.e. follows the normal vector, n, to the element as shown in Fig. 1(a). The only difference between removal and addition is the inwards movement when removing material and outwards movement when adding it. In the case of the control point at the intersection of two B-splines, as depicted in Fig. 1(b), the direction of movement follows an averaged normal vector. This vector, n, is calculated as the resultant of the two normal vectors, nA and nB , to the B-spline curves involved. 2.4. Creation of holes
2.3.3. Distance to move In any iteration, material is either removed or added to the structure by changing its boundary definition. Since control points define the NURBS curves, a set of
Topology changes are carried out by creating internal cavities or holes in low-stressed areas. These holes are shaped according to the internal von Mises (or any other
1906
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
0
e
1
P2
n
Least stressed node
a
Least stressed element
P4
P Pl
nA
b
P0
P0
B-spline B
P3
c
nB
n = n A + nB
B-spline A
P1 (a)
(b)
Fig. 1. Direction of movement when removing material: (a) one B-spline; (b) two B-splines.
selected criterion) stress contours. The free-shapes are created using NURBS curves to define their geo-metry. In a BEA the results at points contained within the boundary of the object are calculated at internal points. Recalling Section 2.2, in the software used in the current study, these internal points are scattered at random (though concentrated near stress concentrations and their distribution smoothed using a Laplacian algorithm). It is the results at these points, and not at the boundary nodes, that form the starting point of the algorithm for the creation of holes in the domain. Nevertheless, the basis to identify the inner low-stressed areas is similar to the one used for the outside boundary (Eq. (1)). The process of creation of holes is summarised in the following steps:
Step 2:
Step 3:
Step 4:
Step 5:
Step 1: Check for any internal point that satisfies Eq. (10). If so, add this internal point to a stack denoted stackIP ð10Þ
rIP 6 RRrmax
where rIP is the internal point von Mises stress or any other selected criterion, rmax is the maximum von Mises stress or any other se-
Step 6:
IP8
IP14
IP14 IP7
IP8
IP1 IP2
P 10
P20 IP0
IP4
IP9
IP6 IP3
IP13
lected criterion and RR is the removal ratio used in Eq. (1). The stackIP containing the selected internal points is sorted in ascending von Mises stress order. The first point of the stackIP is identified as the minimum and removed from the stackIP to another stack called temporarystackIP. This minimum is taken as the initial central point to be enclosed by the new contour. Internal points surrounding this minimum and having a von Mises stress lower than a threshold (generally related to Eq. (10)) are also moved from the stackIP to the temporary-stackIP. A convex polygon is created from the internal points in the temporarystackIP, as displayed in Fig. 2(a). This polygon is split up into two open polygons and the internal points are taken as the control points. Therefore, these polygons together with the control points define the two NURBS curves that create the new cavity (Fig. 2(b)). If the remaining stackIP is not empty, the process is repeated from point 2 until there are no more internal points in the stackIP.
IP13
NURBS 0
P01 P30
IP10 IP11
(a)
Least stressed internal point IP0
P00
NUR NURBS1
IP5 IP12
P41
P 11
P31
P21 IP11
(b)
Fig. 2. (a) Creation of holes from internal points; (b) new hole consists of two NURBS.
IP9
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
Holes, which are described by two NURBS curves per hole, have the same form of boundary representation as the outside boundary. Consequently, they have similar (but not identical) properties and behaviour. Finally, note that if the von Mises stress in most of the internal points surrounding the minimum is higher than the threshold, then the temporarystackIP consists of only a few internal points and in cases such as this no holes are inserted. 2.4.1. Distance to move control points defining holes To calculate the distance to move the control points in internal holes, a similar procedure to the outer boundary is followed. However, the set of control points associated to each specific area to be moved consists of only one control point, the closest point to this area. The reason for moving exclusively the closest control point is that it is common for stresses around cavities to oscillate between high and low stress regions. A good example is the classical case of a circular hole in a plate under uniaxial tension, in which the tangential stress around the hole
1907
varies in a well known cosine form, containing both tensile and compressive stresses and passing through zero in four locations around the circumference. The use of a three control point set would cause geometric changes to be too widespread and would greatly slow down, or probably entirely prevent, convergence to an optimum. Thus, the control point is moved a distance related to the following parameters: • length of the least/most stressed element, Le; • distance, a, of the control point from the low/high stressed node; • the removal factor (RF) if removing material, or the addition factor (AF) if adding it; • thickness of the component, t. For example, in Fig. 3(a) the closest control point to the node of lowest stress is moved a quantity calculated as follows: 1 ð11Þ Movement of P 2 ¼ Le tRF a
Fig. 3. Direction of movement: (a) isolated hole; (b) hole with a close boundary.
P4 A
P3A
P7C
P6C
Hole A and Hole B merge
P4B P5B
P2 A
P5A P3B
P6A
P2B P1B
P8C Hole C
P0A P8 A
dist |P3B - P6A| < min low stress area
P4C
P9C
P7 A
Hole B P0B
Hole A P1A
P3C P0C
P1C
Fig. 4. Merging two holes (A, B) into a new one (C).
P2C
1908
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
Adjacent holes merge only if they satisfy these two conditions. As shown in Fig. 4, merging the holes implies the creation of a new one. Consequently, this follows a similar procedure to the hole creation process detailed in Section 2.4. In this case, stackIP consists of all the control points related to the holes to be merged.
2.4.2. Direction of movement in holes Generally, the direction of movement for a control point in a hole is perpendicular to its nearest least/most stressed boundary element (Fig. 3(a)). Nevertheless, in the case of a close outer boundary (Fig. 3(b)) the direction of movement is determined by this borderline, i.e. parallel to the boundary that is closest to the selected node.
2.4.3.2. Merging hole to outside boundary. The criterion for merging a hole with the outside boundary is similar to the one for merging two holes (Section 2.4.3.1), that is, minimum distance and under stressed area in between. Fig. 5 depicts how a hole is merged to its closest outer boundary. The process implies a change of the boundary definition and also the deletion of the hole. The extreme
2.4.3. Merging holes It appears that the algorithm works towards an optimum by producing multiple holes which later merge to form larger cavities. Moreover, it is found that, without a hole merging algorithm, the holes approach each other to produce multiple thin filaments of material which are not efficient as a use of material and present unnecessary difficulties in manufacture. 2.4.3.1. Merging two holes. This situation occurs when the cavities evolve approaching too close to each other. The criterion for merging holes is based on the following two premises:
P0B
P 0B
Hole P1 B
P 4H P2B
• Minimum distance. If the distance between two holes is less than a minimum, dmin, those holes are likely to be merged. The distances are calculated between control points from each hole. The largest element length, from the elements in the immediate vicinity, is taken as dmin. • Under stressed area in between. Internal points situated in the region between the holes satisfy Eq. (10).
P1 C(u)
P1H
P3B
P0 H
P4B
New Boundary
P6B
low stress area Fig. 5. Merging hole (hole) into boundary (boundary).
P1
Inserted control point P’
C’(u)
P2
P0
C(u) |P1-P2|i > k|P1-P2|0
(a)
(b) Fig. 6. Insertion of control point P 0 .
P2
P2
P4
C’(u)
C(u)
P3
P4
P1 C(u)
P’ P0
P 4B
dist |P3B – P5H| < min
P2
P1
P3B
Hole and Boundary merge P5B
P5 B
d0 = |P1-P2|0
P0
P2B
P2H
P5 H
Boundary
P1 B
P 3H
P’ P0
P3
|P’-P3| < k’ |P2-P’|
Removed control point (a)
(b) 0
Fig. 7. Removal of control point P .
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
The number of control points defining each NURBS curve is not fixed. On the contrary, control points are automatically either inserted or deleted as required throughout the process. This capability has been implemented in order to control geometry changes and keep smoothness. Another important feature related to the geometry control process is the mesh subdivision. This subdivision is performed in areas in which the NURBS curve distortion has caused excessive curvature of an element. 2.5.1. Insertion of control points The insertion of a new control point between two existing control points is done when these control points are moved away from each other more than a constant factor, k, times the initial distance between them jP j P j1 ji > kd 0
ð12Þ
where i is the iteration number, d0 = jPj Pj1j0 and generally k P 1.5. The control point insertion algorithm is implemented to avoid areas whose defining control points are placed too far away, causing loss of localised control over the B-spline curve. This effect is illustrated in Fig. 6. Fig. 6(a) depicts the original curve C(u) with d0 = P1 P2. As the optimisation evolves, the original curve C(u) changes according to the movement of its control points. In this example, P1 is moving away from the curve compared to the other control points. Fig. 6(b) shows the curve C(u) at iteration i for which Eq. (12) is satisfied. Thus, C(u) is updated to C 0 (u) inserting a control point P 0 to the control point set. The new control point is placed at the mid-point of the line joining both control points P1 and P2.
where k 0 < 1, jPj Pj1j is the length of the actual line checked, jPj1 Pj2j and jPj+1 Pjj are the lengths of the previous and next control polygon lines, respectively. Thus, if in Eq. (13), the inequality on the left is satisfied, then Pj1 would be removed. On the other hand, if it is the inequality on the right that is satisfied then Pj is the control point to be removed. For example, Fig. 7(a) displays an original curve C(u) and the distribution of control points defining the curve. Fig. 7(b) shows the curve C(u) some iterations after, when P 0 and P2 have been moved. At this stage P 0 has approached too close to P3 and the distance between the two control points satisfies Eq. (13) since jP3 P 0 j < k 0 jP2 P 0 j. Finally, C 0 (u) is created by removing control point P 0 and so modifying the control polygon definition.
σ y = 10 N/mm 2 120 mm
2.5. Geometry control
2.5.2. Removal of control points The removal of a control point is performed when a set of control points are moved to locations too close together. This gives rise to regions of excessive and unpredictable curvature of the boundary. The criterion to remove the point from this set is based on comparison of the different lengths of the lines of the control polygon P j P j1 < k 0 P j1 P j2 OR P j P j1 < k 0 P jþ1 P j ð13Þ
σ x = 20N/mm 2 B
30 mm
(maximum and minimum x y coordinates) control points defining the hole are identified and added to the set of control points describing the boundary. In addition, control points from the boundary which are between these extreme points are removed.
1909
y x
A
30 mm
120 mm
Fig. 8. Problem definition for a quarter of the plate.
Fig. 9. Initial and final (a) element mesh; (b) control point distribution.
1910
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
2.5.3. Mesh subdivision An essential element of a fast optimisation algorithm is a reliable automatic mesher to discretise each model. Such a mesher is available in the adopted software, but it has been found important to enhance its conventional automesh algorithm to handle the unusual level of distortion of B-splines in the later iterations. In areas of high curvature, an individual element may need to be subdivided in order to maintain suitable solution accuracy. An element is subdivided if its overall length exceeds certain factor of the distance between the ends of this element as stated in Eq. (14) Le > mpj pj1 ð14Þ
criterion. This stopping criterion can be quantified in the following normalised form: iþ1 f f i 6 104 e¼ ð15Þ fi where fi is the value of the monitored quantity at iteration i and fi+1 is the value of the monitored quantity one iteration after. 2.6.1. Stress levelling The failure of structures under service frequently takes place in areas of locally high stresses. Therefore, it is crucial to avoid stress peaks in order not only to
where Le is the length of the element, e, under study. pj and pj1 are the first and last node of the element and jpj pj1j is the distance between both nodes and generally, m P 1.05.
B Numerical optimum b/a = 0.5521 b
2.6. Stopping criteria Different stopping criteria can be adopted depending on the target to achieve. In this paper two different criteria have been considered, one being based on specific stiffness (which has an equivalence with strain energy) and one on achieving a uniformity of stress over the boundary. Monitoring of the required criterion towards, and ultimately beyond, an optimum provides a stopping
Analytical optimum b/a = 0.5 y A
x a
Fig. 10. Comparison between the theoretical and analytical solution.
Fig. 11. (a) Initial and (b) final von Mises stress contour plot.
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
prevent localised yielding but also to increase the fatigue life of the structure and avoid crack initiation. Stress levelling aims to make uniform the stress along a detailed part of the boundary (C). Therefore, the stress levelling objective deals with minimising the following function: Z fC ¼ ðrp rref Þ2 dC ð16Þ C
structure. However, simply maximising stiffness is not very useful in most engineering contexts since this will result in an optimum in which the design space is completely filled with material. Thus, dividing the stiffness (K) by volume (V) a function termed specific stiffness [17] can be defined as fK ¼
where rp is the node stress and rref is a reference stress. 2.6.2. Strain energy criterion A common objective in engineering is to obtain a structure stiff enough to produce a good performance under certain loads and constraints. Such an optimisation problem consists of maximising the stiffness in the
K V
Von Mises Stress (N/mm2)
Initial design Optimum numerical Optimum analytical
B
y
x
A
10
20
30
40 50 60 Angle (degrees)
ð18Þ
The simplest method of computing strain energy is to use the boundary integral Z 1 U¼ Tu dC ð19Þ C 2
θ
0
ð17Þ
Therefore, the overall specific stiffness is maximised which implies that the stiffness is maximised while the volume is reduced. An equivalent expression can be written in terms of strain energy (U). In this case the objective function adopted is related to the minimisation of the specific strain energy fU as follows: fU ¼ UV
90 80 70 60 50 40 30 20 10 0
1911
70
80
90
Fig. 12. Distribution of von Mises stress along the hole. Fig. 14. Initial and final boundary element mesh.
Fig. 13. Problem definition for a short cantilever beam.
Fig. 15. Initial and final control point distribution.
Fig. 16. Optimum designs: (a) Rozvany et al. [31]; (b) Chu et al. [32]; (c) boundary-ESO.
1912
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
Fig. 17. Von Mises stress contour plot.
where T are the tractions over the boundary and u the displacements over the part of the boundary C where the tractions are applied.
1.0E+04
In this section, several examples are presented in order to prove the developed algorithm. The following isotropic material properties are assumed: YoungÕs modulus E = 210,000 N/mm2, PoissonÕs ratio = 0.3 and an arbitrary thickness t = 1 mm. The yield stress for this material is 280 N/mm2. Plane stress conditions are also assumed.
f = U*V(Nmm4)
9.0E+03
3. Examples
8.0E+03 7.0E+03 6.0E+03 5.0E+03
0
5
10
15 20 25 Iteration number
30
35
Fig. 18. Evolution of the strain energy function.
40
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
3.1. Hole in a biaxial stress field This classical example addresses the shape optimisation of a hole in a large plate under a biaxial stress field. Due to the symmetry, only a quarter of the plate is modelled. Fig. 8 shows the initial shape and boundary conditions. The plate is subjected to stress rx = 20 N/mm2 and ry = 10 N/mm2. The hole is modelled by a quadratic NURBS curve (line B–A). The hole is monitored for stress levelling (Eq. (16)) and in this case, the maximum von Mises stress rmax is taken as rref. The theoretical solution [3] is an ellipse with an aspect ratio b/a equivalent to the loading ratio ry/rx, i.e. 1/2, with the stresses uniformly distributed around the hole. This represents the optimum shape condition, in which all the material at the edge of the hole is at the same maximum permissible stress. The problem was also studied by [4] using the BEM and design sensitivity analysis. The optimisation parameters are set to RR0 = 0.01, ERR = 0.01, AR = 0.99, ERA = 0, RF = 0.05 and AF = 0.025. RF and AF are taken to be constant and of small value in order to keep smoothness and thereby do not give rise to artificial stress peaks over the boundary. Initially, the element model contains 44 quadratic elements, a number that increases to only 45 in the optimised design. Fig. 9(a) shows the initial and final boundary mesh in the area of the hole. There is one quadratic NURBS curve defining the central hole. This curve consists initially of six control points. At the end of the process this number has increased to 9. See Fig. 9(b) for the initial and final distribution of the control points. The final design for the hole is reached after 60 iterations in a total CPU time of 14.162 s on a Pentium 4 (2 GHz). Fig. 10 compares the numerical optimum shape with the theory. In order to compare both results, the analytical optimum is calculated taking the major axis a = 47.8 mm. Since a = 2b, for this problem the minor axis is b = 23.9 mm. Thus, the ellipse (analytical optimum) is calculated following the equation for an ellipse being centred at (0, 0): x2 y 2 x2 y2 þ 2¼1! þ ¼1 2 2 a b ð47.8Þ ð23.9Þ2
1913
compared to the analytical 1.5 [3] and other numerical results 1.59 [4]. A reason for this can be that the theoretical results are calculated assuming an infinite plane whereas in our case the geometry is described by finite segments and the effects of the loading and displacement constraints cannot be avoided. It may also be related to the choice of objective function where rref is the maximum von Mises stress in preference to an average value. Fig. 12 shows the distribution of the von Mises stress along the hole for the initial and final designs. In the arbitrary initial design there is a stress peak of 80 N/mm2 that occurs in a localised area at h = 0. Comparing this result to the (numerical) optimum design; in the latter, the stress concentration has been effectively removed and the stress levels are clearly uniform around the hole. The optimum analytical solution [3] rh = (1 + ry/rx)rx, has a constant value rh = 30 N/mm2. Thus, there is a very good general agreement between the analytical and the numerical optimum. 3.2. Short cantilever beam A common benchmark problem is the cantilever beam with a central load and an aspect ratio of 1.6. The definition of the initial problem is shown in Fig. 13. This cantilever of dimensions 160 mm · 100 mm is fixed at the top and bottom corners of the left hand side; whereas it is subjected to a vertical load at the centre of the free end. The load has been chosen arbitrarily to have value 10 N. This value can be regarded to be very small compared to the failure stress; however it produces an appropriate stress distribution to study removal of
ð20Þ
The numerical results show a slight deviation from the theoretical solution. The ratio of the minor axis of the elliptic hole to the major is b/a = 25.29 mm/ 45.81 mm 0.5521 (analytical optimum b/a = 0.5). Comparing the initial (Fig. 11(a)) and final (Fig. 11(b)) stress distribution, as expected, the stresses around the curve become more uniform when approaching the final solution. The ratio between the maximum von Mises stress at the edge of the hole and the larger of the two biaxial components is 1.691 (rmax/rx = 33.828 N/mm2/20 N/mm2 = 1.691), which is high if
Fig. 19. Initial and final design for a bridge under moving load.
1914
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
material as well as the creation and evolution of holes. The optimisation parameters are AF = 1.0, constant AR = 0.98, ERR = 0.01, RR varies from 0.01 to 0.2, and RF changes from 0.8 to 0.5 according to the stress levels within the structure. The boundary domain is initially discretised into 37 quadratic elements (Fig. 14(a)). At the end of the process this number has increased to 165 elements (Fig. 14(b)).
The number and positions of the control points for the initial and final topology design are depicted in the Fig. 15. The initial five non-constrained straight lines are converted into five quadratic NURBS curves. In the final design this number has increased, with the creation of internal holes, to 11 NURBS. The analytical optimum for this case consists of an infinite number of infinitesimally thin members. Fig.
Fig. 20. Von Mises stress contour plots for LC0, LC1 and LC2.
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
16 compares the optimum design obtained with the approach (boundary-ESO) to those approximations using the optimal layout theory by Rozvany et al. [31] and to a solution obtained by Chu et al. [32] using FEESO. Fig. 16(b) shows the ESO solution obtained under stiffness constraints. The optimum design obtained with the current approach (Fig. 16(c)), shows a good agreement with the FE-ESO. Fig. 17 shows, in terms of von Mises stress contours, how the structure evolves towards the final solution. As the optimisation progresses, small cavities are merged into larger cavities. Also cavities that are very close together or to the outer boundary may merge. Specific stiffness is monitored to form a stopping criterion, and the optimum topology design is reached after 38 iterations. The volume reduction is approximately 63% (V/V0 0.37) from the initial design. Although the overall process shows a clear reduction of volume (removal of material) in the structure, it should be noted that material is also added if any stress concentration appears as a result of the movement of the control points. The final design is obtained in a CPU time of 5 min and 27 s (Pentium 4 (2 GHz)). The evolution of the objective function (Eq. (18)) is displayed in Fig. 18. As can be observed, the function decreases smoothly until iteration 16. At this iteration there is an increment in the function due to the insertion of the first hole. A few iterations later, at iteration 21, two new holes are created. The effect of the creation of these latter holes is not as dramatic as the first one because the size of the holes is much smaller. As the process evolves another two holes are created at iteration 27. These last cavities are very small compared to the existing ones and eventually are merged to the larger holes. The final topology is achieved when no further improvements (Eq. (15)) are shown in the objective function which shows a 28% reduction compared to the initial situation. 3.3. Bridge under a moving load This final example tackles the problem of a bridge structure under multiple loading conditions. The initial design for the bridge is shown in Fig. 19(a). The dimensions of the bridge are 160 mm · 50 mm. The bridge is supported at four points in the bottom surface. A point load of 100 N is moving along the top surface and this effect is approximated by five load cases. The logical AND/OR scheme represented by Eqs. (7) and (8) is adopted. Point boundary loads and constraints are difficult to apply in the BEM, and consequently the loads and constraints in this solution are applied over sections of boundary of finite length. The area over which the constraints are applied are evident in the converged optimal geometry.
1915
The top surface is considered to be in the non-design domain as well as the four supports. Moreover, there is a limit for the maximum vertical displacement of 0.02 mm. The rest of the domain is free to change and defined by 5 NURBS curves. The optimum design, shown in Fig. 19(b) that satisfies this problem is found after 42 iterations with RR varying from 0.01 to 0.27. Fig. 20 shows the von Mises stress contours, at the initial and final iteration, for load cases LC0, LC1, LC2. LC3, LC4 are symmetric to LC1, LC0, respectively. As expected, the final solution under multiple load cases is not necessarily a fully stressed design [20]. The ratio of the final volume to the initial volume is 0.54. The maximum vertical displacement (umax) for any load case is below the limit constraint of 0.02 mm. Certainly, these results would have varied if the constraints over the maximum vertical displacement for each load case had been considered differently.
4. Conclusions This paper presents a new algorithm based on the evolutionary structural optimisation (ESO) method. The algorithm exploits the concept that the optimum topology evolves by slow removal and addition of material. The boundary element method is used for the analysis and NURBS curves for describing the changeable lines. The control points defining these curves are selected to be the design variables of the problem. Topology optimisation is accomplished allowing the insertion of internal holes in the inner low-stressed areas of the structure. Holes are also described by NURBS curves and so they have similar behaviour to the outside boundary. To prove the algorithm, several examples have been reproduced showing the effectiveness of the method generating the final topologies as well as rapid solution times. The approach presents the advantage of decreasing the number of design variables compared to FEbased ESO. Also, it has been proved that working directly with the boundary provides a large flexibility in the design. It has been shown how at each iteration the boundary remains smooth without artificial stress concentrations.
Acknowledgments The authors would like to thank Prof. G.P. Steven for inspiring advice and useful discussion. The first author gratefully acknowledges the Durham Research studentship Isabel Fleck awards, Doreen Bretherton studentship and Lindsay scholarship provided during this study.
1916
E. Cervera, J. Trevelyan / Computers and Structures 83 (2005) 1902–1916
References [1] Haftka RT, Grandhi RV. Structural shape optimization, a survey. Comput Methods Appl Mech Eng 1986;57(1): 91–106. [2] Braibant V, Fleury C. Shape optimal design using Bsplines. Comput Meth Appl Mech Eng 1984;44:247–67. [3] Kristensen ES, Madsen NF. On the optimum shape of fillets in plates subject to multiple in-plane loading cases. Int J Numer Methods Eng 1976;10:1007–19. [4] Tafreshi A, Fenner RT. Design optimization using the boundary element method. J Strain Anal 1991;26(4): 231–40. [5] Schnack E, Spo¨rl U. A mechanical dynamic programming algorithm for structure optimization. Int J Numer Methods Eng 1986;23:1985–2004. [6] Mattheck C, Burkhardt S. A new method of structural shape optimisation based on biological growth. Int J Fatigue 1990;12(3):185–90. [7] Bendsøe MP, Kikuchi N. Generating optimal topologies in structural design using a homogenisation method. Comput Methods Appl Mech Eng 1988;71:197–224. [8] Rozvany GIN, Zhou M, Birker T. Generalized shape optimization without homogenization. Struct Optim 1992;4:250–2. [9] Baumgartner A, Harzheim L, Mattheck C. SKO (soft kill option): the biological way to find an optimum structure topology. Int J Fatigue 1992;14(6):387–93. [10] Xie YM, Steven GP. Evolutionary structural optimization (ESO). Springer; 1997. [11] Eschenauer HA, Kobelev VV, Schumacher A. Bubble method for topology and shape optimisation of structures. Struct Optim 1994;8:42–51. [12] Hassani B, Hinton E. A review of homogenization and topology optimization III—topology optimization using optimality criteria. Comput Struct 1998;69:739–56. [13] Cerrolaza M, Annicchiarico W, Martinez M. Optimization of 2D boundary element models using b-splines and genetic algorithms. Eng Anal Boundary Elem 2000;24:427–40. [14] Frangopol DM, Maute K. Life-cycle reliability-based optimization of civil and aerospace structures. Comput Struct 2003;81:397–410. [15] Xie YM, Steven GP. A simple evolutionary procedure for structural optimization. Comput Struct 1993;49(5):885–96. [16] Querin OM, Xie YM, Steven GP. Evolutionary structural optimization (ESO) using a bidirectional algorithm. Eng Computat 1998;15:1031–48. [17] Steven GP, Qing L, Proos K, Xie YM. The role of physical sensitivity in evolutionary topology design optimisation
[18]
[19]
[20]
[21] [22]
[23]
[24]
[25]
[26] [27] [28]
[29]
[30]
[31] [32]
with multi-criteria and multi-physics. In: Proceedings of the Fifth World Congress on Computational Mechanics (WCCM V). Vienna, Austria, 2002. Zhou M, Rozvany GIN. On the validity of ESO type methods in topology optimization. Struct Optim 2001;21: 80–3. Sigmund O, Petersson J. Numerical instabilities in topology optimisation: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Struct Optim 1998;16(1):68–75. Hinton E, Sienz J. Fully stressed topological design of structures using an evolutionary procedure. Eng Computat 1995;12:229–44. Maute K, Ramm E. Adaptive topology optimisation of shell structures. AIAA J 1997;35:1767–73. Hammer VB, Olhoff N. Topology optimization of continuum structures subjected to pressure loading. Struct Multidiscip Optim 2000;19:85–92. Li Q, Steven GP, Xie YM. A simple checkerboard suppression algorithm for evolutionary structural optimization. Struct Multidiscip Optim 2001;22:230–9. Kim H, Querin OM, Steven GP, Xie YM. A method for varying the number of cavities in an optimized topology using evolutionary structural optimization. Struct Multidiscip Optim 2000;19:140–7. Yang XY, Xie YM, Liu JS, Parks GT, Clarkson PJ. Perimeter control in the bidirectional evolutionary optimization method. Struct Multidiscip Optim 2003;24:430– 40. Trevelyan J. Boundary elements for engineers. Computational Mechanics Publications; 1994. Piegl L, Tiller W. The NURBS book. second ed. SpringerVerlag; 1997. Trevelyan J, Wang P. Interactive re-analysis in mechanical design evolution. Part I. Background and implementation. Comput Struct 2001;79(9):929–38. Cervera E. Evolutionary structural optimisation based on boundary representation of B-spline geometry. PhD Thesis. University of Durham, Durham, UK, 2003. Li Q, Steven GP, Xie YM. On equivalence between stress criterion and stiffness criterion in evolutionary structural optimization. Struct Multidiscip Optim 1999;18:67– 73. Rozvany GIN, Bendsøe MP, Kirsch U. Layout optimization of structures. Appl Mech Rev 1995;48(2):41–119. Chu DN, Xie YM, Hira A, Steven GP. On various aspects of evolutionary structural optimization for problems with stiffness constraints. Finite Elem Anal Design 1997;24(4): 197–212.