JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

44, 2, pp. 323-349, Warsaw 2006

APPLICATION OF FINITE VARIATIONS TO TOPOLOGY

AND SHAPE OPTIMIZATION OF 2D STRUCTURES

Dariusz Bojczuk

Faculty of Management and Computer Modelling, Kielce University of Technology

e-mail: mecdb@eden.tu.kielce.pl

Wojciech Szteleblak

PhD student; e–mail: szteleblak@go2.pl

The method of simultaneous topology and shape optimization of 2D
structuresbyfinite topologymodification ispresented in thepaper.Both,
structures in a plane state of stress and bending Kirchhoff’s plates are
analyzed here. Conditions for the introduction of finite topology mo-
dification based on the topological derivative are specified. When the
respective condition is satisfied, finite holes and finite variations of exi-
sting boundaries are introduced into the structure.Next, standard shape
optimization of newholes andvariableboundaries is performed.Twoba-
sic types ofmodification are considered here, namely the introduction of
holes of a prescribed size and shape and the introduction of holes of an
unknown size and shape together with the introduction of finite changes
of other boundaries. A heuristic algorithm for optimal design of topo-
logy and shape is proposed in the paper. Illustrative examples confirm
applicability of the proposed approach.

Key words: 2D structures, optimal topology and shape, topological
derivative, finite modification, structure evolution

1. Introduction

Theoptimal structural design is usually concernedwith the specification of
dimensional parameters, shape (or configuration) parameters and topology of
a structure. For 2D structures, the coupled topology and shape optimization
problem can be stated as optimal determination of external boundaries and
holes within a given domain. Now, at first we can introduce infinitesimally



324 D. Bojczuk, W. Szteleblak

small voids (so called ”bubbles”) into the structure of the optimal shape. The
position of the void, in the case ofminimization of the strain energy functional,
is determinedby applying the bubblemethod (Eschenauer et al., 1994), and in
the case of an arbitrary objective functional, using the topological derivative
approach. For plane elasticity, this problem was studied by Sokołowski and
Żochowski (1999), Burczyński (2002), Burczyński and Kokot (2005).
In thepresent paper, the approachby topological derivative is extended for

finite topology modifications. Now, instead inserting of infinitesimally small
voids, the modification corresponds to introduction of holes of finite dimen-
sions and introduction of finite variations of existing boundaries (cf.Mróz and
Bojczuk, 2003). Next, the boundaries of newholes are describedby some addi-
tional shape parameters, and standard shape optimization of variable external
boundaries and holes is performed.
When the approach by finite topology modification is used, usually the

number of iterations during the optimization process is considerably redu-
ced. Here, this approach is applied to plane elasticity problems and for plate
problems. In the case of plates, former methods are usually based on the
optimization of thickness distribution (cf. Marczewska et al., 2003), or on de-
termination of optimal distribution of twomaterials: weaker and stronger (cf.
Czarnecki et al., 2004). So, the approachpresentedhere differs essentially from
them and develops a newmethodology of optimal design.

2. Topological derivative for 2D structures

2.1. Topological derivative for plane elasticity problems

The concept of the topological derivative for strain energywas analyzed by
Eschenauer et al. (1994). Next, it was generalized for arbitrary displacement
and stress functionals by Sokołowski and Żochowski (1999).
Let Ω ⊂ R2 be a domain occupied by an elastic body with a boundary

Γ = Γu ∪ΓT . The body is subjected to the surface traction T
0 on ΓT and

volume forces p0 in Ω. On the boundary Γu, displacements u0 are given. The
topological derivative is defined as follows

=DG(x)= lim
ρ→0

G(Ω−Bρ(x))−G(Ω)
πρ2

x∈Ω (2.1)

where G is the analyzed functional and Bρ(x) is a circularhole of the radius ρ,
so Bρ(x)= {y∈R2 : |y−x| ¬ ρ}. Moreover, the superscript D denotes the
topological derivative for plane elasticity problems (disks).



Application of finite variations to topology... 325

The topological derivative provides information about infinitesimally small
variation of the functional G if a small hole is inserted into the structure at
the point x. It means that the topological derivative can be used for topology
and shape optimization of structures. To calculate the topological derivative,
an adjoint method is used. In this method, a new structure of the same shape
and material constants as the primary structure but with different boundary
conditions and loads, is introduced.
Consider for a disk of a unit thickness a functional of strains and displa-

cements in the form

G=
∫

Ω

F(ε) dΩ+
∫

Ω

f(u) dΩ (2.2)

where F is a function of strains ε = [ε11,ε22,γ12]> (here γ12 = 2ε12) and
f denotes a function of displacements u. Its variationwith respect to insertion
of a hole of an infinitesimally small radius ρ at the point x, can be written
as follows

δG=
∫

Ω

∂F

∂ε
δε dΩ+

∫

Ω

∂f

∂u
δu dΩ+

2∑

k=1

∫

Γρ

(F +f)nkδϕk dΓρ (2.3)

where n= [n1,n2]> is the vector normal to the boundary Γρ of the hole and
δϕ= [δϕ1,δϕ2]> denotes the function of shape transformation of the hole.
Consider now an adjoint structure of the same shape as the primary struc-

ture, but subjected to initial stresses σai and volume forces pa0, namely

σ
ai =
∂F

∂ε
p
a0 =

∂f

∂u
inΩ (2.4)

The boundary conditions are assumed as follows

T
a0 =0 on ΓT u

a0 =0 on Γu (2.5)

where Ta0 is the surface traction and ua0 is the vector of displacements. The
field of initial stresses induces a field of global stresses σa = [σa11,σ

a
22,σ

a
12]
> in

the form
σ
a =σai+σar (2.6)

where σar is the field of the elastic stresses. Now, taking into account (2.4),
(2.5) and the relationship

∑2
k=1nkδϕk =−δρ, equation (2.3) becomes

δG=
∫

Ω

σ
aiδε dΩ+

∫

Ω

p
a0δu dΩ−

∫

Γρ

(F +f) dΓρδρ (2.7)



326 D. Bojczuk, W. Szteleblak

So, in view of the virtual work equation
∫

Ω

σ
arδε dΩ=

∫

Ω

p
a0δu dΩ (2.8)

the complementary virtual work equation
∫

Ω

σ
aδε dΩ =

∫

Ω

ε
aδσ dΩ =

∫

Γρ

σε
a dΓρδρ−

∫

Γρ

p
0
u
a dΓρδρ (2.9)

and equation (2.6), the first variation of the functional G takes the form

δG=
∫

Γρ

(σεa−F −f−p0ua) dΓρδρ (2.10)

Now, adapting this sensitivity formula to the case of a disk with an arbitra-
ry thickness h, and assuming that the function f and external load p0 are
referred directly to this thickness, equation (2.10) can be rewritten as follows

δG=
2π∫

0

[(σεa−F)h−f−p0ua]ρ dθδρ (2.11)

As for the considered case δG= (∂G/∂ρ)δρ, it is important to notice that if
ρ→ 0, then ∂G/∂ρ → 0. Therefore, let us compute the second derivative at
the point x for ρ=0. Now, we have

∂2G

∂ρ2

∣∣∣
ρ=0
= lim
ρ→0

2π∫

0

{ ∂
∂ρ
[(σεa−F)h−f−p0ua]+

(2.12)

−
∂

∂n
[(σεa−F)h−f−p0ua]+

1
ρ
[(σεa−F)h−f−p0ua]

}
ρ dθ

On the boundary of the hole, in the polar coordinates r, θ with the origin in
its center, we have: σrr =0, σrθ =0, σ

a
rr =0, σ

a
rθ =0. So, only hoop stresses

σθθ, σ
a
θθ may attain non-zero values on this boundary. They can be expressed

as follows (cf. Sokołowski and Żochowski, 1999)

σθθ =(σ1+σ2)−2(σ1−σ2)cos2θ
(2.13)

σaθθ =(σ
a
1 +σ

a
2)−2(σ

a
1 −σ

a
2)cos2(θ−α)



Application of finite variations to topology... 327

where σ1, σ2, σa1, σ
a
2 are the principal stresses in the considered point, re-

spectively for the primary and adjoint structure, and α is the angle between
corresponding principal directions. Taking into account that εaθθ =σ

a
θθ/E and

∂/∂n=−(∂/∂r)|r=ρ, relation (2.12) becomes

∂2G

∂ρ2

∣∣∣
ρ=0
= lim
ρ→0

2π∫

0

{ ∂
∂ρ

[( 1
E
σθθσ

a
θθ−F

)
h−f−p0ua

]
+

+
∂

∂r

[( 1
E
σθθσ

a
θθ−F

)
h−f−p0ua

]∣∣∣
r=ρ
+ (2.14)

+
1
ρ

[( 1
E
σθθσ

a
θθ−F

)
h−f−p0ua

]}
ρ dθ

where E is Young’s modulus.
In the structure with the hole of the radius ρ, for displacements expressed

in the polar coordinates r (r­ ρ) and θ, the following expansion in the point
ρ=0 holds (cf. Il’in, 1992; Sokołowski and Żochowski, 1999)

ur(r,θ)=ur0+
σ1+σ2
4Gr

(1−ν
1+ν

r2+ρ2
)
+

+
σ1−σ2
4Gr

(
r2+

4
1+ν

ρ2−
ρ4

r2

)
cos2θ+ . . .

(2.15)

uθ(r,θ)=uθ0−
σ1−σ2
4Gr

(
r2+2

1−ν
1+ν

ρ2+
ρ4

r2

)
sin2θ+ . . .

where ν is Poisson’s ratio, G = E/[2(1+ ν)] and ur0, uθ0 are displacement
components of the unmodified structure at the point x. Next, using geome-
trical relations

εrr =
∂ur
∂r

γrθ =
1
r

∂ur
∂θ
+
∂uθ
∂r
−
uθ
r

(2.16)

εθθ =
ur
r
+
1
r

∂uθ
∂θ

and Hooke’s law

σrr =
E

1−ν2
(εrr+νεθθ) τrθ =Gγrθ

(2.17)

σθθ =
E

1−ν2
(εθθ +νεrr)



328 D. Bojczuk, W. Szteleblak

in the polar coordinates, the stresses can be expressed in the form

σrr(r,θ)=
1
2
(σ1+σ2)

(
1−
ρ2

r2

)
+
1
2
(σ1−σ2)

(
1−4

ρ2

r2
+3
ρ4

r4

)
cos2θ+ . . .

σθθ(r,θ)=
1
2
(σ1+σ2)

(
1+
ρ2

r2

)
+
1
2
(σ1−σ2)

(
1+3

ρ4

r4

)
cos2θ+ . . . (2.18)

τrθ(r,θ)=−
1
2
(σ1−σ2)

(
1+2

ρ2

r2
−3
ρ4

r4

)
sin2θ+ . . .

Analogous relations can be also obtained for the adjoint structure. After sub-
stitution of expansions of the stresses σθθ, σ

a
θθ into (2.14), it can be easily

proven that the sum of two first terms on the right-hand side of this equation
cancels out. Now, (2.14) can be rewritten in the form

∂2G

∂ρ2

∣∣∣
ρ=0
=
2π∫

0

[( 1
E
σθθσ

a
θθ−F

)
h−f−p0ua

]
dθ (2.19)

Then, taking into account (2.13) and integrating with respect to θ, the topo-
logical derivative at the point x, which is expressed in the form

=DG(x)=
1
2π
∂2G

∂ρ2

∣∣∣
ρ=0

(2.20)

finally, can be written as follows

=DG(x)=h
[ 1
E
(σ1+σ2)(σ

a
1+σ

a
2)+
2
E
(σ1−σ2)(σ

a
1−σ

a
2)cos2α−F0

]
−f−p0ua

(2.21)
Here all quantities are calculated at the point x, and the expression
F0 = (2π)−1

∫2π
0 F dθ should be determined separately for each form of the

function F . In the case of the strain energy U, the adjoint structure is the
same as the primary one. Moreover, assuming that f =0 and p0 = 0 on the
boundary Γρ of a new hole, (2.21) becomes

=DG(x)=
h

2E
[(σ1+σ2)

2+2(σ1−σ2)
2] (2.22)

Now, let us consider the following functional of stresses and reactions

G=
∫

Ω

H(σ) dΩ+
∫

Γu

g(T) dΓu (2.23)



Application of finite variations to topology... 329

where H(σ) is the function of stresses σ= [σ11,σ22,σ12]> and g(T) denotes
the function of reaction forces acting on the boundary Γu. In this case, the
adjoint structure is of the form

ε
ai =
∂H

∂σ
pa0 =0 in Ω

Ta0 =0 on ΓT u
a0 =−

∂g

∂T
on Γu

(2.24)

where εai are the initial strains. The field of the initial strains induces a field
of global strains in the form εa = εai +εar, where εar is the field of elastic
strains. Now, in view of (2.24), the sensitivity of functional (2.23) with respect
to introduction of a hole of an infinitesimally small radius ρ at the point x,
can be expressed analogously to (2.7), namely

δG=
∫

Ω

ε
aiδσ dΩ−

∫

Γu

u
a0δT dΓu−

∫

Γρ

H dΓρδρ (2.25)

where Γρ denotes the boundary of this hole. Next, in view of the virtual work
equation ∫

Ω

ε
arδσ dΩ =

∫

Ω

σ
arδε dΩ =0 (2.26)

and the complementary virtual work equation
∫

Ω

ε
aδσ dΩ=

∫

Γρ

σε
a dΓρδρ+

∫

Γu

u
a0δT dΓu−

∫

Γρ

p
0
u
a dΓρδρ (2.27)

the first variation of functional (2.23) takes the form

δG=
2π∫

0

[(σεa−H)h−p0ua]ρ dθδρ (2.28)

It is easy to notice that for ρ = 0, the first derivative of the considered
functional is equal zero. So, for ρ=0, let us determine the second derivative
of this functional. Then, we have

∂2G

∂ρ2

∣∣∣
ρ=0
= lim
ρ→0

2π∫

0

{ ∂
∂ρ
[(σεa−H)h−p0ua]−

∂

∂n
[(σεa−H)h−p0ua]+

(2.29)

+
1
ρ
[(σεa−H)h−p0ua]

}
ρ dθ



330 D. Bojczuk, W. Szteleblak

Analogously as in the case of functional (2.2), the first two terms on the right
hand side of (2.29) cancel out and finally the topological derivative can be
written as follows

=DG(x)=h
[ 1
E
(σ1+σ2)(σ

a
1 +σ

a
2)+

2
E
(σ1−σ2)(σ

a
1 −σ

a
2)cos2α−H0

]
−p0ua

(2.30)
where

H0 =
1
2π

2π∫

0

H dθ

2.2. Topological derivative for bending plates

Consider a plate occupying, referred to themiddle surface, a domain Ω⊂
R2, with a boundary Γ . Theplate is subjected to the transverse load p0 in Ω,
whereas either the generalized traction T0 ordisplacements are specifiedon Γ .
Here, Kirchhoff’s theory of thin plates is used, and the bending plate can be
treated as a set of layers, each of which is in a plane state of stress.
Now, we confine our analysis to the case when the functional G is expres-

sed by (2.2). Here, taking into account that u = [0,0,w]> on the middle
surface, f(u) can be interpreted as a function of transverse displacements w.
Then, the topological derivative of the analyzed functional with respect to the
introduction of an infinitesimally small, circular hole at the point x=(x1,x2)
of themiddle surface, can be determined analogously as in Section 2.1. So, in
view of (2.21), we have

=PG(x)=

h/2∫

−h/2

[ 1
E
(σ1+σ2)(σ

a
1+σ

a
2)+
2
E
(σ1−σ2)(σ

a
1−σ

a
2)cos2α−F0

]
dx3−f−p

0
u
a

(2.31)
where h denotes the thickness of the plate and x3 is the coordinate normal
to the middle surface. The principal stresses may be expressed as a function
of x3 in the form

σi(x3)= 2
σ̂i
h
x3 σ

a
i (x3)= 2

σ̂ai
h
x3 i=1,2 (2.32)

where σ̂i, σ̂ai are the principal stresses on the upper surface of the plate,
respectively for the primary and adjoint structure.After substitution of (2.32)
into (2.31) and integration, we have



Application of finite variations to topology... 331

=PG(x)=
h

3E
[(σ̂1+ σ̂2)(σ̂

a
1 + σ̂

a
2)+2(σ̂1− σ̂2)(σ̂

a
1 − σ̂

a
2)cos2α]+

(2.33)

−

h/2∫

−h/2

F0 dx3−f−p
0
u
a

Taking into account that the principal moments in the primary and adjoint
structure take the form

Mi =

h/2∫

−h/2

σix3 dx3 M
a
i =

h/2∫

−h/2

σaix3 dx3 i=1,2 (2.34)

the topological derivative can be also presented as follows

=PG(x)=
12
Eh3
[(M1+M2)(M

a
1 +M

a
2)+2(M1−M2)(M

a
1 −M

a
2)cos2α]+

(2.35)

−

h/2∫

−h/2

F0 dx3−f−p
0
u
a

When G coincides with the strain energy and assuming that f = 0 and
p0 = 0 on the boundary Γρ of a new hole, the topological derivative reduces
to the form

=PG(x)=
h

6E
[(σ̂1+ σ̂2)

2+2(σ̂1− σ̂2)
2] (2.36)

It can be also expressed by the principal moments, namely

=PG(x)=
6
Eh3
[(M1+M2)

2+2(M1−M2)
2] (2.37)

Similar considerations can be performed aswell for an arbitrary functional
of stresses and reaction forces expressed by (2.23).

2.3. Topological derivative for cost functional

Next, let us consider the topological derivative of the global cost C. The
cost of the structure, for example, can be assumed as proportional to the
material volume. Then, it can be expressed as follows

C = c(V0−πρ
2h) (2.38)



332 D. Bojczuk, W. Szteleblak

where c is aunit cost, V0denotes the initial volumeof the consideredstructure,
and, as previously, h denotes its thickness and ρ is the radius of the inserted
small hole. Thus, the topological derivative of the cost functional with respect
to the introduction of this infinitesimally small hole, takes the form

=C(x)=
1
2π
∂2C

∂ρ2

∣∣∣
ρ=0
=−ch (2.39)

3. Optimal design method based on finite modification

In this section, conditions for the introduction of finite topologymodifica-
tion based on the topological derivative are derived. Also, sensitivity expres-
sions with respect to the shape modification are presented. These conditions
and formulas can be applied in order to formulate a heuristic algorithm for
simultaneous topology and shape optimization.

3.1. Conditions for introduction of finite topology modifications

We consider now a general optimization problem of the form

minG subject to C−C0 ¬ 0 (3.1)

where G is the objective functional (or function), C denotes the global cost
(or global volume) and C0 is the upper bound on the global cost. When G
expresses the strain energy, the problem corresponds to maximization of the
structure stiffness with the constraint set on the global cost.
At first, we take into account the optimal designwith respect to dimensio-

nal and shape parameters si, i=1,2, . . . ,n. Using Lagrangian

L=G+λ(C−C0) (3.2)

where λ (λ­ 0) is the Lagrangemultiplier, the optimality conditions take the
form

∂L

∂si
=
∂G

∂si
+λ
∂C

∂si
=0 i=1,2, . . . ,n

(3.3)

λ(C−C0)= 0

Therefore, the optimal values of design parameters si and Lagrange’s multi-
plier λ can be determined in the incremental process of optimization.



Application of finite variations to topology... 333

Now, the condition of introduction of an infinitesimally small circular hole,
using the concept of the topological derivative, takes the form

=(x)==G(x)+λ=C(x)< 0 (3.4)

where =(x), =G(x), =C(x) are the topological derivatives at the point x,
respectively of the Lagrangian, objective functional G (in the case of plane
elasticity problems it is denoted by =DG(x), and in the case of plates – by
=PG(x)) and cost functional C.
So, a new hole should be introduced at the point where the topological

derivative of the Lagrangian =(x) attains minimum and condition (3.4) is
fulfilled. This approach is similar to the bubble method, however apart from
the positioning of the hole it additionally formulates the condition of topology
modification.

Fig. 1. Idea of simultaneous finite topologymodification and shape optimization

The shape optimization of the new holes is a usually time-consuming and
complicated process. In order to reduce these difficulties, a concept of finite
modification is proposed (Fig.1). We can distinguish two types of problems
(cf. Bojczuk and Szteleblak, 2003, 2005). The first one corresponds to the
introduction of holes of a prescribed size and shape. Assuming that, approxi-
mately, the introduction of a finite void can be treated as the inserting of a
sum of infinitesimally small voids, the preliminary problem takes the form

∆L(opt) =min
r

∫

Ap

=(x) dAp (3.5)

where Ap denotes a fixed area of the new hole and r is the vector of its
position. The topology modification condition corresponds to negative value
of this increment, namely

∆L(opt) < 0 (3.6)

The second type of problems consists in the introduction of holes of an
unknown size and shape together with the introduction of finite changes of



334 D. Bojczuk, W. Szteleblak

other boundaries.As for the redesign path related to the constant cost C0 and
minimum of the objective function G, themultiplier λ achieves its minimum
at the optimal point, then for finite variations of design we can introduce the
design quality function Λ (cf.Mróz andBojczuk, 2003). Now, the preliminary
problem takes the form

min
A
Λ where Λ=

µλ+∆G
µ−∆C

=
µλ+

∫
A
=G(x) dA

µ−
∫
A
=C(x) dA

(3.7)

Here, A denotes the unknown area of the hole and µ (µ ­ 0) is the scaling
factor. So, the condition of acceptance of thefinite topologymodification takes
the form

Λ<λ (3.8)

When µ=0, problem (3.7) takes the form

min
A
Λ where Λ=−

∆G

∆C
=−

∫
A
=G(x) dA
∫
A
=C(x) dA

(3.9)

and the considered approach becomes analogous to the bubble method.
In order to solve problem (3.5) or (3.7) and check a proper condition for

the introduction of finite topology modification, the finite element discretiza-
tion is used. Then, the solution to the problem is relatively easy to obtain.
However, it gives reduction of the structure cost, so in order to follow the path
of the constant global cost we can, for example, apply the idea of proportional
variation of thickness. The obtained design usually differs only a little from
the global optimum.

3.2. Correction of structure shape

After the finite topologymodification, only an improved structure is obta-
ined, so anadditional standard shapeoptimization shouldbeperformed. It can
be done using optimality conditions (3.3) for optimization problem (3.1) with
an updated vector of design parameters. When an arbitrary gradient method
is used,we need sensitivity expressionswith respect to the shapemodification.
Now, we assume that shape variation of the domain Ω is given in the form

of a transformation f : Ω→Ω∗, namely

x
∗ =x+ δϕ (3.10)



Application of finite variations to topology... 335

where x∗ is the transformed position of the point x and, as previously,
δϕ = [δϕ1,δϕ2]> is variation of the shape transformation function. We al-
so assume that the shape transformation of the boundary ΓT is described by
variation of the shape parameters al (l = 1,2, . . . ,L). Then, variation of the
transformation function can be presented in the form

δϕk =
L∑

l=1

∂ϕk
∂al
δal k=1,2 (3.11)

Now, let us consider particular cases of the shape transformation used in
this paper. First, translation of a hole with a prescribed size and shape will
be taken into account. Then, we have

ϕk = ak k=1,2 (3.12)

where a1, a2 are components of the translation vector a. Finally, variation of
the transformation function takes the form

δϕ1 = δa1 δϕ2 = δa2 (3.13)

In order to provide a proper way of approximate boundary description,
also with sharp edges, B-splines can be used (cf. Gerald andWheatley, 1995;
Kiciak, 2005). The recursive definition of the normalized B-spline function of
the order k is

N
(1)
j (t)=

{
1 for tj ¬ t¬ tj+1
0 otherwise

(3.14)

and
N
(k)
j (t)=

t− tj
tj+k−1− tj

N
(k−1)
j (t)+

tj+k− t

tj+k− tj+1
N
(k−1)
j+1 (t) (3.15)

where t= [t1, t2, . . . , tm]> is the knot vector created as a set of non-decreasing
numbers. Now, a polynomial B-spline curve can be written as follows

ϕi(t)=
n∑

j=0

N
(k)
j (t)xi(j) i=1,2 (3.16)

where xj = [x1(j),x2(j)]
> are coordinates of the jth control point.

Themodification of the boundary can also be described byNURBS (non-
uniform rational B-splines). In this case, we have (cf. Adamski, 1997; Kiciak,
2005)

ϕi(t)=

∑n
j=0N

(k)
j (t)wjxi(j)

∑n
j=0N

(k)
j (t)wj

i=1,2 (3.17)

where wj is the weight of the point xj.



336 D. Bojczuk, W. Szteleblak

The possibility of control of the shape by only several control points and
the possibility of local changes of shape without any influence on the rest of
the structure are themain advantages of B-splines andNURBS. In this paper,
B-splines of the order k=4 (so-called cubic B-splines) are used. Here, in the
shape optimization process, the coordinates of control points are the design
variables. However, also the knot vector andweights of the control points can
be used to modify the shape.
Now, let us consider sensitivity analysis of the functional of strains and

displacements expressedby (2.2).Here, the adjoint structure specifiedby (2.4),
(2.5) is introduced. Let us denote arbitrary shape parameters by al (l =
1,2, . . . ,L). In the case of the consideredhere cubicB-splines, theseparameters
correspond to the coordinates of control points x1(j), x2(j), j = 1,2, . . . ,J,
L=2J. Assuming that the conservative surface traction T 0 is applied on ΓT ,
the variation of functional (2.2) can be written in the form (cf. Dems, 1980;
Kleiber, 1995)

δG=
∫

ΓT

2∑

k=1

L∑

l=1

[
F −σaε−2(f+T0ua)K+

∂

∂n
(f+T0ua)

]
nk
∂ϕk
∂al
δal dΓT +

(3.18)

+
∫

Ω

δp0ua dΩ

where, as previously, ∂/∂n denotes a derivative in the direction normal to
themodified boundary, K is the curvature of this boundary and δp0 denotes
local variation of volume forces referred to the initial shape.
Next, let us consider the functional of stresses and reaction forces expres-

sed by (2.23). Here, adjoint structure (2.24) with initial strains and non-zero
boundary conditions on Γu is introduced. Assuming the conservative surface
traction T0 on ΓT , the variation of functional (2.23) can be written in the
form (cf. Dems, 1980; Kleiber, 1995)

δG=
∫

ΓT

2∑

k=1

L∑

l=1

[
H−σεa+

∂

∂n
(T0ua)−2T0uaK

]
nk
∂ϕk
∂al
δal dΓT +

(3.19)

+
∫

Ω

δp0ua dΩ

For a particular case, when only the unloaded part Γ0 of the boundary ΓT
can be modified and f = 0 on Γ0, sensitivity expressions (3.18) and (3.19)
take the simplified forms, namely



Application of finite variations to topology... 337

δG=
∫

ΓT

2∑

k=1

L∑

l=1

(F −σaε)nk
∂ϕk
∂al
δal dΓT (3.20)

and

δG=
∫

ΓT

2∑

k=1

L∑

l=1

(H−σεa)nk
∂ϕk
∂al
δal dΓT (3.21)

However, sensitivity expressions (3.18)-(3.21) are formulated for plane elasti-
city problems, analogously as in Section 2.2, they can be easily adapted to
plate structures.
Moreover, let us consider the cost functional of the form

C =
∫

V

c dV (3.22)

where analogously to (2.38), c is a unit cost, and V denotes volume of the
considered structure. Then, variation of this functional can be written as fol-
lows

δC =
∫

ΓT

2∑

k=1

L∑

l=1

cnk
∂ϕk
∂al
δal dΓT (3.23)

Presented here sensitivity expressions can be used in an arbitrary gradient
method of shape optimization (shape correction).

4. Optimization algorithms and numerical examples

4.1. Optimal design by introduction of holes of prescribed size and shape

Let us assume that a hole (holes) with the prescribed shape and fixed area
Ap should be introduced into a plane structure. Now, the problem of the form
analogous to (3.1) is considered. It can be stated as a searching of such a
position of the hole for which the functional expressed by (2.2) or (2.23) with
constraints imposed on the global cost, attains minimum (maximum). Also,
other constraints, for example geometrical constraints, can be additionally
used. To solve the problem, a heuristic algorithm composed of two steps is
proposed here.
In the first step, the hole is introduced in a place where the integral of

the topological derivative over the domain of the hole is minimal (maximal).



338 D. Bojczuk, W. Szteleblak

So, in order to find the best position r of the hole, preliminary problem (3.5)
should be solved.
In the second step, standard optimization procedures are employed to cor-

rect the position of the hole. When a gradient method is used, the sensitivity
analysis with respect to horizontal and vertical translation of the hole is per-
formed. Now, sensitivity expressions presented in Section 3.2 can be applied.

4.1.1. Example 1: Optimization of square hole position

Let us assume that for some technological reasons a square hole of dimen-
sions 1m×1m should be inserted into the rectangular domain (4m×6m) of
a plane structure (Fig.2). This structure is made of steel (Young’s modulus
is E =2.1 ·105MPa and Poisson’s ratio ν =0.3). It is clamped on the lower
edge, loaded by linearly distributed shear forces on the upper edge, and its
initial thickness is h=30mm (h= const).

Fig. 2. Optimization of the position of a square hole (a) primary structure;
(b) adjoint structure

The optimization problem considered here is analogous to (3.1). The aim
is to find the hole position which minimizes the horizontal displacement uA
of the point A with a condition imposed on the global cost proportional to
thematerial volume.Also additional geometrical constrains,which assure that
the hole lies inside the domain of the structure, are used here. A dashed line
located 0.3m from the edges (Fig.2a), denotes this admissible domain.
Now, the objective functional can be written in the form

G=
∫

ΓT

δ(x−xA)u1(x) dΓT =uA (4.1)



Application of finite variations to topology... 339

where δ(x−xA) denotes Dirac’s delta and u1(x) is the field of horizontal
displacements. In order to solve the considered problem, an adjoint system
of the same geometry and boundary conditions as the primary structure is
introduced. Moreover, it is loaded at the point A by a unit force Pa = 1
acting in accordance with the direction of displacement uA (Fig.2b).
The optimal position of the hole, obtained directly after solution of preli-

minary problem (3.5), is shown in Fig.2a. Now, both vertical and horizontal
geometrical constraints are active, and correction of the hole position is not
necessary. The ratio of the horizontal displacements uA of the initial (witho-
ut hole) and optimal design is u

(init)
A /u

(opt)
A = 1.002, and the final thickness

arising from the cost (volume) condition equals h = 31.304mm. So, the in-
troduction of the technological hole has not caused the effect of increase of
the analyzed displacement and – on the contrary – even a small decrease has
appeared here.

4.1.2. Example 2: Optimization of circular hole position

Letusassumethat for technological reasonacircularhole of 0.5mdiameter
should be introduced into the rectangular plate (3m×2m) shown in Fig.3.
The structure is made of aluminium (Young’s modulus is E = 75GPa and
Poisson’s ratio ν = 0.35). The initial thickness of the plate is 15mm. The
structure is clampedon three edges and the fourth (upper) edge is free (Fig.3).
The transverse load changes linearly along height of the plate.

Fig. 3. Geometry of the plate

The aim of the optimization process is to find the hole position whichmi-
nimizes the strain energy of the structure with a condition imposed on the



340 D. Bojczuk, W. Szteleblak

global cost proportional to the material volume. Additionally, geometric con-
straints,whichassure that thehole lies inside thedesigndomain, are usedhere.
A dashed line situated 0.15m from the edges (Fig.4) denotes this admissible
domain. Here, the objective functional, which expresses the strain energy, can
be written in the form

G=
1
2

∫

Ω

Mκ dΩ (4.2)

where Ω is the middle surface of the plate, M = [M1,M2,M12]> is the
moment vector and κ= [κ1,κ2,κ12]> denotes the curvature vector. Now, the
adjoint structure is the same as the primary one. Two types of modification
are analyzed here.

Fig. 4. Optimal position of the circular hole: (a) with unchanged load; (b) with
removed load

In the first case, it is assumed that the overall load does not change. Howe-
ver, its part initially applied to the removed hole is replaced by the equivalent
load distributed around it. It corresponds to the situation where the intro-
duced hole is closed by a circular cover simply supported on the boundary.
Now, in order to find the optimal position of the hole, at first preliminary pro-
blem (3.5) is solved, where the topological derivative expressed by (2.36) or
(2.37) is used. The obtained design is presented in Fig.4a by the dashed line.
Next, correction of the hole position is performed. The optimal hole position
is denoted by the continuous line, where coordinates of the hole center are
x
(opt)
1 = 0.792m, y

(opt)
1 = 1.548m, and now geometrical constraints are not

active. The ratio of the strain energy of the initial (without hole) and optimal
design is G(init)/G(opt) = 1.073 and the final thickness arising from the cost
condition equals h = 15.507mm. Due to symmetry of the structure and lo-
ading, also the second, an equivalent solution with the symmetric position of
the hole, namely x

(opt)
2 =−0.792m, y

(opt)
2 =1.548m exists.

In the second case, it is assumed that the load is removed together with
insertion of the hole. Now, the topological derivative expressed by (2.36) or



Application of finite variations to topology... 341

(2.37) should be completed by an additional term, namely (−p0u). Then,
solving preliminary problem (3.5), we get the initial position of the hole,
which is presented in Fig.4b by the dashed line. The optimal position of
the hole, denoted by the continuous line, is obtained by standard optimiza-
tion. The coordinates of the hole center, which lies on the symmetry axis, are
x
(opt)
3 = 0.000m, y

(opt)
3 = 0.785m, and the ratio of the strain energy of the

initial and optimal design is G(init)/G(opt) =1.274.

4.2. Optimal design by introduction of finite topology and shape modifi-

cations

Let us consider the optimization problem of form (3.1). Here, it can be
treated as the search for such a size and shape of a hole (holes) and shape of
variable boundaries for which the functional expressed by (2.2) or (2.23) with
constraints imposed on the global cost attains minimum (maximum). Also,
other constraints, for example geometrical constraints, can be additionally
used. In order to solve the problem, a heuristic algorithm for simultaneous
topology and shape optimization composed of two mutually interacted stages
is proposed here.

Fig. 5. Flowchart of the optimization process by finite topology and shape
modification

At thefirst stage, thefinitemodificationby introductionof ahole or change
of a variable boundary is applied. In order to choose the bestmodification, the



342 D. Bojczuk, W. Szteleblak

preliminary problem of form (3.7) should be solved. Quantity of the removed
boundary is controlled by the parameter µ (µ­ 0).When the objective func-
tional after finitemodification does not decrease, the value of the coefficient µ
can be reduced.This procedure is repeated until it is possible to introduce any
finite modification.
At the second stage, standard shape optimization procedures are used to

smooth the variable boundaries. The boundaries of new holes are described
by some additional shape parameters. So, the problem of form analogous to
(3.1), but with an updated vector of design variables should be solved. It can
be done using an arbitrary gradient method where optimality conditions are
expressed by (3.3) and respective sensitivity formulas are presented in (3.18)-
(3.21) and (3.23). Next, especially when huge changes are inserted during the
shape optimization process, we can return to the first stage.
When any modification is not possible to introduce, the process of simul-

taneous topology and shape optimization is finished. A general flowchart of
the optimization process is shown in Fig.5.

4.2.1. Example 3: Optimization of topology and shape of plane bridge structure

Let us consider optimal design problem (3.1) for a plane structure shown
in Fig.6a, where G corresponds to the strain energy U and the cost C is
proportional to the material volume. Also additional geometrical constraints
are used here, namely it is assumed that the domain above the dashed line
can not be removed.

Fig. 6. Optimal design of the bridge structure (µ=0.6): (a) initial structure;
(b)-(e) every fifth successive finite modification; (f) final design



Application of finite variations to topology... 343

The structure is simply supported on a part of the lower edge and lo-
aded by vertical forces on the upper edge. It is assumed that it has a con-
stant thickness in the whole domain, which can change during the opti-
mization process. The material properties are as follows: Young’s modulus
E = 1.5 · 104MPa, Poisson’s ratio ν = 0.167. The structure is divided into
3000 finite elements. Two cases are considered here.
In the first case, the scaling factor µ = 0.6 is assumed. The history of

optimization is shown in Fig.6. The optimal structure (Fig.6f) is obtained
after 18 finite modifications and final shape optimization. The ratio of the
strain energies of the initial and optimal design is U (init)/U(opt) =1.406.
In the second case, the scaling factor µ = 3.0 is assumed. The pro-

cess of optimization is shown in Fig.7, and now the optimal structure
(Fig.7e) is obtained only after 7 finite modifications and final shape opti-
mization. The ratio of the strain energies of the initial and optimal design is
U(init)/U(opt) =1.456.

Fig. 7. Optimal design of the bridge structure (µ=3.0): (a) initial structure;
(b)-(e) every second successive finite modification; (f) final design

It is important to notice that the assumption of a bigger value of the
scaling factor µ accelerates the optimization process. However, when a too
large value of µ is taken, it may lead to some difficulties – for example to the
loss of structure connectivity. The coefficient µ should be chosen andmodified
by the user.Moreover, the optimal solution is not unique (compare Fig.6f and
Fig.7e), and we may expect local optima.

4.2.2. Example 4: Optimization of topology and shape of plane, beam-like structure

The next application is the optimal design of a plane, beam-like structure
shown in Fig.8a. The optimization problem is the same as in the previous



344 D. Bojczuk, W. Szteleblak

example, i.e. to minimize the strain energy with a condition imposed on the
cost, where the cost is proportional to the material volume. The structure is
made of steel with Young’s modulus E = 2.1 · 105MPa and Poisson’s ratio
ν =0.3. It is divided into 2880 finite elements.

Fig. 8. History of optimization: (a) initial structure; (b)-(e) successive finite
modifications; (f) final design

The history of optimization is shown in Fig.8. The optimal structure, pre-
sented in Fig.8f, is obtained after 5 finitemodifications and final correction of
the shape. The ratio of the strain energies of the initial and optimal design is
U(init)/U(opt) =1.275.

4.2.3. Example 5: Optimization of topology and shape of simply supported plate

The rectangular plate (200mm×100mm) shown in Fig.9 is analyzed. The
structure is made of steel. It is simply supported on each edge and loaded by
a uniformly distributed load, non-symmetrically located near the center of the
structure. The initial thickness of the plate is 5mm.

Fig. 9. Geometry of the plate, boundary and loading conditions

Now, let us consider again optimal design problem (3.1). Here G corre-
sponds to the strain energy U and the cost C is proportional to thematerial
volume. Themaximum value of the plate thickness is limited to 20mm.



Application of finite variations to topology... 345

Fig. 10. History of optimization: (a)-(c) every second iteration of thickness reduction

Fig. 11. History of optimization: (a)-(c) every second iteration of middle layer
removal; (d) final design obtained after shape correction

When, in the optimization procedure thematerial is removed all along the
thickness of the plate, often non-connectivity of the structure appears. This
situation, for example, takes place near simply supported edges. In order to
avoid these difficulties, the process of material removal is divided into two
phases. In the first phase, using condition (3.7) based on the topological deri-
vative, the thickness of the plate is symmetrically reduced by 2/3, see Fig.10,
where the lighter colour denotes the smaller thickness. It is a pre-selection of
the domain, which may be completely removed. In the second phase, using
condition analogous to (3.7), but based on the shear energy, only thematerial
fromthe thin layer is successively reduced (Fig.11a,b,c). Next, the thickness of
thenon-removeddomain is equalized. Finally, shapeoptimization is performed
and the optimal structure is shown in Fig.11d. Constraints imposed on thick-



346 D. Bojczuk, W. Szteleblak

ness of the structure are active. During the optimization process, the strain
energy was reduced from 2.107J down to 0.0811J (U (init)/U(opt) =26.292).

4.2.4. Example 6: Optimization of topology and shape of clamped plate

Let us consider optimal design problem (3.1) for a plate structure shown in
Fig.12. As previously, G corresponds to the strain energy U and the cost C is
proportional to the volume of the structure. It is assumed that areas bounded
by the broken line and denoted as the ”passive domain” can not be removed.
Moreover, a condition on the maximal thickness of the plate is imposed.

Fig. 12. Geometry of the plate

The rectangular plate (150mm×100mm) is made of aluminium with
Young’s modulus is E = 75GPa and Poisson’s ratio ν = 0.35. Its initial
and maximal thicknesses are respectively 5mm and 15mm. The structure is
uniformly loaded on all edges of the external rectangle by transverse forces
of intensity q = 5kN/m. The plate is clamped along the non-symmetrically
located internal rectangle, marked by the broken line in Fig.12.
Here, analogously as in the previous example, the process of material re-

moval is divided into two phases. In the first phase, thickness of the plate is
symmetrically reducedby 2/3. Figure 13 illustrates this step,where the darker
colour denotes the total thickness. In the second phase, thematerial from the
central layer is successively reduced (Fig.14a,b,c). Next, shape optimization of
the structure with equalized thickness is performed. The optimal structure is
shown in Fig.14d. The constraints imposed on themaximum thickness of the
structure are active. The rtio of strain energies before and after optimization
is U(init)/U(opt) =6.585.



Application of finite variations to topology... 347

Fig. 13. History of optimization: (a)-(d) every second iteration of thickness reduction

Fig. 14. History of optimization: (a)-(c) every second iteration of thin layer removal;
(d) optimal design obtained after shape correction

5. Concluding remarks

A heuristic algorithm of simultaneous topology and shape optimization
which uses finite topologymodifications is presented in the paper. It is applied
to optimal design of 2D structures working in a plane state of stress and
beingKirchhoff’s plates. It is important to notice that the application of finite
modifications essentially reduces computation time required for generation of



348 D. Bojczuk, W. Szteleblak

improved or optimal designs. Another advantage of this method arises from
the natural way of evolution of the optimal design. Here, the optimization
process can be stopped at any level of the structure complexity, and usually
the objective functional only slightly differs from the global minimum.
Numerical examples shown in the paper confirm the applicability and use-

fulness of the approach.

Acknowledgement

Theworkwas partially supported by the StateCommittee for ScientificResearch;
KBN grant No. 4T07E02328.

References

1. Adamski W., 1997, Application of NURBS in numerical modeling of objects,
Proceedings of the XIII Polish Conference on ComputerMethods inMechanics,
Garstecki A., Rakowski J. (edit.), Politechnika Poznańska, 1, 89-97

2. Bojczuk D., Szteleblak W., 2003, Design optimisation for plane elasticity
problems using finite topology variations,Proceedings of the 15th International
Conference on Computer Methods in Mechanics, Wisła, CDROM

3. Bojczuk D., Szteleblak W., 2005, Topology and shape optimization of
plates using finite variations,Proceedings of the 16th International Conference
on Computer Methods in Mechanics, Częstochowa, CDROM

4. Burczyński T. (edit.), 2002, Computational sensitivity analysis and evolu-
tionaryoptimizationof systemswithgeometrical singularities,ZeszytyNaukowe
Katedry Wytrzymałości Materiałów i Metod Komputerowych Mechaniki, Poli-
technika Śląska, Gliwice

5. Burczyński T., Kokot G., 2003, Evolutionary algorithms and boundary
elementmethod in generalized shape optimization, J. Theor. and Appl. Mech.,
41, 2, 341-364

6. Czarnecki S., Dzierżanowski G., Lewiński T., 2004, Topology optimi-
zation of two-component plates, shells and 3D bodies, Optimal shape design
and modelling, Selected papers presented at WISDOM 2004, Lewiński T., Sig-
mund O., Sokołowski J., Żochowski A. (edit.), Akademicka OficynaWydawni-
cza EXIT,Warszawa, 15-29

7. DemsK., 1980,Wieloparametrowaoptymalizacjakształtu konstrukcji,Zeszyty
Naukowe, 371, Politechnika Łódzka, Łódź

8. Eschenauer H.A., KobelevV.V., SchumacherA., 1994, Bubblemethod
for topology and shape optimization, Struct. Optim., 8, 42-51



Application of finite variations to topology... 349

9. GeraldC.F.,Wheatley P.O., 1995,Applied Numerical Analysis, Addison-
Wesley Publishing Company

10. Il’in A.M., 1992,Matching of asymptotic expansions of solutions of boundary
value problems,Translations of Mathematical Monographs, 102, AMS

11. Kiciak P., 2005,Podstawy modelowania krzywych i powierzchni, WNT, War-
szawa

12. Kleiber M. (edit.), 1995, Komputerowe metody mechaniki ciała stałego,
PWN,Warszawa

13. Marczewska I., Sosnowski W., Marczewski A., Bednarek T., 2003,
Topology and sensitivity – based optimization of stiffened plates and shells,
Short Papers of the Fifth World Congress of Structural and Multidisciplinary

Optimization, 271-272

14. Mróz Z., Bojczuk D., 2003, Finite topology variations in optimal design of
structures, Struct. Multidisc. Optim., 25, 153-173

15. Sokołowski J., Żochowski A., 1999, On topological derivative in shape
optimization, SIAM J. Control and Optimiz., 37, 4, 1251-1272

Zastosowanie modyfikacji skończonych w optymalizacji topologii

i kształtu konstrukcji dwuwymiarowych

Streszczenie

Wpracy rozpatrywana jestmetoda jednoczesnej optymalizacji topologii i kształtu
konstrukcji dwuwymiarowychprzy użyciu skończonychmodyfikacji topologii. Rozwa-
żania dotyczą zarówno konstrukcji tarczowych pracujących w płaskim stanie naprę-
żenia, jak i płyt Kirchhoffa pracujących w stanie zgięciowym. Przy wykorzystaniu
pochodnej topologicznej wyprowadzonowarunki wprowadzania skończonychmodyfi-
kacji topologii. Gdy spełniony jest odpowiedni warunek modyfikacji, do konstrukcji
wprowadzane są otwory o skończonychwymiarach oraz ewentualnie skończonemody-
fikacje pozostałych brzegów. Następnie wykonywana jest standardowa optymalizacja
kształtu otworów i brzegów zewnętrznych.Analizowane są dwa podstawowe typymo-
dyfikacji, a mianowicie wprowadzanie otworów o zadanej wielkości i kształcie oraz
wprowadzanie otworówo nieznanej wielkości i kształcie wraz z ewentualną skończoną
zmianą pozostałych brzegów. W pracy sformułowano odpowiedni algorytm heury-
styczny optymalizacji topologii i kształtu rozpatrywanychkonstrukcji. Przedstawione
przykłady ilustracyjne potwierdzają przydatność zaproponowanego podejścia.

Manuscript received November 14, 2005; accepted for print February 6, 2006