close

Вход

Забыли?

вход по аккаунту

?

295

код для вставкиСкачать
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
Документ
Категория
Без категории
Просмотров
2
Размер файла
303 Кб
Теги
295
1/--страниц
Пожаловаться на содержимое документа