INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 2000; 48:185–209 A simple approach towards fully stressed design in hyperelasticity exploiting the isoparametric nite element concept§ Eckhard Kirchner1; ∗; †; ‡ and Rainer Immel2 1 2 Department for Machine Elements and Machine Acoustics; Darmstadt University of Technology; Magdalenenstrasse 4; 64289 Darmstadt; Germany Adam Opel AG; International Technical Development Center; IPC 81-90; 65423 Russelsheim; Germany SUMMARY A novel approach towards fully stressed designs in hyperelasticity is discussed leading to closed-form expressions for the sensitivities of the objective and displacements with respect to design variations. The key idea is the modication of the classical approach coupled with a so-called design element method oering a lot of parallelism to standard nite element methods. We bypass implicit constraints on dependent quantities and derive an explicit linearly constrained optimization problem solved by means of rst-order procedures. The results obtained with the proposed method are adequate from an engineering point of view though being computed with a simple method. Copyright ? 2000 John Wiley & Sons, Ltd. KEY WORDS: shape optimization; sensitivity analysis; fully stressed design 1. INTRODUCTION In the last years a sequence of approaches towards optimal shape design has been published. The majority is either devoted to problems with small deformations or to elastic behaviour. For engineers it is more and more desirable to include all kinds of non-linearities in computational optimization processes to satisfy the requirements addressed to new designs. Large deformations and inelastic material behaviour shall be considered to achieve both light weight and structural rigidity and reliability. These aims are frequently contradictory. Usually, engineers do not want their designs to be plastically deformed. But when a certain behaviour under extreme (inelastic) loading is desired, the structural optimization has to be amplied to conciously account for inelastic problems associated with large deformations. As a rst step ∗ Correspondence to: Eckhard Kirchner, Department for Machine Elements and Machine Acoustics, Darmstadt University of Technology, Magdalenenstrasse 4, D-64289 Darmstadt, Germany † E-mail: eckhard.kirchner@de.opel.com ‡ Current address: Adam Opel AG, International Technical Development Center, Powertrain Simulation, IPC 81-34, 65423 Russelsheim, Germany § Dedicated to Professor Franz G. Kollmann on the occasion of his 65th birthday Copyright ? 2000 John Wiley & Sons, Ltd. Received 9 March 1999 Revised 29 September 1999 186 E. KIRCHNER AND R. IMMEL we therefore transform the implicitly constrained classical approach towards fully stressed linear elastic structures into an explicitly constrained optimization problem. We thus need to allow for large hyperelastic deformations associated to stresses larger than the classical limits from linear engineering analysis. In such extremely loaded states it is then only of minor importance whether customary limit states are exceeded or not. This fact motivates us to modify the classical approach of fully stressed designs, cf. Wei et al. [1] or Wu and Arora [2], to overcome the (algorithmic) shortcomings of the ordinary stress constrained approaches. We propose a simple modication of the classical approach for fully stressed designs and show that it is convenient from an engineering point of view, and also with regard to a manageable realization within the framework of the nite element method. Our method is not only applicable to problems with unsharp yield limits, but also easy to implement and requires only a limited number of changes in the nite element program to be run successfully. However, we aim at applications to inelastic materials with rate-dependent eects where a well-dened yield limit does not exist. For details concerning applications to inelastic anisotropic materials we refer to the paper by Kirchner [3]. But here we restrict ourselves to hyperelasticity and focus on the modication of the classical approach and some algorithmic details. The general treatment of inelastic rate-dependent problems is shown, for instance, by Tsay and Arora [4] in a variational sensitivity analysis context without addressing special objectives or respective constraints. Mahnken and Stein [5] present details on parameter identication methods using nite elements. The concept outlined in this work is well applicable to viscoplasticity. In these cases, the sensitivity analysis covers also the computation of recursively dened expressions for the sensitivities of the stress eld according to the chosen time integration scheme. Here we restrict ourselves to giving simple but nevertheless exact expressions for hyperelastic problems. Nowadays, shape optimization follows two dierent paths. The rst one addresses the classical eld of form nding problems, often referred to as shape optimization. The second branch is referred to as topology optimization. An example for a combination of both branches is given by Ramm and co-workers [6; 7]. Shape optimization focuses on the shape of a body itself often including cross-sectional parameters or element thicknesses as design variables. By contrast, topology optimization considers the material distribution as the major problem. It is obvious that a clear distinction between these two branches is not available in a general form. Furthermore, Barthold et al. [8; 9] rely on the parametrization of the boundary contours using design nodes, -lines and -surfaces while dening the interior of the body independent of the boundaries. In addition, a hierarchy of design variables may be introduced to derive a model-adaptive type shape optimization procedure (see e.g. References [9–11]). Transferring and extending concepts of the nite element method—or more precisely its basic ideas—yields a simple shape optimization procedure. Kirchner [12] presents a chart focusing on the parallels between the nite element method and the so-called design element method which is utilized herein. Our motivation is to demonstrate an ecient way to achieve highly loaded designs in nite strain hyperelasticity which does not require the application of mathematic programming procedures. Numerous numerical schemes in that eld are devoted to the solution of constrained and unconstrained optimization problems, among them are for instance the methods by Schittkowski [13] or Bertsekas [14]. Most of the recent approaches are based on sequential linear and quadratic approximations (SLQ and SQP) of the objective function or its associated Lagrangian if constraints are included. (We refer for example to Wei et al. [1] or Reitinger et al. [15]. The majority of modern approaches uses some approximation of the Hessian of the objective function, or of a suitable Lagrange function. For simplicity, we use in this paper a rst-order Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 ISOPARAMETRIC FINITE ELEMENT CONCEPT 187 procedure of predictor-corrector type coupled with an adaptive step length estimate in the design parameter space of the optimization problem. This adaptive step length estimate is in parallel to the method of moving asymptotes (MMA) by Svanberg [16] or Bletzinger [17] which adapts the step length in accordance to suitable approximations of the objective function. A rigorous computation of the required sensitivities in a framework of isoparametric nite elements enables us to establish closed-form expressions reducing the additional numerical eort. We do not apply numerical approximation schemes such as nite dierences, rather we use explicit expressions only. This paper is structured as follows: subsequent to the motivation and introduction we review some results from continuum mechanics and give a brief description of the nite element method applied. In Section 3 we explain an optimization approach towards highly loaded hyperelastic structures and discuss the objective functions as well as the type of constraints under consideration. Also, analytic discrete expressions for the sensitivities are given. The applicability of the proposed method is shown by three examples in Section 4. Finally, we draw some conclusions. 2. CONTINUOUS AND DISCRETE DESCRIPTION OF FINITE STRAIN HYPERELASTICITY We briey review the most important items of continuum mechanics needed for an appropriate description of the optimization problem discussed here. Furthermore, we give the information about the nite element scheme needed. A fully non-linear notation is used throughout. The following index convention is applied: Plain roman indices, such as i, denote usual Cartesian components of vectors and tensors. Einstein’s summation convention is assumed on repeated indices from 1 to ndm , with space dimension ndm of the Euclidian space Rndm . Indices in parenthesis ( j) are used for iteration counters, roman letters in square brackets [l] , [m] (l; m = 1; : : : ; nnle ) and [L] (L = 1; : : : ; nng ) refer to local and global nodal quantities in the nite element method. The number of element nodes of a typical structural element e is nnle , the total number of nodes of the structure is nng . 2.1. Continuum mechanics background Consider a continuous body B in a Euclidean space Rndm which deforms under loading. Let e be the motion in time t x = e(X; t) (1) As usual the vector X denotes the reference conguration of a material point, the vector x is its current position, x = X + u, where u is the displacement vector. Using convected curvilinear co-ordinates i , i = 1; 2; 3, with associated (contravariant) base vectors Gi the deformation gradient F is F = Grad e = Grad x = @x ⊗ Gi @i (2) dening a mapping of (innitesimal) material line elements from the reference to the actual placement. In the sequel, we use the following transformation of the gradient operator from the material (Grad) to the spatial (grad) setting, i.e. grad (•) = Grad (•) F−1 Copyright ? 2000 John Wiley & Sons, Ltd. (3) Int. J. Numer. Meth. Engng. 2000; 48:185–209 188 E. KIRCHNER AND R. IMMEL Furthermore, we need the Gateaux derivative [18] of the deformation gradient, D F, with respect to an incremental change u of the displacement eld u D F · u = Grad u · u ⇒ F(u) = Grad u (4) For the description of the material behaviour we quantify the deformation of a material point in terms of objective strain measures being invariant under superposed rigid-body motions. For convenience, we choose the symmetric left and right Cauchy Green tensors b = FFT and C = FT F (5) For the derivation of an appropriate nite element formulation we need the variation and linearization of the right Cauchy Green strain C expressed in terms of variations of the displacement or incremental displacements, w and u, respectively C(u) = 12 (FT F + FT F) = Sym ((Grad u)T F) (6) C(w) = 12 (FT F + FT F) = Sym ((Grad w)T F) (7) Finally, the linearization of the variation (or vice versa) of C reads C(u; w) = Sym ((Grad u)T Grad w) (8) We introduce two symmetric contravariant stress measures, the Kirchho and the second Piola Kirchho stress tensors = ij gi ⊗ gj and S = S ij Gi ⊗ Gj (9) refering to instantaneous and reference conguration, respectively. Both are related to each other by means of push-forward and pull-back operations [18; 19] S = ∗ () = F−1 F−T (10) For our computations we use simple non-linear hyperelastic material models of the Ogden [20] family; for details we refer to Reese [21]. It is well known that use of a spatial formulation is advantageous with regard to an ecient nite element formulation for these material models. We use the following Neo-Hookean material model in a spatial framework = 2 (J − 1)i + (b − i) 2 (11) where J is the volumetric deformation J = det(F). In (11) we use a Cartesian unit tensor i of rank two, and are Lame parameters, cf. (52). For the example presented in Section 4.3 an Ogden type model is applied. It is derived from the following strain energy function W, expressed in terms of principal stretches i2 associated to the left Cauchy Green tensor b, i = 1; 2; 3 3 X r r r r (1 + 2 + 3 − 3) − r ln J + (J 2 − 1 − 2 ln J ) (12) W= r 4 r=1 Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 ISOPARAMETRIC FINITE ELEMENT CONCEPT 189 The material parameters i and i , i = 1; 2; 3, are given in Section 4.3 for a rubber-like incompressible material having a large bulk modulus compared to the shear-type moduli i . The Kirchho stress tensor follows from (12) with a Doyle Erickson-type relation =2 @W b @b (13) Since the stored energy function (12) is dened in terms of the principal stretches, the following relation using the eigenvector basis ni? , i = 1; 2; 3, of the left Cauchy Green tensor b (which coincide with the eigenvectors of by the restriction to isotropy) for the stress computation is more convenient = 3 @W P ni? ⊗ ni? i=1 @i (14) For specic details of this model concerning the derivation of the associated tangent arrays we refer to Reese [21] or Barthold [8]. Anyway, the i are denoted as principal stretches in this paper since they are the principal values from right and left stretch tensor following from the polar decomposition theorem applied to the total deformation F (see e.g. [22]). 2.2. Weak form of equilibrium und linearization Here, the nite element method applied to geometrically non-linear hyperelasticity is derived from the weak form of equilibrium in its standard form Z Z Z 1 t̂ · w dA − S : C(w) dV − 0 b · w dV = 0 (15) G= 2 B0 @B0t B0 with C(w), cf. (7), being the variation of the Cauchy Green strain due to a variation of the displacement eld w, also interpreted as a virtual displacement. The part of the boundary of B with given tractions t̂ is denoted as @B0t , b is the body force per unit mass. Application of the concept of the Gateaux derivatives with respect to incremental displacements u to the weak form of equilibrium yields, according to (15) and (6), (7) as well as (8), the material and spatial linearizations 1 DG · u = 2 = 1 2 Z B0 Z dS : C + S : C dV C : dC (16) d : sym (grad w) db B0 + sym (grad u) : sym (grad w ) dV sym (grad u) : (17) Use of a material framework for the computation of the linearization in (16) is advantageous because its push forward (17) leads to a simple and ecient element formulation with the same structure as in the geometrically linear regime [23]. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 190 E. KIRCHNER AND R. IMMEL 2.3. Discrete nite element formulation We discretize the reference domain B0 by splitting it into nel nite elements, nel [ @B0; i ∩ @B0; j neighbouring elements i and j B0; e with B0; i ∩ B0; j = B0 = ∅ else (18) e=1 and use a standard isoparametric framework for two- or three-dimensional elements. Then, we introduce discrete approximations uh and Xh of the displacement eld and the reference co-ordinate uh = nP nl e l=1 N[l] (W) u[l] and Xh = nP nl e l=1 N[l] (W) X[l] (19) with u[l] and X[l] being the nodal displacements and co-ordinates, N[l] are shape functions of Lagrangian type for a typical element e, e = 1; : : : ; nel . For the computation of both structural response and sensitivities we need the Jacobian Jeiso of the isoparametric mapping from the parent domain with a co-ordinate W to the reference conguration of an element e with nnle nodes, Jeiso = nl e dN[l] dXh nP = X[l] ⊗ dW dW l=1 (20) usually = [−1; 1] × [−1; 1] for two-, and = [−1; 1]3 for three-dimensional volume elements. Using (19) and (20), the material gradient operator appears as d dW d d d d (Jiso )−1 ⇔ = = = (J iso )−1 dXh dW dXh dW e dXih dj ji e (21) Consequently, the discrete deformation gradient Fh follows from (2) and (19) Fh = i + Grad uh = i + nP nl e l=1 h u[l] ⊗ nP nl e dN[l] dN[l] iso −1 h (J ) =i + u[l] ⊗ h dX dW e l=1 (22) Application of standard arguments [24; 25] leads to the following discrete balance equation derived from (15): Pi ((u); u) − Pa (t) = 0 (23) with Pi and Pa being vectors of internal and external nodal forces, respectively, written as an assembly of the contributions of all nel elements P = i Pa = nel S Z e=1 nel S Z e=1 (ue ) Be (ue ; Xe ) det(Jeiso ) dW (24) 0 NeT b det(Jeiso ) dW (25) In (24), Be is the discrete strain displacement operator for element e which employs the shape functions Ne and represents the symmetric spatial gradient operator according to (3) and applied to the displacement eld u and its variation w in (17). The shorthand notations ue and Xe represent all nodal displacements u[l] and co-ordinates X[l] , l = 1; : : : ; nnle appearing in (19). The integral over the parent domain is computed numerically by means of Gauss quadrature rules. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 ISOPARAMETRIC FINITE ELEMENT CONCEPT 191 Figure 1. Classical objective with stress constraint and self-limiting modication for fully stressed designs. To solve the non-linear system of algebraic equations (23) Newton’s method is applied using the consistent linearization [18] of the balance equation, i.e. " u(hj+1) = u(hj) − #−1 dPi def : (P(ij) − Pa ) = u(hj) − (K(Tj) )−1 (P(ij) − Pa ) du ( j) (26) where ( j) denotes the iteration index. Equation (26) is solved repeatedly until a certain tolerance value TOL is attained, P(ij) − Pa( j) 6TOL for a suitable maximum or Euclidean norm. In (26) we assume conservative loads as in (15), i.e. Pa 6= Pa (u). 3. AN ALTERNATIVE APPROACH TO FULLY STRESSED DESIGNS We discuss a simple modication of the classical approach of fully stressed design leading to an optimization problem which does not require the introduction of implicitly dened constraints formulated in terms of (dependent) displacements or stresses. Furthermore, we discuss the computation of the gradient of the objective function based on an isoparametric nite element concept. The equations given herein apply to both two- and three-dimensional nite element formulations. 3.1. Modication of the classical objective function So far, most approaches to obtain fully stressed designs have been based on the following objective function class , which relates an equivalent stress eq to a given limit stress y . This function is illustrated in Figure 1 and referred to as ‘classical’ in this paper ( 2 ) Z eq (Y) 1− dv (27) min [class (Y)] with class = y B(Y) For simplicity, only small strains are considered and the notation from the geometrically linear regime is applied. The vector Y contains the nodal design parameters introduced in the sequel. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 192 E. KIRCHNER AND R. IMMEL The integrand of the objective functions such as in (27) is also denoted as specic quality in this paper. To name only a few examples for the objective function (27) we refer to Wei et al. [1] or Wu and Arora [2]. To avoid rapidly growing equivalent stresses eq → ∞, and thus → − ∞, it is necessary to introduce additional constraints on the stress eld g(Y) = eq (u(Y)) − y 60 ∀x ∈ B(Y) (28) The number of quadrature points exceeding the limit (28) in a computation generally varies from iteration to iteration. It is, therefore, necessary to take into account this variation by means of so-called active constraints of the optimization problem. For each quadrature point that violates a constraint g¿0 a dual variable is introduced in the Lagrangian associated to the objective (27). The numerical and computational disadvantages are obvious. Having applications of viscoplasticity in mind the following modication mod of the classical objective function (27) is applied leading to a locally optimal solution with vanishing variations “near” a well-dened target stress t , Z mod (Y) = B(Y) 1− eq t 2 !2 Z dv = B(Y) ( 1−2 eq t 2 + eq t 4 ) dv (29) For a single material point, both objectives functions (27) and (29) are equivalent for a target stress equal to the yield limit, y = t = eq . However, a gradient-based optimization method would increase the equivalent stress eq locally for the classical objective function (27) to obtain a smaller value of the objective function. This is due to the negative gradient pointing in the direction of fastest decrease of the objective function. These shortcomings are avoided by using the modied objective function (29). From an engineering point of view, this modication makes sense if local stresses above a target stress t are permitted, resulting in small plastic deformations of ductile materials when t = y . Hereby, the strengthening eect due to the development of plastic zones is exploited to increase the structural load carrying capacity. The objective functions (27) and (29) can—by the restriction to linear isotropic elasticity—be interpreted mechanically as linear and quadratic forms of the distortional energy density W dist = 2 eq =(2), respectively. The term 1 − (eq =t )2 represents a specic distortional energy deviation compared to the yield or target stress distortional energy density, 1− 2 2 t2 =(2 ) − eq =(2 ) Wtdist − W dist (eq (Y)) eq = = 2 t t =(2 ) Wtdist (30) Such an interpretation of both objective functions is in agreement with considerations of the distortional energy also used in standard J2 -plasticity models, cf. Lubliner [22]. Hence, it is obvious that constraints for the classical objective function (27) being a linear function of the energy deviation (30) are necessary to avoid innitely large negative deviations from the target state for large equivalent stresses. On the other hand, the quadratic form of the distortional energy deviation (30) corresponding to the modied objective function (29) is self-limiting (see Figure 1) and thus treatable without additional implicit constraints. Expanding the simple one-dimensional formulation (29) to a three-dimensional nite strain approach, we replace the equivalent engineering stress eq by the Kirchho stress and obtain the Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 193 ISOPARAMETRIC FINITE ELEMENT CONCEPT following objective function to be minimized: ( 2 ) Z Z : Pvm : : Pvm : def : 1−2 dv = + = 2d 2d B(Y) B0 (Y) J dV (31) Here, the target equivalent stress is d . The Kirchho equivalent stress of von Mises type is obtained by Pvm = 23 (I − 13 i ⊗ i) with I being the unit tensor of rank four. The computation of the objective function by integrating over the reference conguration creates additional terms in the gradient of the objective function representing the sensitivities of the volumetric deformation. It should be mentioned that, for time-dependent problems, the reference conguration can be used to formulate incremental recursive formulas for the stress computation using the material co-ordinate system rather than convected spatial frames. This oers some advantages when calculating the stress sensitivities for history-dependent constitutive models. However, explicit constraints which are expressed in terms of the independent optimization variables Y should be added to have certain bounds for the design node positions. We restrict the discussion here to linearly explicit constraints dening upper and lower limits for the nodal co-ordinates by Ymin 6Y6Ymax (32) The numerical treatment of these conditions is based on a projection procedure known as the method of feasible directions [26; 27]. In the case that volume conservation constraints need to be included additonally, a penalty type approach can be used as published by Kirchner [12]. Formally, a penalty-based approach considering the nodal position bounds (32) and corresponding volume constraints represents a least-squares method. 3.2. Design element approach As mentioned before, we want to apply the basic ideas of the nite element method to establish a modular shape optimization routine. In particular, the following characteristic properties of the nite element method can be directly utilized in shape optimization problems: Variational principle: The minimization of a hyperelastic potential transfers to the minimization of the objective function. Solution strategy: Newton or quasi-Newton methods are applied to solve the global equilibrium condition (23) representing the minimization of a hyperelastic potential while generalized xed point iteration methods are used to nd stationary points of the objective function with vanishing gradients. Discretization approach: Subdivision of the structure into nite elements and assembly of element contributions for structural analysis is in parallel to the shape approximation with design elements and corresponding assembly of design element contributions for objective function and gradient computation. Kirchner [12] gives a more comprehensive overview of the parallelisms between the nite element method and the design element approach in shape optimization applied herein. However, to explain the shape description we introduce an overall vector of shape parameters Y covering the nodal co-ordinates Y[A] of all design nodes A = 1; : : : ; ndn by T T T ; Y[2] ; : : : ; Y[n } YT = {Y[1] dn ] Copyright ? 2000 John Wiley & Sons, Ltd. (33) Int. J. Numer. Meth. Engng. 2000; 48:185–209 194 E. KIRCHNER AND R. IMMEL Figure 2. Design element concept: (a) reference domain for the description of the design element geometry, (b) design element in two-dimensional space being subdivided into individual nite elements. The position of the structural nodes X[L] , L = 1; : : : ; nng , is computed by a formula similar to (19) that are as demonstrated in Figure 2. The design element is provided with shape functions N[A] (W) dened in an appropriate reference domain using non-dimensional co-ordinates W: X[L] = ndn P A=1 Y[A] N[A] (W[L] ) (34) [A] and are mapped onto Nodes in the reference domain W of the design element are denoted by H Y[A] with ndn being the number of design nodes of a Lagrangian type design element. Dierentiation of (34) and assembly of the respective contributions of all design nodes yields a transition relation. This relation transforms the derivatives with respect to (global) design variables from the (ndm · ndn )-dimensional design space to locally computable derivatives with respect to the element node positions in Euclidian space R ndm . It is given by nS nS nS nS ng ng dn dn d d dX[L] d = = N[A] (W[L] ) dY A=1 L=1 dX[L] dY[A] A=1 L=1 dX[L] (35) The transition to the elementary nodes (35) is independent of actual design Y( j) during the optimization iteration with index j. It depends solely on the choice of the reference nodal points W[L] in (34) and is, therefore, economically computed only once in a session. Bletzinger [11] and Falk [9; 10] introduce a hierarchy of approximations using one-, two-, and three-dimensional design parameters (nodes–lines–surfaces) which does not oer the advantages of the simple method discussed above. 3.3. Sensitivity analysis using isoparametric nite elements In nite element computations the objective function (31) can be computed in analogy to the chosen spatial discretization (18) by integration over all element domains. Hence, = nel Z [ {(u(Y)); d } det{F(u(Y); Xe )} det Jeiso (Xe ) dW (36) e=1 Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 195 ISOPARAMETRIC FINITE ELEMENT CONCEPT where the quality per unit volume is dened in (31). In (36) we already pointed towards the dierent types of dependencies of the entries. The nodal co-ordinates Xe on an element e depend explicitly on the design variables Y via (34) while the displacement eld u is, as to be seen in the sequel, only implicitly dependent on design in a global sense. Thus we use the notation det(F) = det{F(u(Y); Xe )}, etc. The discrete gradient ∇ = d=dY of the objective function follows from (36) by total dierentiation with respect to the vector of design parameters Y, cf. (33). Applying product and chain rule yields @ d @ du = + dY @u dY @Y # "n Z el [ du iso @ det(Je ) { {(u(Y)); d } det{F(u(Y); Xe )}} dW = @u dY ∇ = e=1 + nel Z [ e=1 " @ { {(u(Y)); d } det{F(u(Y); Xe )} det{Jeiso (Xe )}} dW @Y nel Z [ det(Jeiso ) = det(F(u; Xe )) e=1 @ d @ du # @ + (; d ) {det(F(u; Xe ))} @u {z | dW } du dY global contribution nel Z [ d det(Jeiso ) (; d ) det(F(u; Xe )) + dY e=1 | @ {det(F(u; Xe ))} + (; d ) det(Jeiso ) @Y {z dW } (37) local contribution with ∇ ∈ R N and N being the total number of design parameters. The derivatives in (37) are computed by total dierentiation with respect to the element nodal co-ordinates exploiting the connection (35) between derivatives with respect to design variables Y and element parameters X[L] , L = 1; : : : ; nng . Afterwards, the element contributions are assembled accordingly. The local contributions to the gradient in (37) depend on the displacement eld u only via the stress eld and the total deformation F. They are computable locally for each element without having to account for the displacement sensitivities of the global part. Making use of the denition of the isoparametric map (20) we have the sensitivity of the element volume dN[l] @ det(Jeiso ) = det(Jeiso ) @X[l] dX Copyright ? 2000 John Wiley & Sons, Ltd. (38) Int. J. Numer. Meth. Engng. 2000; 48:185–209 196 E. KIRCHNER AND R. IMMEL with d(det(A))=dA = det(A) A−T , cf. Gurtin [28]. Furthermore, we need the derivative of the inverse Jacobian of the isoparametric map (20) with respect to the nodal co-ordinates iso d(Jijiso )−1 dJop dN[l] iso −1 = −(Jioiso )−1 (Jpj ) = −(Jikiso )−1 dXk [l] dXk [l] dXj (39) Using (39) we obtain the sensitivity of the volumetric deformation with respect to a variation of the element geometry dN[m] @(Jeiso )−1 u[m] ⊗ dW @X[l] m=1 dN [l] = −det(F h ) (Grad u h )T (F h )−T dX @ det(F h ) d det(F h ) = : @X[l] dF h nnl Pe (40) which completes the local contribution of the gradient of the objective function. In a similar way we obtain the derivative of the volumetric deformation with respect to the nodal displacements dN[l] @ det(F h ) = det(F h ) (F h )−T @u[l] dX (41) Since we assume a spatial element formulation throughout this paper it is convenient to interpret the term (F h )−T dN[l] =dX in (40) and (41) by comparison with (3) as spatial derivatives of the respective shape functions N[l] already known from the tangent stiness computation. The derivative d =d follows directly from (31) and the expression d=du equals the consistent tangent operator of the constitutive model. For the calculation of the displacement sensitivity we dierentiate the vector-valued balance equation of internal and external nodal forces (23) globally with respect to the design vector d [Pi {(u(Y)); u(Y); Y} − Pa {Y}] dY = @Pi (; u; Y) + @Y @Pi (; u; Y) @Pi (; u; Y) d du @Pa (Y) + − =0 @u @ du dY @Y | {z } KT (42) Here, we can identify the expression in brackets as the tangent stiness matrix KT already computed and inverted during the equilibrium iteration, see (26). Hence, it covers all changes of the internal forces due to a change in displacement u including the change in stress. Resolving (42) leads to i @Pa @P du = −(KT )−1 − (43) dY @Y @Y which is a matrix with dimension (ndof · nng ) × (ndm · ndn ). The indices ndof , nng , ndm and ndn denote the number of nodal degrees of freedom, number of nodes, number of spatial dimension Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 ISOPARAMETRIC FINITE ELEMENT CONCEPT 197 and number of design nodes, respectively. Note again that we assume conservative external forces in this paper, thus Pa 6= Pa (u). Tsay and Arora [4], for instance, comment on non-conservative external forces, Smith et al. [29] discuss design-dependent boundary conditions. In the terminology of Mahnken and Stein [5] the displacement u is the primal unknown of the direct problem (23) and considered as a dependent variable in the optimization process while the design vector Y serves as the key variable of the shape optimization problem. It is thus consistent to formulate the change in displacement in terms of changes of design via (43). Finally, we compute the sensitivities of the vectors of internal and external forces, Pi and Pa . The derivative of the latter follows directly from its denition (25) with (38) while the expression @Pi =@Y requires further argumentation. From the notation adopted in (42) it is clear that all inuences of the (dependent) displacement eld on the stress eld are contained in the consistent tangent matrix KT . We thus need to care for the dependency of the element volume and of the discrete strain displacement operator Be (see equation (24)) on the design vector Y or—expressed in terms of element nodal co-ordinates—X[m] , m = 1; : : : ; nnle . Interpreting the operator Be as a symmetric spatial gradient operator according to (3) and applying the denitions (19), (21) and (22), using the chain rule, we obtain @Pii[l] @Xk[m] Z = @N[l] @ det(Jeiso ) @N[l] @N[m] det(Jeiso ) + ij − dW @xk @xj @xj @Xk [m] (44) where i and k are co-ordinate indices, i; k = 1; : : : ; ndm , and l and m are indices associated to element nodes, l; m = 1; : : : ; nnle . Equation (44) thus denes the sensitivity of the local internal force vector of node l with respect to the position of node m in Euclidian space R ndm . The global sensitivity @Pi =@Y follows again with (35) and appropriate assembly of all element contributions. 3.4. Approximation procedure for iterative design updates The optimization algorithm used herein is a line search method using a third-order-accurate approximation of the objective function in the parameter space R N . Its description is split into two parts: The approximation itself and an accompanying scheme for ecient treatment of so-called linearly explicit constraints (32). These constraints simply constitute a feasible domain in the parameter space with a piecewise linear boundary. 3.4.1. Conguration update. The cubic objective approximation allows for a line search method that starts from a given initial conguration and improves the objective function with every design iteration. A trial conguration is evaluated for each design step in order to determine the optimal step length in direction of the steepest descent, expressed in terms of the negative gradient of the objective function, cf. (37). Box 1 displays the dierent steps inherent in the design iteration process. A conguration update is executed after the nal load level has been attained. According to Tsay and Arora [4] it is not necessary to account for intermediate load levels less than maximum load for stable hyperelastic problems. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 198 E. KIRCHNER AND R. IMMEL Box 1. Steps of the cubic objective approximation method. 1. Set initial conguration Y(0) and initial step length (0) 2. Compute structural response and simultaneously the objective function (0) and d = ∇(0) , set iteration counter j = 0 its gradient dY (0) 3. Set trial conguration with a trial step length ( j) , cf. (46) Y(∗j+1) = Y( j) − ∇( j) ( j) |∇( j) | (54) 4. Compute structural response, objective function ∗( j+1) and its gradient ∇(∗j+1) 5. Check termination criterion, e.g. stop if |∇(∗j+1) |6TOL holds. Minimum conguration is Y opt = Y(∗j+1) ˜ 6. Determine cubic approximation polynomial () with ∈ [0; 1], see Box 2 7. Set new conguration Y( j+1) that improves the objective 8. Increase iteration counter j = j + 1, perform next trial step (3) The values of the objective function ( j) and ∗( j+1) for the starting and the trial conguration as well as the gradients for the respective states are utilized to determine the four parameters that ˜ constitute the cubic approximation (), see Box 2. As discussed by Freed and Iskovitz [30] such a cubic approximation is accurate to the fourth order for a given step length ( j) . (See Shampine [31] for an error analysis for such an interpolation scheme.) Box 2. Determination of the cubic approximation polynomial. ˜ Approximation polynomial: () = ˜ 0 + ˜ 1 + ˜ 2 2 + ˜ 3 3 ˜ ˜ 2 = −30 + 31 − 200 − 01 Coecients: 0 = 0 0 ˜ 1 = 0 ˜ 3 = 20 − 21 + 00 + 01 Abbreviations: 0 = (Y( j) ) 00 = −( j) |∇j | ∇( j) · ∇∗( j+1) 1 = (Y(∗j+1) ) 01 = −( j) |∇( j) | The trial conguration has to be located in the search direction in a reasonable distance from the actual (starting) conguration. The distance in parameter space R N is expressed in terms of a step length parameter ( j) that eects the eciency of the algorithm. If the trial step is too large, the approximation polynomial becomes inaccurate and the quality of the conguration that improves the objective function suers. If the step length is too small the optimized conguration ( j) (45) Y( j+1) = Y( j) − opt ∇( j) k∇( j) k will be located at opt = 1. In the latter case, the trial step Y(∗j+1) is used as a new improved conguration. Thus, the algorithm tends to become inecient because too many small trial steps need to be executed for a satisfying conguration of design variables. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 ISOPARAMETRIC FINITE ELEMENT CONCEPT 199 ˜ = 0 with extr ∈ [0; 1] occur. If the step length ( j) is suciently large, then roots extr of d =d ˜ If such a root corresponds to a (local) minimum of the cubic approximation polynomial (), i.e. 2 ˜ 2 opt d =d ¿0 at = extr , then we set = extr . So, if the current step was larger than necessary to come to an optimum point of the approximation ˜ in direction of −∇ in design parameter space, a correction is performed. As mentioned before, this correction is based on a third-order-accurate approximation of the objective funtion. This advantage of providing an estimate for a correction in case of an unsuccessful step is not oered by other simple gradient-based optimization methods. It is desirable that the trial step length is just below the step length that results in the optimal conguration in the respective direction. By that the trial step is accepted every time. An adaptive trial step length control is applied to overcome the shortcomings of a constant trial step length this situation and to increase the eciency of the predictor step. The most simple choice for a step length adaptivity is opt ( j) if opt ¡1 (46) ( j+1) = if opt = 1 d ( j) In (46), the parameter d is a second parameter used in the cubic approximation algorithm. It determines the increase of step length after an accepted predictor step. The performance of this simple line search algorithm is demonstrated in Section 4. 3.4.2. Treatment of explicit linear constraints. As described before, the major benet of the present formulation of the fully stressed design problem is the absence of implicit constraints on the design variables. All constraints considered here, cf. (32), are linear and decoupled in the design variables. They simply dene upper and lower bounds for the design variables. This allows for a very simple treatment of the constraints for each design variable individually. The idea is the one behind the methods of feasible directions described by Zoutendijk [26] or Haftka and Gurdal [27]. Any search direction in the design space R N of the cubic objective approximation algorithm that violates a constraint is projected into a feasible direction that satises all constraints (32). According to Zoutendijk, an optimization problem with (linear) constraints on Y ∈ R N can be written in the form min[(Y) | AY6b; EY = e] N N (47) N with matrices A; E ∈ R × R and vectors b; e ∈ R . In general, the matrices A and E represent any kind of dierential algebraic operator that can be applied to the design parameters to obtain explicit non-linear constraints, which may couple the design variables. If these matrices are constant, A 6= A(Y), we have linear constraints. If only upper and lower bounds for the design variables are to be considered A and E are diagonal matrices. The vector b represents exactly these bounds in the case of linear decoupled constraints. Equality constraints EY = e can be reformulated as inequality constraints in (47) and thus they do not need to be treated separately. The set of inequality constraints must be divided into the active or binding subset A1 Y = b1 and the inactive subset A2 Y¡b2 where AT is decomposed into (A1T ; A2T ) and bT = (b1 T ; b2 T ). In the case of active constraints, the corresponding design variables in the actual design conguration Y are located on the boundary of the feasible domain, cf. Figure 3(a). Constraints that are inactive for the actual conguration might be active for the trial conguration which has to be dealt with separately in the sequel, cf. Figure 3(b). Figure 3(a) depicts the rst case in an example with two design variables, Y1 and Y2 . The constraint Y2 6Y2;max is active and the negative gradient −∇( j) would cause a violation of this Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 200 E. KIRCHNER AND R. IMMEL Figure 3. Treatment of linear constraints. constraint for the trial step. The projection simply sets the respective gradient entry—@=@Y2 for the example used in Figure 3(a)—to zero. For the computation of the trial and the improving ˜ ( j) obtained this way is used rather than the full actual conguration, the feasible gradient ∇ gradient ∇( j) . Figure 3(b) illustrates a case where the trial conguration Y(∗n+1) would violate an inactive constraint for design state Y( j) . It makes no dierence whether the gradient pointing towards the trial conguration has been modied yet to meet feasibility or not. In this case, the trial step length needs to be scaled down by a factor ∈ [0; 1] for the trial conguration to stay within the bounds. So we have to nd the value min = mini=1; : : : ;N [i ] of all design parameters Yi , i = 1; : : : ; N , with @( j) ( j) 6Yimax Yimin 6 Yi ( j) − i @Yi k∇( j) k {z } | ∗ Yi( j+1) (48) ∗ to produce a feasible trial conguration Ỹ( j+1) . So far, we have summarized all steps necessary for a shape optimization procedure seen from an engineer’s point of view. These comprise denition of objective functions, shape description, sensitivity analysis, and treatment of explicit linear constraints during design iteration. However, even though the discussion herein is restricted to non-linear hyperelasticity, it captures the key problems of shape optimization in general non-linear solid mechanics. Kirchner [12; 3] presents further examples using anisotropic viscoplastic materials at small strains where an objective of self-limiting type as in (31) is succesfully implemented. In principle, the algorithmic problems connected to shape optimization are all covered, but the time-dependent stress computation and the consistent sensitivity analysis taking the whole history dependency of the stress eld into account is much more complex. 4. NUMERICAL EXAMPLES Three examples demonstrate the applicability, capability, and performance of the presented objective and of the shape optimization approach. The rst one uses a simple one-dimensional truss where the cross-section is optimized. A comparison with an exact solution shows illustratively the Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 ISOPARAMETRIC FINITE ELEMENT CONCEPT 201 Figure 4. Geometry of the truss bar example. benets of the discussed modication of the objective (29). Furthermore, two three-dimensional applications to elastic materials at nite strains are computed in a nite element context. The method discussed is implemented in the nite element programme FEAP by Zienkewicz and Taylor [23; 24] and employs a fully consistent linearization (17) of the weak form of equlibrium (15). 4.1. Assessment of the modied objective function This brief example illustrates the feasibility and capabilities of the modied version of the classical approach of fully stressed designs from an engineering point of view. We consider a vertical truss bar, clamped at its top end, that is loaded by its own weight and an additional force F, cf. Figure 4. The stress distribution in the bar is given by Z L 1 F+ gA(x) dx (49) (x) = A(x) 0 The objective is to achieve a homogenous stress state (x) = y ∀x ∈ [0; L] by determining the optimal cross-section distribution A(x). The exact solution of this problem is given as gx F exp (50) A(x) = y y which yields a constant stress (x) ≡ y ∀x ∈ [0; L]. Introducing three design variables {A0 ; A1 ; A2 } we approximate the cross-sectional area by a second-order polynomial A(x) = A0 + A1 x + A2 x2 (51) This simple problem was investigated with the set of material and geometry parameters shown in Figure 4. A xed point search was used to compute the congurations given in Table I for three dierent optimization problems, for a more comprehensive discussion including convergence behaviour of dierent solution strategies we refer to Immel [34]. Table I shows a list of computed shape parameters for the truss bar (51) within three dierent optimizations. The resulting cross-sectional area is depicted in Figure 5 and the corresponding stresses are plotted in Figure 6. Design 1 in Table I represents the classical formulation class (27) including nonlinear implicit constraints on the design variables. The resulting maximum stress in the bar is y = 250 MPa. This value reached two times at positions x = 0:35L and x = L. The stress at the lower end (x = 0) is Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 202 E. KIRCHNER AND R. IMMEL Table I. Truss bar results for dierent objective functions, y ; t in MPa, volume in mm3 : Design 1 2 3 Objective A0 (mm) A1 (mm=m) A3 (mm=m2 ) max (N=mm2 ) class (y = 250) 0.434637 −1.464481 mod (t = 250) 0.048558 −0.020335 mod (t = 200) 0.076476 −0.206634 3.548807 0.553001 1.514208 249.9997 285.4438 248.9262 Volume (%) 265 165 100 66 769 25 143 292 54 Figure 5. Cross-sectional areas for dierent objective functions. Figure 6. Stress distributions for dierent objective functions. only about 10 per cent of y which is extremely small. Figure 5 explains the reason for this. The cross-sectional area is completely over-dimensioned at the bottom compared to the analytical solution. The volume of design 1 is used as a reference for the modied objective. Design 2 in Table I presents the results of the modied objective mod (29) which correspondingly uses t = 250 MPa as target stress. The stress exceeds the former limit y for x=L ∈ [0:07; 0:32] Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 ISOPARAMETRIC FINITE ELEMENT CONCEPT 203 Figure 7. Three-dimensional elastic beam submitted to end forces F. and x=L ∈ [0:70; 1:00] as indicated in Figure 6. This is an important feature of the modied approach. Figure 5 indicates that the cross-sectional area is always smaller than that of the analytical solution, which explains both the high stress values and the small volume which is only 25 per cent of the reference volume. This is a rather desirable fact that greatly supports the alternative approach. If partial plastic deformation accompanied with elastoplastic hardening were permitted this would be a very good design. If the target stress t is an absolute upper limit, the alternative objective function can still be applied with a slight modication. This leads to the third design which is obtained with a pragmatic engineering approach. The maximum stress in design 2 exceeds the target value by about 14 per cent This suggests to simply decrease the desired stress value t by a similar amount and rerun the optimization process. In this design the target stress was reduced by 20 per cent to a value of t = 200 MPa. The result is very convincing since the maximum stress stays just below the original limit of y = 250 MPa and the volume is only 54 per cent of the reference volume. Figure 5 illustrates that the cross-sectional area of this design stays closest to the analytical solution. This example justies the self-limiting modication of the objective for the following reasons: All implicit nonlinear constraints collapse to explicit linear constraints in the design variables which simplies their treatment signicantly. A simple adaption of the objective parameters can correct deviations from the original objective. After a ‘preprocessing’ type optimization we are able to detect a solution of type (51) for this problem that is by far better than the reference solution obtained with the classical approach for fully stressed design based on the objective function (27) but meets the same requirements. 4.2. Shape optimization of a slender beam To demonstrate the abilities of the modied approach of fully stressed design and of the cubic approximation algorithm we consider a beam with varying cross-sectional parameters, cf. Figure 7. Empty circles in that gure denote design nodes with variable position in the X –Y plane, lled circles are xed. The material parameters for the Neo-Hookean material (11) are chosen as = 1 000 MPa and = 2 000 MPa (52) We start with cross-sections 1 and 2 being bi-unit squares. For the optimization algorithm sketched in Box 1 we use the initial values displayed in Figure 8 for ( 0) and an increase parameter d = 1:2 in (46). The desired stress is d = 25 MPa. The algorithm presented here preserves the symmetry of the problem as expected which is indicated in Figure 7. Neither unsymmetric deformations nor design modications were observed. To justify the assumptions made when modifying the objective function in Section 3.1, we monitored the minimum and maximum equivalent stresses of initial and optimized congurations for three dierent element numbers nel , see Table II. The maximum stresses are always obtained in Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 204 E. KIRCHNER AND R. IMMEL Figure 8. Objective function versus iteration number for three dierent meshes and dierent initial step lengths (0) . Table II. Comparison of initial and optimized equivalent stress extrema. No. 1 2 3 4 5 nel (0) 16 128 128 1024 1024 0.50 0.05 2.00 0.05 2.00 initial eq; min 8:72 6:57 6:57 1:76 1:76 MPa MPa MPa MPa MPa initial eq; max 10:81 21:26 21:26 38:99 38:99 optimized eq; min MPa MPa MPa MPa MPa 14:83 10:49 11:30 1:17 0:99 MPa MPa MPa MPa MPa optimized eq; max 26:02 33:64 35:21 59:71 55:39 MPa MPa MPa MPa MPa the rst element layer at the clamping, Z = 0. With increasing element numbers some quadrature points come arbitrarily close to this stress singularity. Constraints applied to these stress maxima would prevent an optimization algorithm from approaching the typical double-T-shaped cross-section known from basic engineering mechanics. It is clear from Table II that a dependency of the optimized structural behaviour on the nite element model exists. In our computations we observed a moderate dependency of the optimization result on the chosen mesh renement, but for increasing element numbers the solution rapidly converges against the solution of computation number 5, cf. Table II, with 1024 elements. However, we expect the classical objective (27) with constraints (28) to be more sensitive to changes in the computational mesh than our approach because the constraints (28) need to be considered for each quadrature point separately. The better resolution of the peak stresses near the clamping leads to increasing maximum stresses and additional constraints become active. This shortcoming drops out in our method. A quadrature point which exceeds the desired equivalent stress eq d penalizes the objective (27) and thus diminishes the quality but does not impose additional constraints. If we consider the decrease of the objective function in Figure 8 for ve characteristic computations, the dependency on the choice of the initial step length ( 0) in design parameter space becomes obvious. For each curve we performed up to 15 structural computations, within which Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 ISOPARAMETRIC FINITE ELEMENT CONCEPT 205 Examples 1, 3 and 5 obviously attained a stationary point. For this solution, we found 4 out of 12 design degrees of freedom located at the boundary of the feasible domain in design space, the remaining 8 slightly oscillated around a xed value. The algorithm approaches this extremal conguration in the interior of the feasible domain arbitrarily close. The cubic approximation polynomial (Box 2) increases the eciency of the standard steepest opt descend method (i.e. ≡ 1) as expected. In Figure 8, we added symbols where a reduction of the step length was indicated by the proposed cubic approximation scheme. Unfortunately, even a third-order-accurate approximation of the objective does not guarantee an increase of design quality in each step. This is due to the mild tendency of a third-order polynomial approximation to predict non-physical changes of the approximated objective function, a well-known result in approximation theory. 4.3. Fully stressed design of a notched bar The third example has been taken from Barthold [8] where the problem was solved using the objective ‘minimum weight’ subjected to stress constraints (28). The material is of Ogden type governed by the stored energy function (12) and assumed as (almost) incompressible. To avoid locking phenomena appearing for incompressible materials modelled by displacement-based elements we use a nite bulk modulus instead of demanding det(M) = J ≡ 1 during our simulations. Precisely, we use the following material and optimization parameters = 10:00 MPa; 1 = 0:63 MPa; 2 = −0:01 MPa; 3 = 0:0012 MPa d = 9 MPa; 1 = 1:3; 2 = − 2:00; 3 = 5:00 (53) We consider a notched plate as shown in Figure 9. The thickness is t = 1 mm and the problem is assumed symmetric with respect to the X - and Y -axis, respectively. The specimen is subjected to a vertical elongation v0 = 36 mm leading to an overall deformation of 200 per cent. The upper and lower boundary surfaces perpendicular to the Z-axis, see Figure 9, are assumed stress-free. We started the optimization procedure with the design model shown in Figure 9 using four design elements with a total of 30 design nodes, one layer of 15 nodes at z = 0 and a second layer at z = t = 1 mm. In each layer 9 design nodes are xed and 6 considered as real optimization variables. Again, we investigated four dierent initial step lengths. The convergence of the algorithm is shown in Figure 10. Bullets denote corrections of respective trial congurations indicated by the cubic approximation. The approximation algorithm performs well. Several large steps were succesfully detected and lead to an acceleration of the step length search procedure. As in the previous example, small initial step lengths (0) need more design iterations which of course could be compensated by an improved choice of the increase parameter d . We used d = 1:2 for all runs displayed here. The distribution of the equivalent Kirchho stress for the initial and optimized conguration (cf. Figure 11), allows for two conclusions. First, for the optimized design, the material is more uniformly loaded as expected. Second, the peak stresses occuring in the notch ground are reduced despite the slendering of the notch ground region. The desired stress d is not exceeded for the optimized conguration although the modied objective function (29) allows for local stress peaks higher than d . A comparison with Barthold’s [8] results (cf. Figure 12) shows a qualitatively good agreement. The dierences between the two optimized congurations are mainly caused by the restriction to contours with vertical tangents on the right-hand side corners posed on Bartholds computations. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 206 E. KIRCHNER AND R. IMMEL Figure 9. Initial conguration and design model for the notched tension strip with four design elements (thick lines) and 30 design nodes arranged in two layers at z = 0 and 1 mm. Figure 10. Development of objective during iteration for dierent initial step size (0) : However, despite the simplicity of the proposed method, the behaviour of the optimized structure is remarkable from an engineer’s point of view. Not only a weight reduction but also a decrease of the maximum stress occurred for this example although both properties are often considered as contradictory. 5. CONCLUSIONS We propose a simple shape optimization approach that is fully on closed-form expressions for all sensitivities involved. The transition from the underlying design model to a computational model is simple and does not require large numbers of design parameters. Implementation of the method Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 ISOPARAMETRIC FINITE ELEMENT CONCEPT 207 Figure 11. Von Mises-equivalent stress for initial and optimal conguration. Figure 12. Comparison of the present result (a) with Barthold’s optimized conguration (b). is comparatively easy and is restricted to adding some new macro commands to an existing FEM code. The results obtained with the optimization approach are very much usable as shown in the examples presented, but we cannot guarantee that a given maximum stress is not exceeded. Avoiding implicit constraints on dependent quantities leads to an explicitly constrained optimization problem. This problem is eciently solved by a rst-order method combined with a feasible direction algorithm by taking into account the bounds imposed on the design parameters. We do believe that the presented method is not only simple to implement but also capable of solving real-world engineering problems. The rate of convergence could be increased by application Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 208 E. KIRCHNER AND R. IMMEL of more sophisticated solution procedures for the optimization problem. The price for this speedup, however, would be an increased exchange of information between the nite element algorithm and the optimization routines leading to a higher eort in coding. The main goal of this paper was to propose a method for design optimization of highly loaded structures which does not require the introduction of implicitly dened constraints formulated in terms of dependent quantities of the optimization problem such as stresses or displacements. The denition of a design independent target state for the objective function is preferable for future applications involving inelastic materials where well dened yield limits often do not exist. Furthermore, if one considers a certain mechanical behaviour under loading as the reference or target state, it is physically justied to permit variations above and below this state as demonstrated in this paper, cf. Figure 1. ACKNOWLEDGEMENTS The authors are grateful to Dr. F.-J. Barthold, Hannover, for providing his optimization results in the fully stressed design of a notched plate. REFERENCES 1. Wei X, Chandra A, Leu L-J, Mukherjee S. Shape optimization in elasticity and elasto-viscoplasticity by the boundary element method. International Journal of Solids and Structures 1994; 31:533–550. 2. Wu CC, Arora JS. Design sensitivity analysis and optimization of nonlinear structural response using incremental procedure. AIAA Journal 1986; 25:1118–1125. 3. Kirchner, E. Modelling single crystals: time integration, tangent operators, sensitivity analysis and shape optimization. International Journal of Plasticity 1999; Submitted for publication. 4. Tsay JJ, Arora JS. Nonlinear structural design sensitivity analysis for path dependent problems. Part 1: General theory. Computer Methods in Applied Mechanics Engineering 1990; 81:183–208. 5. Mahnken R, Stein E. Identication of parameters for visco-plastic models via nite-element methods and gradient methods. Modelling Simulations Materials Science Engineering 1994; 2:597–616. 6. Ramm E, Maute K. Transition from shape to topology optimization. In Second ECCOMAS Conference on Numerical Methods in Engineering, Paris, 9-13 September 1996, Desideri J-A. et al. (eds.). Wiley: New York, 1996. 7. Bletzinger K-U, Maute K, Reitinger R, Ramm E. Layout of linear and nonlinear structures by shape and topology optimization. In IUTAM Symposium on Optimization of Mechanical Systems, Bestle D, Schielen W (eds). Kluwer Academic Press: Dordrecht, 1996; 49–56. 8. Barthold F-J. Theorie und Numerik zur Berechnung und Optimierung von Strukturen aus isotropen, hyperelastischen Materialien. Forschungs- und Seminarberichte F 93=3, Institut fur Baumechanik und Numerische Mechanik der Universitat Hannover, 1993. 9. Falk A, Barthold FJ, Stein E. Hierarchical modelling in shape optimization. In Proceedings of 1st World Congress of Structural and Multidisciplinary Optimization, May 28–June 2, Goslar, Olho N, Rozvany GIN (eds). Elsevier Science: Oxford, 1995; 371–376. 10. Falk A. Adaptive Verfahren fur die Formoptimierung von Flachentragwerken unter Berusichtigung der CAD-FEMKopplung. Forschungs- und Seminarberichte F 95=1, Institut fur Baumechanik und Numerische Mechanik der Universitat Hannover, 1995. 11. Bletzinger, K-U. Formoptimierung von Flachentragwerken. Bericht No. 11, Institut fur Baustatik der Universitat Stuttgart, 1990. 12. Kirchner E. Ein Beitrag zur Formoptimierung in der nichtlinearen Strukturmechanik. Ph.D. thesis, Fachgebiet Maschinenelemente und Maschinenakustik der Technischen Universitat Darmstadt, 1999. 13. Schittkowski K. The nonlinear programming method of Wilson, Han, and Powell with an augmented lagrangian type line search function. Numerische Mathematik 1981; 38:83–114. 14. Bertsekas, D. Projected newton methods for optimization problems with simple constraints. SIAM Journal on Control and Optimization 1982; 20(2):221–246. 15. Reitinger R, Bletzinger K-U, Ramm E. Shape optimization of buckling sensitive structures. Computing Systems in Engineering 1994; 5:65–75. 16. Svanberg K. The method of moving asymptodes—A new method for structural optimization. International Journal for Numerical Methods in Engineering 1987; 24:359–373. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209 ISOPARAMETRIC FINITE ELEMENT CONCEPT 209 17. Bletzinger K-U, Extended method of moving asymptotes based on second-order information. Structural Optimization 1993; 5:175–183. 18. Wriggers P. Konsistente Linearisierung in der Kontinuumsmechanik und ihre Anwendung auf die Finite-ElementMethode. Forschungs- und Seminarberichte F 88=4, Institut fur Baumechanik und Numerische Mechanik der Universitat Hannover, 1988. 19. Marsden E, Hughes TJR. Mathematical Foundations of Elasticity. Prentice-Hall: Englewood Clis, 1983. 20. Ogden RW. Large deformation isotropic elasticity: on the correlation of theory and experiment for compressible rubberlike solids. Proceedings of the Royal Society of London, Series A 1972; 326:567–583. 21. Reese S. Theorie und Numerik des Stabilitatsverhaltens hyperelastischer Festkorper. Ph.D. thesis, Fachbereich Mechanik der TH Darmstadt, 1995. 22. Lubliner J. Plasticity Theory. Macmillan Publishing Company: New York, 1990. 23. Wriggers P. Nichtlineare Finite-Element-Methoden. Institut fur Mechanik IV, Technische Hochschule: Darmstadt, 1995. 24. Zienkiewicz OC, Taylor RL. The Finite Element Method (4th edn), Vol. 1. McGraw-Hill Book Company: London, 1989. 25. Zienkiewicz OC, Taylor RL. The Finite Element Method (4th edn), Vol. 2. McGraw-Hill Book Company: London, 1989. 26. Zoutendijk G. Method of Feasible Directions. Elsevier: Amsterdam, 1960. 27. Haftka RT, Gurdal Z. Elements of Structural Optimization (3rd edn). Kluwer Academic Publishers: Dordrecht, 1992. 28. Gurtin ME. An Introduction to Continuum Mechanics (1st edn). Academic Press: New York, 1982. 29. Smith DE, Tortorelli DA, Tucker CL. Optimal design for polymer extrusion. Part 1: Sensitivity analysis for nonlinear steady state systems. Computer Methods in Applied Mechanics and Engineering 1998; 167:283–302. 30. Freed AD, Iskovitz IS. Development and applications of a Rosenbrock integrator. NASA Memorandum 4709, NASA Aeronautics and Space Administration, 1996. 31. Shampine LF. Lipschitz constants and robust ODE codes. In Computational Methods in Nonlinear Mechanics, Oden JT (ed.). North-Holland Publishsing Company: New York, 1980; 427– 449. 32. Barthold F-J, Stein E. A continuum mechanical-based formulation of the variational sensitivity analysis in structural optimization. Part I: Analysis. Structural Optimization 1996; 11:29– 42. 33. Barthold F-J, Stein E. Implementing variational design sensitivity analysis. In Proceedings WCSMO-2 “Second World Congress of Structural and Multidisciplinary Optimization, Mroz Z, Gutkowsky W (eds). Zakopane, Poland, 1997; 181–186. 34. Immel R. Entwicklung und Implementierung eines Algorithmus zur Formndung in der nichtlinearen Strukturmechanik. Diplomarbeit, Fachgebiet Maschinenelemente und Maschinenakustik der Technischen Universitat Darmstadt, 1998. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 2000; 48:185–209

1/--страниц