Meccanica31: 143-161,1996. @ 1996 Kluwer Academic Publishers.

Printed in the Netherlands.

On Optimum Design of Structures and Materials* NIELS OLHOFF Aalborg Universiv,

Institute of Mechanical

Engineering;

DK-9220 Aalborg East, Denmark

(Received: 7 July 1995) A surveyof problems of optimum design of structuresand materialsis presentedwitb the main empbasis on fundamental aspectsand on current methods and capabilitiesfor topology and shapeoptimization. The methods are selectedfrom conditions of versatility and suitability for integration into an engineering design optimization system which realizes the design process as an iterative solution procedure of a multicriterion optimization problem based on the concept of integration of finite element analysis, sensitivity analysis, and optimization by mathematical programming. A picture of current possibilitiesand the presentstatusof the field is given through a number of examples. Abstract.

Sommario. Si passanoin rassegnavari problemi di progetto ottimale di strutturee materiali, ponendo in particolare lute gli aspettidi baseed i piti attuali procedimenti e strumentiperl’ottimizzazione della topologia e della forma. I procedimenti presentati sono scelti tra quelli the possonoessereintegrati in un sistema di ottimizzazione the realizzi il process0della progettazione ingegneristica come una procedura di soluzione per successiveapprossimazioni di un problema di ottimizzazione multicriterio, basato sull’integrazione dell’analisi ad elementi fin&i, analisi di sensibilita e ottimizzazione per mezzo della programmazionenumerica. Si presentano infine alcuni esempi volti a dare una panoramicadelle attuali possibilita dell’ottimizzazione strutturale. Key words: Optimization, Mathematical programming, Structural mechanics,Mechanicsof materials.

1. Introduction

Engineeringactivity has always involved endeavourstowardsoptimization, andthis particularly holds true for the field of engineeringdesign.Earlier, engineeringdesignwas conceived as a kind of ‘art’ that demandedgreat ingenuity and experienceof the designer,and the development of the field was characterizedby gradual ~V&&WZ in terms of the continual improvement of existing types of engineeringdesigns.The designprocessgenerally was a sequential‘trial and error’ processwherethe designer’sskills and experiencewere the most important prerequisitesfor successfuldecisionsfor the ‘trial’ phase. In contrast,today’s strong technologicalcompetition which requiresreduction of design time and using of products with high quality and functionality with resulting high costs, togetherwith the currentemphasison saving of energyandreuseof materials,consideration of environmentalproblems,etc., often involves the creationof new productsfor which prior engineeringexperienceis totally lacking. The developmentof suchproductsmust naturally lend itself towardsscientific methods. Hence,during recentdecades,engineeringdesignhaschangedfrom ‘art’ and ‘evolution’ to scientifically basedmethodsof rational designand optimization. This developmenthasbeen strongly boostedby the advent of reliable generalanalysis methodssuch as finite element analysis, design sensitivity analysis, and methods of mathematicalprogramming, together * General invited lecture presentedat the 12thItalian Congresson Theoretical and Applied MechanicsAIMETA ‘95, Naples, 3-6 October 199.5.

144

Niels Olhoff

with the exponentially increasingspeedand capacityof digital computers.Thus, methodsof rational design and optimization arenow finding widespreaduse in aeronautical,aerospace, mechanical,nuclear,civil, and off-shoreengineering.In materialsscience,the techniquesare now being usedin researchdevotedto tailoring of materialswith specific properties. The developmenthasbeensupportedby vigorousresearchin thefields of designsensitivity analysis and optimum design; see,e.g.,the review papersby Haug [25], Schmit [61], Olhoff and Taylor [41], Haftka and Grandhi [22], Ding [14], Rozvany et uZ. [59], the monograph by Bendsoe[4], and the proceedingsfrom various conferencesand symposiapublished by Haug and Cea [24], Morris [33], Eschenauerand Olhoff [15], Atrek et uZ. [2], Bennet and Botkin [lo], Mota Soares[34], RozvanyandKarihaloo [57], Eschenaueret UL [16], Cinquini and Rovati [13], Bendsoeand Mota Soares[6], Haug [23], Pedersen[46], Rozvany [53], Herskovits [26], Gilmore et uZ. [21], andOlhoff andRozvany [40]. It is the objectiveof the currentpaperto presentbasicconceptsfor problemsof engineering designtreatedas multicriterion optimization problems,andto give a brief surveyof selected enablingmethodsillustratedby examples.Someof thetechniqueshavebeendevelopedwithin the Danish TechnicalResearchCouncil’s Programmeof Researchon ComputerAided Design in fruitful collaborationbetweenresearchersfrom the Institute of Mechanical Engineeringof Aalborg University, and the Mathematical Institute and the Departmentof Solid Mechanics of the Technical University of Denmark.Many colleaguesandfriends from abroadhave also participatedin the work, and their inspiring cooperationis gratefully acknowledged.Most of the examplesin the paperhavebeenobtainedusing the Qptimum mign sstem ODESSY (see[54], [56], [38], [39] and [29], [37]), which is being developedat Aalborg University. 2. Basic Concepts

The label engineeringdesign optimization identifies the type of design problem where the setof structuralparametersis subdividedinto so-called‘preassignedparameters’and ‘design variables’, and the problem consistsin determining optimum values of the designvariables suchthat they muximize or minimize a specific function termed the objective (or criterion, or cost) finction, while satisfying a set of geometricaland/orbehaviouralrequirementswhich are specifiedprior to design,and arecalled construints. The designvariableswill be denotedby ui, i = 1, ..., 1, and are assembledin the vector a. They canbe categorizedas follows. - Geometricul design vuriubles: l Sizing design vuriubles: describecross-sectionalpropertiesof structuralcomponents

l

l

l

suchasdimensions,cross-sectionalareasor momentsof inertiaof bars,beams,columns and arches,or thicknessesof membranes,platesand shells. ConjigurutionuZ design vuriubles: describethe coordinatesof the joints of discrete structuressuchas trussesand frames;or the form of the center-lineor mid-surfaceof continuousstructureslike curvedbeams,archesandshells. Shape design vuriubles: governthe shapeof externalboundariesand surfaces,or of interior interfacesof a structure.Examples are the cross-sectionalshapeof a torsion rod, column or beam; the boundary shapeof a disk, plate, or shell; or the shapeof interfaceswithin a structuralcomponentmadeof different materials. Topologicul design variubles: describethe type of structure,numberof interior holes, etc., for a continuous structure.For a discretestructurelike a truss or frame, these

On Optimum Design of Structures and Materials

145

variablesdescribethe number,spatial sequence,and mutual connectivity of members andjoints. representconstitutive parametersof isotropic materials, or, e.g.,stackingsequenceof lamina, andconcentrationandorientationof fibersin composite materials.

- MateriaZ design variabk

- Support design variables:

describethe support(or boundary)conditions,i.e., thenumber, positions and types of supportfor the structure.

- Lhading

design variabzes: describethe positioning and distribution of external loading which in somecasesmay be also at the choice of the designer.

- Manufacturing

design variables: parameterspertainingto the manufacturingprocess(es), surfacetreatment,etc., which influencethe propertiesand cost of the structure.

The objective function and constraintfunctions that enter the definition of a given optimization problem must be expressedmathematicallysuchthat their valuescan be determined for any combination of admissiblevaluesof the designvariables.In orderto allow for a broad range of problem formulations, with easy interchangeof typical objective and constraint functions from oneproblem to another,it is convenientto apply the common notion ProbZem Variables (PVs) for thesefunctions.Distinction may be here madebetweengZobaZ cost PVs and gZoba1 and1ocaZ behaviouralPVs, of which the following examplesmay be given. l

GZobaZ cost PVs: structural volume or weight, material cost, manufacturing cost, life

cycle cost, etc. l

PVs: compliance, integral stiffness,buckling load, plastic collapse

Global behavioural

load, eigenfrequency,dynamic response,fatigue life, etc. l

Local behavioural

PVs: stresscomponents,vonMises or Trescareferencestresses,strains,

displacements,etc. 3. Mathematical FormulaGon for Muhicriterion Optimization In practical probl%msof optimum engineeringdesign, it is necessaryto take into account several performance and failure criteria in the problem formulation, i.e., a multicriterion approachmust be adopted.Here, the multicriterion problem will be cast in scalarform by stating it as minimization of the maximum of a weighted set of global criterion functions fj, j = l,... J and local criterion functions f;, k = 1, . . . A’. Thus the designobjective is consideredin the initial form

where tij and ti; are given weighting factors for the separatecriterion functions, and where it is assumedthat eachof the criteria havebeensuitably preconditionedfor minimization. In Equation (1), x designatesthe coordinatesof any point within the structuraldomain 0 over

146

Niels Olho#

which the maximum valuesof the weightedlocal criterion functions fl, k = 1,. . . K (e.g., von Mises stressessubjectto different loading cases)are sought. Notice that Equation (1) only presentsthe objective of the multicriterion optimization problem. The full problem consists of constraintsas well, along with finite element state equationsand other equationsthat are necessaryfor calculation of all the PVs that enter in the objective function andthe constraintsof theproblem. Theseconstraintsandfinite element equationsof the full problem arelisted below, asEquations(5)-(g). The objective (1) of the multicriterion optimization problem is equivalentto

To circumvent the inherentdifficulty thatthis min-max objectiveis generallynot diflerentiuble with respectto the designvariablesui, i = 1, . . . , I, we developa bound formulation of the problem; seeBendsoeet al. [9], andOlhoff [35]. By this techniquewe introducean additional variable, @,which is termed the bound purameter, andwrite Equation (2) in the form of the di#erentiubZe probZem definedby Equations(3) and (4) below. Thus, we simply transform the statement(2) into the constraints(4) and minimize the common upperbound p on theseconstraintsaccordingto Equation (3). Notice that the upper bound parameterp is an additional variable which replacesthe original non-differentiable objective (1) or (2) and is to be minimized over a constraint set in an enlargedspace.The original points of non-differentiability correspondto ‘comers’ in the constraintset (4) of the enlargedspace,and arisefrom intersectionsof differentiableconstraints. In addition to Equations(3) and (4) that representthe multicriterion design objective (1) and (2), we now write as Equations(5)-(8) the original constraints,etc., of the full problem as mentioned earlier, and arrive at the following bound formulation of the muZticriterion optimization

problem:

(3) Subject to: Variable bound constraints:

(4) Original

constraints:

(global) p=l,...,P gz 5 Gz , s = 1,. . . , S, V x G 0 (local).

&&$I,

(5)

Side constraints for design variables:

(7)

On Optimum Design of Structuresand Materials

147

Finite elementstate equations: (KD = F)t , t = 1, . . . , T, etc.

@I

In Equation (5), &‘p and ez denoteprescribedupper values of global and local constraint functions gr, p = 1,. . . , P, and gz, s = 1,. . . , S, respectively,and ai and pi are given lower andupper sideconstraintvaluesfor the designvariablesui, i = 1, , . . ,I. Equation (8), which may alsoinclude eigenvalueproblems,etc.,representsall the setsof finite elementstate equationsnecessaryfor computationof the different displacementvectorsD requiredfor the determinationof the various global andlocal PVs in Equation(7). Notice that while the local behaviouralPVs dependon the designvariablesa;, the pertinentdisplacementvector D, and the coordinatesz of the structure,then global behaviouralPVs are independentof Z, and global cost PVs are moreoverindependentof displacementsD. Actually, in practice,there is no needto make a distinction betweenthe global and local PVs. Thus, the local PVs arejust evaluatedat a number of nodal points of the finite elementdiscretizedstructureaccordingto some suitableactive setstrategy,andtherebyresult in a numberof inequality constraintsthat have the samemathematicalform as thosefor the global PVs in Equations(4) and(5). While Equations (7) and (8) representthe necessarytools for finite element analysis, Equations (3)-(6) constitute a standardproblem of mathematicalprogramming which may be solved iteratively using either a SIMPLEX algorithm, the CONLIN optimizer [19] or the Method of Moving Asymptotes [67], which require calculation of first order design sensitivities of all the PVs. For this purpose,a sensitivity analysis is performed for each design variableat eachstepof iteration. 4. Design Sensitivity

Analysis

In finite elementbasedmethodsof design sensitivity analysis,the sensitivitiesof the various behaviouralPVs enteringin the variableboundconstraints(4) andoriginal constraints(5) are determinedfrom designsensitivitiesof the samepertinentdisplacementfields D (obtainedin Equation (8)) from which the PVs arethemselvescomputedby useof Equation(7). Therefore, when the displacementsensitivities are known, the sensitivities of the behaviouralPVs are generally easily computed(see,e.g., [37]). Assume now that after factorization of the global stiffnessmatrix K, we have solved one of the stateequations(8) for the displacementvector D that is relevantfor computation of a particular PV. The direct approachto obtain designsensitivities of the displacementvector is basedon impZicit diflerentiation of thefinite eZementstate equations.If Equation (8) is differentiated with respectto a design variable u; and the terms arerearranged,the following expressionis obtainedfor computationof the displacementsensitivities i3D/& ,

Equation (9) is of the sameform asEquation(8),so the factorized stiffnessmatrix K can be re-used,and only the new right hand side, which is termed the pseudo Zoadvector, needs to be calculated before the sensitivities ~D/IYu~ can be found by forward and backward substitution.The derivativesaF/%i of the force vector in Equation (9) are generallyeasily calculated(notethat they vanish for designindependentloads),andthenthe determinationof

148

Niels Olhoff

in Equation (9) only requirescalculation of the designsensitivities aK/&zi of the global stiffnes matrix. The derivativesarenormally calculatedat the elementlevel, i.e.

~D/&I;

where k is the element stiffness matrix and 7~~is the number of finite elements. If the design sensitivities i3k/& are determinedanalytically before their numerical evaluation, the approachis called unulyticd design sensitivizyundysis, and if they are determinedby numerical differentiation, the method is called s~~Avuz~~Gcu~ designsensitivi~ undysis, cf. [69], [17]-[18] and [12]. The method of undyticd design sensitivity analysis is generally quite cumbersometo implement in a general optimum design system; a large amount of analytical work and programming is requiredin order to developanalytic expressionsfor derivatives of various stiffness matrices with respectto all possible types of design variables.It is much more attractiveto usethe methodof semi-unulyticuldesignsensitivityundysis in this context, asit is easyto implement for different kinds of designvariablesand finite elementtypes,because simple and computationally inexpensivefirst order finite differencesareused.Therefore,the method of semi-analytical sensitivity analysisis popular and, in general,very efficient and reliable. It has beenrecently shown by way of examples[3] that the method may yield erroneous shapedesign sensitivities for plates, shells, and long slenderstructures.However, a simple, easily implemented and computationally inexpensivemethod of ‘exact’ semi-analytical sensitivity analysisthat completely eliminates the deficiency is now available (see[38] and [3W 5. Shape Optimization

of Discrete Structures

While optimization of discrete2-D or 3-D trussor frame structureswith cross-sectionalareas of bar or beam membersas designvariablesis termed sizing optimizution, and optimization with positions of joints asdesignvariablesis called conjgurution optimizution, the combined problem of optimizing both sizing andconfigurationis labeledshupeoptimizution of the discretestructure.Minimization of the structuralweight is most oftenthe designobjectivein such problems which usually encountermultiple loading conditionswith constraintsthat include stresses,displacements,and elastic as well as plastic buckling of members in compression (see,e.g., [43-45], [36] and [62]. Problems of this kind are seento be well within the scopeof the generalformulation in Equations(3)-(8), anda simple exampleis presentedlater in Figure 8(B). For theseproblems, schemesof analytical designsensitivity analysismay be developedwithout excessivework. 6. Shape Optimization

of Continuum

Structures

In shapeoptimization of continuum structures,the goal is to determine the shapeof the structural domain, i.e., the problem is defined on a domain that is unknown u priori. Like optimization problems in general,shapeoptimization is a highly non-linearproblem where it is necessaryto employ iterative numerical solution schemesand to determinethe optimum designthrougha sequenceof redesignandreanalysis.This implies that the structuralgeometry must be repeatedlyconvertedinto a finite elementmodel with the properloadsand boundary

On Optimum Design of Structures and Materials

149

Figure I, Design model of turbine disk: (A) problem model, (B) variable design model.

conditions, andthat the variable structuraldesignmust be described(parameterized)in terms of a finite number of geometricalvariables. Particularly for problems of shapeoptimization of continuum structures,it is necessary to make a clear distinction betweenthe analysis model as representedby the finite element model, andthe parameterizedgeometricmodel of the variable structurewhich is termed the design modeZ (cf. [17-181, [l 11,[53] and [36]). The designmodel is endowedwith additional significancebecauseit canbe closelyconnectedwith a CAD model asdescribedby Rasmussen [53], [%I, and Olhoff et al. [38]. The design model may consist of so-calleddesign elements as presentedby Braibant and Fleury [ 111.The boundaries(or surfacesin the caseof a three-dimensionalmodel) of the design elementscan be curves of almost any character,i.e., piecewise straight lines, arcs, b-splinesof any given degreeof continuity, Bezier curves,Coonspatches,etc. It is therefore very simple to generaterelatively complicated geometrieswith a small number of design elements.The shapesof the boundariesare controlled by a number of control points, also often termed master nodes. For exemplification, let us considera problempresentedby Lund [29], namely theaxisymmetric model shown in Figure l(A) of a turbine disk which rotates at 2094 rad/s (= 20 000 rev/min). The bladeshavebeenreplacedby a uniformly distributedload at the rim of the disk which representsthe centrifugal forcesfrom the blades. Furthermore,as is also indicated in Figure l(A), the turbine disk is subjectedto different temperatures.Thesetemperaturesderive from the hot exhaustgaswhich drives the disk. Part of the boundary is subjectedto a convection boundarycondition which is also due to the exhaustgas.At theseboundariesthe temperatureof the environmentis specifiedto be 45O’C and the convectioncoefficient is 0.0012W/(mm2K). The design model as shown in Figure l(B) consistsof only two designelements.There are two design boundaries,i.e., boundarieswhose shapesare allowed to change.Each of these shapesare defined by the positions of a number of masternodes,and this createsan evidentconnectionbetweenthe designvariables(themovementsof the masternodes)andthe shapeof the geometry.In this example, the direction of the movementof each masternode is constrainedto follow some predefinedmove directionsspecifiedby the designeras shown in Figure l(B). Such translationaltransformationsare very common in definitions of shape designmodels.Otherpossiblemodifiers may be point scaling,line scaling,rotation,etc. Thus,

150

Niels Olho3Y ollK~

z.~.o|...o z.ttT..... i.m.l.n~

I~pr~f

1.1411.11 i

..ll~.Mi

s,zl=x-~s

VN i , , l l l l S l J l t l t

~lllllll41tl.lU'l.t~lllll I ,I.!

v~ ~,,.l..i,..~,~ylM~tl~l.~llllllllU~l.l..l~lllllll

I*l,i

(A) F i g u r e 2. Finite element discretization of the analysis model, with indication of local magnitudes of ratio between von Mises and temperature dependent yield stress in Equation (11): (A) initial analysis model, 03) final model corresponding to the optimum design.

in this example, the design variables ai, i = 1 , . . . , I, are simply taken to be the magnitudes of the movements of the master nodes along the associated move directions. In addition, the distribution of finite element nodes on the boundaries and the desired finite element type for each design element must be defined. All necessary specifications including loading conditions are assigned to the design model and a preprocessor with automatic mesh generation can automatically convert it into a finite element analysis model. The discretization of this model is shown in Figure 2(A), where the preprocessor has meshed the design elements with a mixture of 6 and 9 node isoparametric 2-D axisymmetric finite elements. All load specifications are automatically converted into consistent nodal loads on the analysis model. The design model has now been converted into an analysis model, and by changing the values of the design variables, the geometry can be changed parametrically into other shapes. We now perform optimization of the turbine disk. The objective will be to minimize the mass moment of inertia of the disk. At the same time, we specify a temperature dependent non-linear yield stress constraint on the von Mises stresses. This constraint is realized by normalizing the von Mises stress trvm,d in any nodal point j of the analysis model with the yield stress given by try = 550 (MPa) if Tj < 623~ and try = 1484.5 - 1.5 Tj (MPa) if Tj > 623~ see Equation (11). It should be noted that at given points of the structure this non-linear stress constraint changes during the optimization process as the temperature field changes with design. Furthermore, we specify the manufacturing constraint that the radius of curvature of the boundary of the disk may nowhere be less than 3 mm. Thus, the definition of the optimization problem is Minimize the mass moment of inertia { ~rvmd[MPa]( 1 550

Subject to stress constraints ~r. . . .

--

if Tj < 623~

i[MPa]

(11)

1484.5_1.5Td ~ 1 if Td > 623~ minimum boundary radius of curvature, 3 mm. The ratio between the von Mises stress and the temperature dependent yield stress is displayed for the initial structure in Figure 2(A) which shows that from the outset the stress constraint is violated by 34%. A shape optimization of the disk is now performed, whereby

On Optimum Design of Structures and Materials

151

the mass moment of inertia is reducedfrom 3.12 . 104to 2.33 . 104kg.mm2. The optimum geometryis shown in Figure 2(B) which illustratesthat the stressconstraintis now satisfied. The final minimum value of the radiusof curvatureof the boundaryis found to be 13.5 mm, i.e., this constrainthasnot becomeactive. 7. Topology Optimization It is characteristicthat the solution to a shapeoptimization problem (cf. Sections 5 and 6) will always have the same topology as that of the initial design. Thus, the topology of a mechanical structureor componentto be optimized cannotbe changedby using methodsof shapeoptimization. As the choice of the besttopology hasconsiderableimpact on the gain to be achievedby designoptimization, the developmentby BendsoeandKikuchi [8] of a method for topology optimization was a remarkablebreak-throughin the field of optimum design, The readeris referredto an exhaustivemonographby Bendsoe[7] for recentdevelopments and publications. Contrary to shapeproblems, a problem of topology optimization is defined on a jked domain of space, and the structureis consideredas a spatial sub-domainwith high density of material. Basically the problem of topology optimization is one of discreteoptimization, but this difficulty is avoided by introducing relationshipsbetweenstiffnesscomponentsand density,basedon the physical modeling of porous,periodic microstructureswhoseorientation and density are describedby continuousvariablesover the admissible design domain, The solution of this problem is basedon a finite elementdiscretizationof the admissibledomain, andthe optimum valuesof the designvariables(densityandorientationof themicrostructures) are determinediteratively. More precisely,for the topology optimization we minimize compliancefor a fixed, given volume of material, andusea densityof materialasthedesignvariable.The densityof material andthe effective material propertiesrelatedto thedensityis controlledvia geometricvariables which govern the material with microstructurethat is constructedin order to relate correctly material density with effective material property. The problem is thus formulated as: Minimize L(w) Subject to u~(w, w) = L(v) Volume 5 V

Vv E H

w

where

Here,B; andTi arethe bodyforcesandsurfacetractions,respectively,andQ denotelinearized strains.R is the set of kinematically admissibledisplacements.The problem is definedon a fixed referencedomain 0 andthe componentsof the tensorof elasticity ,?$jkl dependon the design variables used. For a so-calledsecondrank layering constructedas in Figure 3, we have a relation

152

Niels Olhoff

MACRO-SCALE MICRO-SCALE

1 MICRO-SCALE 2

Figure 3. Construction of a layering of secondrank.

Figure 4. (A) Periodic microstructurewith squareholes rotated the angle @,(B) Squarecell with a squarehole.

where p and 7 denote the densities of the layering and 0 is the rotation angle of the layering. The relation (14) can be computed analytically (see [4], [7]) and for the volume we have Volume

=

/-

p

+ 7 - P7) cm*

Layered materials is just one possible choice of microstructure that can be applied. The important feature is to choose a microstructure that allows density of material to cover the complete range of values from zero (void) to unity (solid), and that this microstructure is periodic so that effective properties can be computed (numerically) through homogenization (theory of cells). This excludes circular holes in square cells, while square or rectangular holes in square cells, see Figure 4, are suitable choices of simple microstructures. For the case of a rectangular hole in a square cell, the volume is also given by Equation (15), with /J and 7 denoting the amount of material used in the directions of axes of the cell and hole; for this microstructure the angle 0 of rotation of the unit cell becomes a design variable. The optimization problem can now be solved either by optimality criteria methods as in [7], or by duality methods [20], where advantage is taken of the fact that the problem has just one constraint. The angle 19of layer or cell rotation may be controlled via results on optimum rotation of orthotropic materials as presented by Pedersen [48-49], see also Sacchi et al [60]. As stated above, the optimum topology is determined from the condition of minimum compliance subject to a bound on the total structural volume. In optimum design, on the other hand, there is a need to handle a much larger variety of formulations of which stress and volume minimizations are the most frequent. However, in spite of the incompatible formulations, compliance optimized topologies tend to perform well also from a stress minimization point of view. This is due to the fact that a relatively high amount of energy is stored in possible areas of stress concentration which then become undesirable also in the compliance minimization.

On Optimum Design of Structures and Materials

153

Distributed

Figure 5. Available design space,loading, and foundation for bearing pedestal.

The initial topology optimization can therefore in many cases lead to substantial improvements of the final result even though the actual aim is to perform an optimization of a different type. The topology optimization thus results in a prediction of the structural type and overall lay-out, and gives a rough description of the shape of outer as well as inner boundaries of the structure. This motivates an integration of topology and design optimization, and in order to gain the full advantage, it is necessary that they be integrated and implemented in a flexible, user-friendly, interactive CAD environment with extensive computer graphics facilities (see, e.g., [531, [361, [381, WI, and VI. It will now be demonstrated via two examples that, depending on the amount of material available, the generated topology will basically define either the rough shape of a continuum structure, possibly with macroscopic interior holes (which shape optimization procedures cannot create), or the skeleton of a truss- or frame-like structure with slender members. Thus, topology optimization can be used as a preprocessor for refined shupe optimization carried out as discussed in Section 5 for discrete structures and in Section 6 for continuum structures. As a first example [53], we aim at minimizing the maximum von Mises stress in the bearing pedestal in Figure 5. The initial geometry completely occupies the available design space. We prescribe that 20% of this space be an upper limit for the structural volume, and start the topology optimization with an evenly distributed density of 20% material (in square microscale cells with rectangular holes) throughout the structure. An exception is made for the rim of the hole of the bearing, where we prescribe 100% density of material. Solving now the topology optimization problem with symmetry taken into account, we obtain the optimum structural topology shown in Figure 6 where black areas represent solid material and white areas represent void. In this example the interfaces between solid material and void are relatively clear and leads us to define a design model for shape optimization of the bearing pedestal as shown in Fig 7(A). Here, out of practical considerations, the initial shape has been slightly modified in comparison with the topology optimized geometry of Figure 6; the additional material provides a basis for attachment to the foundation by, e.g., a bolt joint, and the shape design model is defined such as to ensure a minimum thickness of this region. We now formulate the subsequent shape optimization problem as minimization of the structural volume subject to given upper bounds on the maximum von Mises stress in the structure and on the vertical displacement of the loaded inner surface, The final finite element model visualizing the optimum design is shown in Fig 7(B).

154

Niels Olhojf

Figure 6. Optimized topology of bearing pedestal(symmetry assumed).

Figure 7. (A) design model, (B) final finite element model of shapeoptimized bearing pedestal.

As a secondexample[36], we considerFigure 8(A) which presentsthe result of optimizing the topology of a structure within the rectangulardesign domain shown. The structureis required to carry a vertical load P, see Figure 8(A), and is offered support along the two hatchedpartsof the left handsideof the designdomain. The designdomain is subdividedinto 2610 finite elements,andthe solid volume fraction is given to be 16%. The result in Figure 8(A) canbe clearly interpretedasa trusswith a well-definedtopology in terms of number of joints and connectivetyof bars.We now perform shape opimizdon (i.e., optimization of sizing and configuration, cf. Section 5) of a truss with this topology, and obtain the result shown in Figure 8(B). Thus, subject to the given amount of available material, the given value of the force P and its distancesfrom the supportingjoints for the truss,the trussis not only optimal in termsof the numberandpositionsof joints, but also with respectto the connectivity andthe cross-sectionalareasof the bars. S. Optimization

of Fiber Composites

Optimization methods can be employed to tdur certain material propertiesof fiber/matrix composites by varying the stacking sequence,lamina thicknesses,orientation and density

On Optimum Design of Structures and Materials

155

(A) Figure8. (A) solution to topology optimizationproblem, (B) truss interpretationof the result with bar areas and positions of joints determined by sizing and configurationoptimization,respectively. of fibers in the lamina, etc. In this section, the problem of tailoring a linearly elastic fiber reinforced composite plate or single-ply lamina in plane stress (see [68]), shall be briefly discussed. The plate or lamina is assumed to have variable local orientation and density of two mutually orthogonal fiber fields embedded in the matrix material. The thickness and the domain of the plate or lamina is assumed to be given, together with prescribed boundary conditions and in-plane loading. The problem consists in determining throughout the structural domain the optimum orientations and densities of the fiber fields in such a way as to maximize the integral stiffness of the composite plate or lamina under the given loading. The optimization is performed subject to a prescribed bound on the total cost or weight of the composite that for given unit cost factors or specific weights determines the amounts of fiber and matrix materials in the structure. The mathematical formulation of the problem is identical to that of topology optimization, i.e. Equations (12)-(15), except that the expression for the variable tensor of elasticity Eqkl is different, and that a generalized form of the cost function in Equation (15) is adopted. As an example, consider a rectangular single-ply lamina whose domain is indicated in Figure 9(A). The lamina has one of its sides fixed against displacements in both the x and y directions while the opposite side is subjected to a parabolically distributed shear loading as indicated in Figure 9(A). The result of maximizing the integral stiffness of the lamina for given total amounts of the fiber and matrix materials is shown in Figure 9(B). Here, the direction and density of the hatching (which is shown on top of the applied finite element mesh) illustrate the fiber orientation and density, respectively. As is to be expected, the solution contains no fiber reinforcement in lowly stressed sub-domains, only one-way fibers in sub-domains with dominance of a single principal stress component, and a cross-ply fiber arrangement in subdomains with dominance of shear (i.e., two large principal stresses of opposite signs).

9. Design of Materials for Prescribed Elastic Properties The problem of tailoring materials with specified elastic properties has attracted considerable interest over recent years, and been considered by researchers from various fields (see, e.g., [1], [27], [31], [32], [5], [63]-[66], [28]).

136

Niels CUhoJ-

Figure 9. (A) Domain, loading, and boundary conditions for a composite single-ply lamina. (B) Optimum distri-

bution and orientation of fibers in the matrix material of the lamina.

Milton and Cherkaev[32] developeda mathematicalmethod for the designof materials which havean elasticity tensormatchinganarbitrarily prescribedpositive semi-definitefourth orderelasticity tensorcompatiblewith thermodynamics.A material with Poisson’sratio close to - 1 was constructedby Milton [31]. Thesedevelopmentswere basedon the useof layered materialmicrostructuresconstructedfrom a rigid andaninfinitely weakmaterial.The materials obtainedthis way serveto substantiatethat prescribedas well as extremeelastic properties canbe achieved,but they aremerely mathematicaltools thanpracticalcompositesdueto their widely differing length scales[27]. Sigmund [63], [66], on the other hand,treatsthe material designproblem as a topology optimization problemof a periodic microstructurerepresentedby a basecell that only involves a single length scale.The topology optimization problem is formulated as minimization of the concentrationof material in the basecell, subjectto equality constraintsthat expressthe prescription of the elastic properties.Some of the resultswill be presentedin the sequel.In Sigmund [66], the microstructureis discretizedby continuum-typefinite elementswhere the materialdensitiesin the individual elementsaretakenasdesignvariables.The effective elastic propertiesof the discretizedmicrostructurearefoundby a numericalhomogenizationmethod which leads to the definition of the material designproblem as an inverse homogenization problem.

157

On Optimum Design of Structures and Materials

l---CC) 1

r-1

I

0 i L--

-1

-I 2-J

I .75 0 (,q = o.o!l .75 I 0 0 0 .I25 1 t Q = 0.28 Figure 10. Honeycomb microstructures with Poisson’sratio equal to 0.75. Three different base cells, with (A)

21 x 12 elements, (B) 63 x 36 elements, and (C) 15 x 15 elements,are used.

The optimization problem can be statedas Minimize pcett= Jo pP dQ

0 < Pmin

5

P

5

Pmax

where pcell is the density of material in the basecell, the local material densitiesp are the designvariablesconstrainedby pminandpmax,fl denotesthe designdomain (basecell), andp is a penalizationexponentusedto ensurethat eachelementin thediscretizedcell will represent either solid or void. ,?Z$$denotesthe effective elasticity tensordeterminedby the standard numerical homogenizationprocedureas describedin [66] or [4], and E&[ is the prescribed elasticity tensor. The microstructurespresentedin the following [66] areall obtainedfrom rectangularbase cells discretizedby four nodebilinear finite elements,and the basematerial for all examples has Young’s modulus 0.91 andPoisson’sratio 0.3 suchthat a purely solid basecell will have Et t tr = 1.O(assumingplanestress). As a first example,the optimization algorithm is usedto ‘reinvent’ the ‘perfect’ honeycomb. Prescribing the elastic propertiesof a material with E 1111= 0.09 andPoisson’sratio 0.75 in a rectangularand quadratic domain for the basecell, the results in Figure 10 are obtained. The optimized microstructuresin Figs. 10(A),(B) are seento display great similarity with the ‘perfect’ honeycomb.If the rectangularbasecell domain in Figs. 10(A),(B) is changed to a quadratic one discretizedby 15 x 15 elements,the optimization yields the ‘octagonal honeycomb’ shownin Figure 10(C).As the densityof the materialis very nearly the samefor the ‘perfect’ andthe ‘octagonal’honeycomb,preferencemay be madefrom manufacturability considerations. As a secondexample, we considerthe designof materialswith negativePoisson’sratio, i.e., materials that expandtransverselywhen subjectedto elongation. It is well known that Poisson’sratio is confinedto the interval from 0 to 0.5 for solid materials,but for porous(e.g., cellular) materials Poisson’s ratio may approachthe value -1 without violation of positive semi-definitenessof the elasticity tensor.Negative Poisson’s ratio materials have attracted considerableinterestin recent years,and use of such materials may, for example, be advan-

158

Niels Olhoff

-] I

I

I

I

(A)

(!"i/

{E") -- 0.02 -.

0.

0 = 0.25

(B)

I 'l- 006 -0I - . 6

i)

0 = 0.59

Figure 11. Negative Poisson's ratio materials [(A): -0.8; (B): -0.6] obtained from base cells with 40x40 elements and enforcementof vertical symmetry.

i

3 mm

Figure 12. Testrig for experimental investigations of negative Poisson's ratio materials in Figure 11. The testrig has been producedusing micromachining techniques. (By courtesyof J.M. Guedes).

tageous for hydrophones and mechanical fasteners (e.g., 'Rawlplugs'). However, numerical experiments in [66] have shown that only low overall stiffness can be obtained if extreme elastic properties are prescribed. The reader is referred to [63] for potential applications and for a survey of recent research. Figure 11 (left) displays topologies of base cells obtained by [66] subject to specification of elastic properties of materials with a negative Poisson's ratio. The illustrations on the right hand side of Figure 11 are obtained by repeating the base cell periodically, and offer a clear picture of the aggregate microstructures; it is easily seen that horizontal elongation of the

On Optimum Design of Structures and Materials

159

microstructureswill give rise to expansionin the vertical direction, and thereby a negative value of the Poisson’sratio. It should be finally mentioned that a collaborative experimental project devoted to the manufacturing of the negative Poisson’s ratio materials in real micro-scale (size of base cell: 30 pm) and subsequenttesting of their elastic properties,has been initiated between 0. Sigmund and researchersof the Microelectronics Centerin Lyngby, Denmark, see [28]. Figure 12 shows a test rig with a test specimencomposedas an array of 30 x 30 pm2 basecells of the type shown in Figure 1l(A). The test rig, which is manufacturedby laser micro machining, is subjectedto four point bending, and elastic propertiesof the specimen are determinedby studying the displacementfield through a microscopeand measuringthe eigenfrequencyof the specimenfor different loadings. 10. Summary

Basic concepts of optimization were discussedin this paper, and a unified mathematical formulation presentedfor multicriterion design of structuresand materials. Fundamental aspectsand genericexamplesof optimum shape,sizing, andtopology designwere presented, and it has beenpointed out that use of topology optimization as a proprocessorfor shapeor sizing optimization allows thedesignerto arrive at muchbetterfinal results.The paperhasalso demonstratedthat optimization tools canbeusedto greatadvantagein the designof composite materials, and that the techniques open the possibility of tailoring material microstructures subjectto prescribedmacro-mechanicalproperties. Acknowledgements

This paper to a large extent reflects a coordinatedresearcheffort supportedby the Danish TechnicalResearchCouncil (Programmeof Researchon ComputerAided Design).The author is indebtedto many Danish and foreign colleaguesfor inspiring cooperation. References 1. Almgren, RF., ‘An isotropic three-dimensionalstructurewith Poisson’sratio = -l’, J. Elasticity, 15 (1985), 427-430. 2. Atrek, E. R., Gallagher, R.H., Ragsdell, K.M. and Zienkiewicz, O.C., (eds.), New Directions in Optimum Structural design, John Wiley and Sons,New York, 1984. 3. Barthelemy,B., Haftka, R.T., ‘Accuracyanalysisof the semi-analyticalmethod for shapesensitivityanalysis’, AIAA Paper 88-2284, Proc. AIAAlASME/ASCElASC 29th Structures, Structural Dynamics and Materials Conference, Part 1, 1988, pp. 562-58 1, Also Mechanics of Structures and Machines, 18 (1990), 407432. 4. Bendsoe, M.P., Optimization of Structural Topology, Shape, and Material. Springer-Verlag,Berlin, 1995. 5. Bendsoe, M.P., Dfaz, A.R., Lipton, R. and Taylor, J.E. ‘Optimal design of material properties and material distribution for multiple loading conditions, To appear in Int. J. Num. Meth. Engng. 6. Bendsqie,M.P. and Mota Soares,CA., (eds.), Topology Design of Structures, Kluwer, Dordrecht, The Netherlands, 1993. 7. Bendsoe, M.P., ‘Optimal shape as a material distribution problem’, Structural Optimization, 1 (1989). 193202. 8. Bendsoe, M.P. and Kikuchi, N., ‘General optimal topologies in structural design using a homogenization method’, Comp. Meth. in Appt. Mech. and Engrg., 71(1988), 197-224. 9. Bendsoe, M.P., Olhoff, N. and Taylor, J.E., ‘A variational formulation for multicriteria structural optimization’, J. Struct. Mech., 11(1983), 523-544. 10. Bennett, J.A. and Botkin, M.E., The Optimum Shape. Automated Structural Design. Plenum Press,New York, 1986. 11. Braibant, V. and Pleury,C., ‘Shapeoptimal designusing B-splines’,ComputerMethods in Applied Mechanics and Engineering, 44 (1984), 247-267.

160

Niels Olhofs

12. Cheng, G. and Liu, Y., ‘A new computational schemefor sensitivityanalysis’,Engineering Optimization, 12 (1987), 219-235. 13. Cinquini, C. and Rovati, M., ‘Problemsand methods of structuraloptimization’, In: Rovati, M. and Contra, R. (eds), Lecture notes for COMETS Course ‘Computer Aided Optimal Design of Structures - Advanced Materials and Shape Optimization’, September 1991, University of Trento, Italy. 14. Ding, Y, ‘Shape optimization of structures:a literature survey’, Computers and Structures, 24(b) (1986), 98.5-1004. Methods in Structural Design, B.I.15. Eschenauer, H.A. and Olhoff, N., (eds.), Optimization Wissenschaftsverlag,Mannheim, 1983. 16. Eschenauer,H.A., Mat&reck, C. and Olhoff, N., (eds.), Engineering Optimization in Design Processes. Springer-Verlag,Berlin, 1991. 17. Esping, B.J.D., Minimum WeightDesign ofMembrane Structures, Ph.D. Thesis,Dept. Aeronautical Structures and Materials, The Royal Institute of Technology, Stockholm, Report 83-1, 1983. 18. Esping, B.J.D., ‘Minimum weight design of membrane structuresusing eight node isoparametricelements and numerical derivatives’, Computers and Structures, 19(4) (1984), 591-604. 19. Fleury, C., ‘CONLIN: an efficient dual optimizer based on convex approximation concepts’, Structural Optimization, 20.

Journalfor 21.

1 (1989), 81-89.

Fleury, C. and Braibant, V., ‘Structural optimization: a new dual method using mixed variables’,International Numerical

Methods in Engineering,

23 (1986), 409428.

Gilmore, B.J., Hoeltzel, D.A., Azram, S. and Eschenauer,H.A., ‘Advancesin design automation’, Proc. 1993 ASSME Design Technical Conferences - 19th Design Automation

22. 23. 24. 25. 26. 27.

Conference, Albuquerque,

New Mexico,

The American Society of Mechanical Engineers,New York, 1993. Haftka, R.T. and Grandhi, R.U., ‘Structural shapeoptimization - a survey, Computer Methods in Applied Mechanics and Engineering, 57 (1986), 91-106. Haug, E.J., (ed.), Concurrent Engineering: Tools and Technologiesfor Mechanical System Design, SpringerVerlag, Berlin, 1993, pp. 998. Haug, E.J. and Cea, J., (eds.), Optimization of Distributed Parameter Structures, Sijthoff and Nordhoff. The Netherlands, 1981, vol. 1, pp. l-842, ~012, pp. 843-1609. Haug, E.J., ‘A review of distributed parameter structural optimization literature’, In: Hang, E.J. and Cea, J. (eds), Optimization of Distributed Parameter Structures, vol. 1, Sijthoff and Nordhoff. The Netherlands, 1981, pp. 3-74. Herskovits,J., (ed.), Proc. Structural Optimization ‘93-The World Congresson Optimal Design of Structural Systems1 - 2. Federal University of Rio de Janeiro,Brazil, 1993. Lakes, R., ‘Design considerationsfor materials with negative Poisson’sratio’, J. Mech Design, 115 (1993), 696-700.

Larsen, U.D., Sigmund, 0. and Bouwstra, 8, Fabrication of Compliant Micromechanisms and Materials with Negative Poisson> Ratio. Manuscript, 1995. 29. Lund, E., Finite Element Based Design Sensitivity Analysis and Optimization, Ph.D. Thesis, SpecialReport No. 23, Institute of Mechanical Engineering, Aalborg University, Denmark, 1994, p. 240. 30. Lund, E. and Olhoff, N., ‘Shapedesign sensitivity analysisof eigenvaluesusing ‘exact’ numerical differentiation of finite element matrices’, Structural Optimization, 8 (1994), 52-59. 31. Milton, G.V., ‘Composite materials with Poisson’sratio close to -l’, J. Mech. Phys Solids, 40(5) (1992), 1105-l 137. 32. Milton, G.V. and Cherkaev, A.V., ‘Materials with elastic tensorsthat range over the entire set compatible with thermodynamics’, Presentedat the joint ASCE-ASME-SES MettN, June 1993, University of Virginia, Charlottesville, Virginia, USA. 33. Morris, A.J., (ed.), Foundations of Structural Optimization: A UntfiedApproach. Chichester,England, 1982. 34. Mota Soares,CA., (ed.), Computer Aided Optimal Design: Structural and Mechanical Systems. SpringerVerlag, Berlin, 1987. 3.5. Olhoff, N., ‘Multicriterion structural optimization via bound formulation and mathematical programming’, 28.

Structural 36. 37. 38.

Optimization,

1 (1989), 11-17.

Olhoff, N., Bendsoe,M.P. and Rasmussen,J., ‘On CAD-integrated structural topology and design optimization’, ComputerMethods in Applied Mechanics and Engineering, 89 (1991), 259-279. Olhoff, N. and Lund, E., ‘Finite element basedengineering design sensitivity analysisand optimization’, In: Herskovits, J. (ed.) Advances in Structural Optimization, Kluwer, The Netherlands, 1995, pp. l-45. Olhoff, N., Lund, E. and Rasmussen,J., ‘Concurrentengineeringdesign optimization in a CAD environment’, In: E.J.Hang (ed), ConcurrentEngineering: Toolsand TechnologiesforMechanicalSystemDesign, SpringerVerlag, Berlin, 1993, pp. 523-586.

On Optimum Design of Structuresand Materials

161

39. Olhoff, N., Rasmussen,J. and Lund, E., ‘A method of ‘Exact’ numerical differentiation for error elimination in finite element based semi-analytical shape sensitivity analysis’,Mechanics of Structures and Machines, 21(1993),

l-66.

40. Olhoff, N. and Rozvany, G.I.N., Structural and Multidisciplinary Optimization, Pergamon,Oxford, 1995. 41. Olhoff, N. and Taylor, J.E., ‘On structural optimization’, Journal of Applied Mechanics, 50 (1983), 11391151. 42. Papalambros,F?Y.and Chirehdast, M., ‘An integrated environment for structural configuration design’, J. Engrg. Design. 1 (1990), 73-76.

43. Pedersen,P., ‘On the minimum masslayout of trusses’,Proc. AGARD Symposium 1970, Istanbul, AGARDCP-367., 1970. 44. Pedersen,P., ‘On the optimal layout of multi-purpose trusses’, Computers and Structzwes. 2 (1972), 695-7 12. 45. Pedersen,P., ‘Optimal joint positions for spacetrusses’,J. Struct. Div., ASCE 99 (1973), 2459-2476. 46. Pedersen,P., (ed.), Optimal Design with Advanced Materials. Elsevier,The Netherlands, 1993. 47. Pedersen,P., ‘Concurrent engineeringdesign with and of advancedmaterials’,In: E.J.Hang (ed.), Concurrent Engineering: Tools and Technologiesfor MechanicalSystem Design, Springer-Verlag,Berlin, 1993, pp. 62748. 49.

670. Pedersen, P., ‘Bounds on elastic energy in solidsof orthotropic materials’, Structural Optimization, 2 (1990), 55-63. Pedersen,P., ‘On optimal orientation of orthotropic materials’, Structural Optimization, 1 (1989), 101-106.

50. Pedersen,P., Cheng, G. and Rasmussen,J., ‘On accuracyproblems for semi-analyticalsensitivity analysis’, Mechanics of Structures and Machines, 17(3) (1989), 373-384.

51. Pedersen,P., ‘Design with severaleigenvalueconstraintsby finite elementsand linear programming’, Journal of

Structural

Mechanics, lO(3) (I 98 1), 243-271.

52. Rasmussen,J., Development of the Interactive Structural Shape Optimization System CAOS (in Danish), Ph.D. thesis, Institute of Mechanical Engineering, Aalborg University, Denmark, Speck1 Report No. la, 1989. 53.

54. 55. 56. 57.

Rasmussen,J., ‘The structural optimization systemCAOS’, Structural Optimization 2 (1990), 109-l 15. Rasmussen,J., Lund, E., Birker, T. and Olhoff, N., ‘The CAOS system’, International Series of Numerical Mathematics, 110 (1993), 75-96, BirkhauserVerlag, Basel. Rasmussen,J., Lund, E. and Olhoff, N., ‘Integration of parametric modeling and structural analysis for optimum design’, Gilmore et ul. (eds), Proc. Advances in Design Automution, The American Society of Mechanical Engineers,Albuquerque, New Mexico, USA, 1993. Rasmussen,J., Olhoff, N. and Lund, E., ‘Foundations for optimum design systemdevelopment’, In: Lecture Notesfor: Advanced TEMPUS Course on NumericalMethods in Computer Optimal Design, 1l-l 5 mai 1992, Zakopane, Poland, 1992. Rozvany, G.I.N. and Karihaloo, B.L., (eds.),Structurul Optimization, Kluwer, Dordrecht, The Netherlands, 1988.

58.

59. 60. 61. 62. 63. 64.

Rozvany,G.I.N., (ed.), OptimizationofLargeStructuralSystems, Kluwer, Dordrecht,TheNetherlands, 1993, vol. I, pp. l-622, vol. 2, pp. 623-l 198. Rozvany, G.I.N., Bendsoe, M.P. and Kirch, U., ‘Layout optimization of structures’, Applied Mechanics Reviews. 48 (1995), 41-l 19. Sacchi Landriani, G. and Rovati, M., ‘Optimal design for two-dimensional structuresmade of composite materials’, J. Engineering Materials und Technology, 113 (1991), 88-92. Schmit, L.A., ‘Structural synthesis- its genesisand development’, A1AA Journal, 20 (1982), 992-1000. Sergeyev,0. and Pedersen,P., On design ofjoint positions for minimum mass 3Dfiames Manuscript, 1995 p. 21. Sigmund, O., ‘Design of Material Structures Using Topology Optimization’, DCAMM Special Report No. 69, Ph.D. Thesis,Technical University of Denmark, 95 p. +22 p. (appendix), December 1994. Sigmund, O., ‘Materials with prescribedconstitutive parameters:an inverse homogenization problem’, Znt. J. Solids Structures, 31(17) (1994), 2313-2329.

Sigmund, O., ‘Tailoring materials with prescribedelastic properties’, Mech. Muterials, 20 (1995), 351-368. 66. Sigmund, O., ‘Design of material structuresusing topology optimization’, In: Olhoff, N. and Rozvany, G. (eds.) Structural andMultidisciplinary Optimization, Pergamon,Oxford, 1995. 67. Svanberg,K., ‘The method of moving asymptotes- a new method for structural optimization’, Internutional Journalfor Numerical Methods in Engineering, 24 (1987), 359-373. 68. Thomsen, J. andolhoff, N., ‘Optimization of fiber orientation and concentrationin composites’, Control and

65.

Cybernetics, 19 (1990), 327-341.

69. Zienkiewicz, O.C. and Campbell, J.S., ‘Shapeoptimization and sequentiallinear programming’, In: Gallagher, R.H. and Zienkiewicz (eds.),Optimum StructuralDesign, Theory andApplications, Wiley and Sons,London, 1973, pp. 109-126.