PLACE IN RETURN BOX to remove th'
To AVOID FINES return on

MAY BE RECALLED with earlier

IS checkout from your record.
or before date due.
due date if requeSted.

 

DATE DUE

DATE DUE

DATE DUE

 

NI

MW ’-

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mo amour-cum“

SOLVING TOPOLOGY OPTIMIZATION PROBLEMS USING
WAVELET-GALERKIN TECHNIQUES

By

Giuseppe C. A. DeRose Jr.

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mechanical Engineering

1998

ABSTRACT
SOLVING TOPOLOGY OPTIMIZATION PROBLEMS USING
WAVELET-GALERKIN TECHNIQUES
By
Giuseppe C. A. DeRose Jr.

The topology optimization problem in two— and three-dimensional elasticity is
solved with a meshless, wavelet-based solution scheme. A fictitious domain approach
is used to embed the design domain into a simple regular domain: a square in two
dimensions and a cube in three dimensions. The material distribution and displace-
ment field are discretized over the fictitious domain using fixed-scale, Shift-invariant
wavelet expansions. A discrete form of the elasticity problem is solved using a wavelet-
Galerkin technique during each iteration of the topology optimization sequence. Ap—
proximate solutions are found with an efficient Preconditioned Conjugate Gradient
(PCG) solver using a new class of non-diagonal preconditioners. Elastostatic analysis
and topology optimization examples Show convergence rates that are insensitive to
the level of resolution. The convergence and memory management properties of the
FCC algorithm make the analysis of large-scale problems possible. Several topology

optimization examples in two and three dimensions are included.

To my parents

iii

ACKNOWLEDGMENTS

I would like to express my sincere gratitude to my advisor Dr. Alejandro Diaz. Your
guidance made this research possible. I would also like to thank the members of my
committee: Dr. Ronald Averill, Dr. André Benard, and Dr. Michael Frazier. Thank
you for the time spent in answering my questions and reviewing this document.

I am deeply indebted to Dr. Ronald Rosenberg; your support is truly appreciated.
A special thanks is extended to Dr. Craig Somerton; your time for discussions about
attending graduate school was invaluable.

I would also like to thank my fellow graduate students: Charles Birdsong, Brooks
Byam, Michael Hales, Carlos LOpes, and Mark Minor. Your support during this
adventure is duly noted.

Lastly, I thank my family. I thank you for patiently waiting for the completion of

this work, understanding my struggles, and supplying my main source of motivation.

iv

TABLE OF CONTENTS

LIST OF TABLES viii
LIST OF FIGURES ix
1 Introduction 1
1.1 Goals ..................................... 2
1.2 Approach ................................... 3
1.3 Layout of Dissertation ............................ 4
2 TOpology Optimization 5
2.1 Background .................................. 5
2.2 Formulations of Topology Optimization ................... 8
2.2.1 Macroscopic TOpology Optimization ................... 13
2.2.2 Topology Optimization Using Homogenization .............. 14
2.2.3 TOpology Optimization Using A Fictitious Material Model ....... 18
2.3 Solution Methods ............................... 20
2.4 Summary ................................... 21
3 Wavelet Approximations in Numerical Analysis 23
3.1 Introduction .................................. 23
3.2 Multiresolution Analysis ........................... 25
3.3 Defining Wavelets ............................... 27
3.4 Defining Wavelet Bases ............................ 28
3.5 Examples of Orthogonal Wavelets ...................... 31
3.5.1 Haar Wavelets ............................... 31
3.5.2 Daubechies D6 Wavelets .......................... 33
3.6 Examples of Non-Orthogonal Wavelets ................... 34
3.6.1 Linear Spline Wavelets ........................... 35
3.6.2 Quadratic Spline Wavelets ......................... 36
3.7 Wavelet Representations and Algorithms .................. 37
3.8 Wavelets In Multiple Dimensions ...................... 41
3.9 Wavelet Applications ............................. 43
3.9.1 Function approximations .......................... 43
3.9.2 Image Compression ............................. 43
3.9.3 Wavelet-Based PDE Solvers ........................ 47
3.10 Summary ................................... 50

4 Fictitious Domain Approach Using Wavelet-Galerkin Techniques

For Solving Partial Differential Equations 51
4.1 Background .................................. 51
4.2 Solving Two-Dimensional Elasticity Problems Using a Fictitious Domain

Approach and Wavelet Bases ....................... 52
4.2.1 Problem Statement ............................. 53
4.2.2 Wavelet Discretization ........................... 58
4.2.3 Constructing the Discrete System of Equations ............. 66
4.2.4 Alternative Fictitious Domain Problem Formulation .......... 79
4.2.5 Preconditioned Conjugate Gradient Solution Technique ......... 82
4.3 Post Processing Wavelet-Based Solutions .................. 86
4.4 Summary ................................... 88
5 Wavelet-Based TOpology Optimization in Two Dimensions 89
5.1 Introduction .................................. 89
5.2 Elastostatic Analysis Using Wavelet-Based Analysis ............ 90
5.2.1 Accuracy Of The Wavelet-Based Fictitious Domain Solution ...... 93

5.2.2 Performance Of The Wavelet-Based Fictitious Domain Solution Method 98
5.3 Solving Topology Optimization Problems Using Wavelet-Based Analysis 103

5.4 Adjusting Wavelet-Based Analysis For Topology Optimization ...... 104
5.4.1 Effect Of Wavelet Used For Approximating The Displacement Field:
Checkerboard Patterns ...................... 105
5.4.2 Effect Of The Strength Of The Weak Material .............. 109
5.4.3 Effect Of The Extent Of The Fictitious Domain ............. 112
5.5 Performance Of Wavelet-Based Topology Optimization .......... 113
5.5.1 Effects Of Varying Resolution ....................... 115
5.5.2 Effects Of Varying The Extent Of The Fictitious Material ....... 116
5.6 Wavelet-Based Topology Optimization Examples ............. 117
5.6.1 Single Load Case Problems ........................ 118
5.6.2 Multiple Load Case Problems ....................... 125
5.7 Summary ................................... 132
6 Wavelet-Based Topology Optimization in Three Dimensions 133
6.1 Introduction .................................. 133
6.2 Fictitious Domain Problem Statement ................... 133
6.3 Problem Discretization In Three Dimensions ................ 134
6.3.1 Approximation Of The Material Distribution And Displacement Field . 135
6.3.2 The Discretized Equations ......................... 138
6.4 Post Processing Three-Dimensional Solutions ............... 143
6.5 Performance Of The Fictitious Domain Approach ............. 144
6.5.1 Effect Of Varying Resolution ....................... 145
6.5.2 Effect Of Varying The Amount Of Fictitious Material .......... 149
6.6 A More Memory Efficient Solution Technique ............... 150
6.7 Three-Dimensional Topology Optimization Examples ........... 153
6.7.1 Single Load Case Examples ........................ 154

vi

6.7.2 Multiple Load Case Examples ....................... 158

6.8 Summary ................................... 162
7 Conclusions 165
7.1 Summary ................................... 165
7.2 Areas Of Future Work ............................ 167
BIBLIOGRAPHY 168

vii

LIST OF TABLES

5.1 Discretization levels for analysis of FCC performance ........... 101

6.1 Discretization levels for analysis of FCC performance ........... 146

viii

2.1
2.2
2.3
2.4

2.5
2.6
2.7
2.8

3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9

LIST OF FIGURES

Example of a sizing optimization problem .................
Example of a fixed-topology Optimization problem .............
Example of a generalized topology optimization problem .........
Example of a topology optimization problem using a structure character-

istic function x(x) ............ ' ................
Example of forming an optimal design when approximating x(x) .....
Two-dimensional microcell in macroscopic materials ............
Effect of material exponent p ........................
Effect of discretization level on optimal shape resolution .........

Sample Daubechies wavelets .........................
Effects of dilation and translation of wavelets ...............
Example Haar scaling function and wavelet ................
Example Daubechies D6 scaling function and wavelet ...........
Example linear spline scaling function and wavelet .............
Example quadratic spline scaling function and wavelet ..........
Multiresolution using haar wavelets .....................
Multiresolution using Daubechies D6 wavelets ...............
Original component image ..........................

3.10 Component image created from 5% of the wavelet coefficients ......
3.11 Wavelet based approaches for solving ordinary and partial differential

4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8
4.9

equations .................................

An example topology optimization problem statement ...........
Package space embedded into a fictitious domain .............
Level 4 discretization of an example topology optimization problem . . .
Example of using Haar approximations of a circle at various levels . . . .
Example of edge pixel formation of Aw ...................
Structure of A when using Daubechies D6 fixed scale expansions .....
Structure of A when using quadratic spline S4 fixed scale expansions
Structure of A when using linear Spline S3 fixed scale expansions
Structure of A], when using Daubechies D6 fixed scale expansions . . . .

4.10 Structure of Ah when using quadratic spline 84 fixed scale expansions . .
4.11 Structure of Ah when using linear spline 83 fixed scale expansions . . . .

5.1
5.2

Example elastostatic problem definition ...................
Example problem placed in a periodic fictitious domain ..........

RIG)

9
12
15
19
21

25
3O
32
33
36

44
45
46

48

53
55
58
60
63
69
70
71
74
75
76

90
91

5.3 Wavelet-based discretization ......................... 92

5.4 Finite element discretization ......................... 93
5.5 The areas of comparison for the wavelet—based and finite element dis-

cretizations ................................ 94
5.6 Norm of the traction along the top and right edge for D6 wavelet-based

and finite element solution techniques .................. 96
5.7 Norm of the traction along the right edge for a 192 x 120 pixel discretization

and D6 wavelet-based solution ...................... 98
5.8 Norm of the traction along the top and right edge for various wavelets and

E2 = 10”4 ................................. 99

5.9 The performance Of the PCG solver using various wavelet approximations 100
5.10 Performance of the PCG solver for wavelet-based methods using S3 spline

wavelets and traditional iterative finite element methods ....... 102
5.11 Optimal topologies using various wavelet approximations for the displace-

ment field ................................. 106
5.12 Material mixture in layer form for two approximations .......... 107
5.13 Optimal topologies found when approximating the material field one level

lower than the displacement field .................... 110
5.14 The effect of the strength of the fictitious material on the optimal topology 111
5.15 Bracket problem definition .......................... 113
5.16 Optimal topologies for varying amounts of fictitious material ....... 114

5.17 Solutions of the 8x5 design domain problem for varying levels of resolution 115
5.18 Number of PCG iterations required for solution of the 8x5 design domain

at varying levels of resolution ...................... 116
5.19 Performance of the PCG solver for various design domain Sizes ...... 117
5.20 Problem statement for truss example 1 ................... 118
5.21 Results of topology Optimization for truss example 1 ........... 119
5.22 Problem statement for truss example 2 ................... 121
5.23 Results Of topology optimization for truss example 2 ........... 122
5.24 Problem statement for the L domain example example .......... 123
5.25 Results of topology optimization for the L domain example ........ 124
5.26 Problem statement for the five bar example ................ 126
5.27 Results of topology optimization for the five bar example ......... 127
5.28 Problem statement for the bridge example ................. 128
5.29 Results of topology optimization for the bridge example .......... 129
5.30 Problem statement for the dual load 8x5 design domain example ..... 130
5.31 Results of topology optimization for the 8x5 domain dual load example . 131
6.1 Example voxel discretization ......................... 135
6.2 Example problem statement ......................... 145
6.3 Optimal topologies found for various discretizations ............ 147
6.4 PCG convergence properties for various discretizations .......... 148
6.5 PCG convergence for the initial topology Optimization iteration ..... 149
6.6 Problem statement for fictitious material investigation .......... 150
6.7 Results of topology optimization for various cube discretizations ..... 151

6.8 PCG solver performance for various cube discretizations ......... 152

6.9 Problem statement for truss example 1 ................... 154
6.10 Results of topology optimization for truss example 1 ........... 155
6.11 Problem statement for truss example 2 ................... 156
6.12 Results of topology optimization for truss example 2 ........... 157
6.13 Problem statement for the X-truss example ................ 159
6.14 Results Of topology Optimization of the X—truss example for various volume
constraints ................................. 160
6.15 PCG performance for various volume constraints for the X-truss example 161
6.16 Problem statement for the bridge example ................. 162
6.17 Results of the topology optimization for the multiple load case bridge
example .................................. 163

xi

Chapter 1

Introduction

The goal Of topology or layout optimization is to determine the shape Of an optimal
structure that supports a given load, constrained only by the amount of material
available and restricted spatially to be within a prescribed package space. One of the
main benefits of topology optimization is that the general Shape and connectivity Of
the design are not specified a priori. Because of this, topology Optimization has be-
come an integral part of conceptual design in several engineering contexts. Engineers
use topology optimization to aide in the design of mechanical components creating
stronger, lighter, and more cost effective designs.

The use of topology optimization often requires the use of large-scale models.
Topology optimization problems are usually solved using an iterative Optimization
algorithm. During each iteration of the optimization sequence, the governing equa-
tions are solved to determine the performance of the current design. Typically, finite
element methods are used for this purpose. In these approaches, one finite element is

used per “shape” design variable. Thus, when a high resolution representation of the

shape is desired, a very large number of elements is required. This implies that when
topology optimization is used in a detailed design setting, very large-scale finite ele-
ment models must be solved numerous times (at least once per topology Optimization
iteration). As the problem size increases, the memory requirements and CPU time
necessary to obtain finite element solutions dominate the overall process. For large
scale problems, high computational demands Often make the problems impractical to
solve using reasonable computational resources.

One approach to reduce some of the computational resources needed to solve
large-scale topology optimization problems is to use an iterative solution technique.
Iterative solvers need more CPU time than direct methods but use substantially less
in-core memory. Unfortunately, when finite elements are used, the condition number
of the system matrices increases with the size of the problem and eventually iterative
schemes fail to converge due to poor conditioning. This implies that the amount
of shape detail available in the design domain is limited by the analysis technique.
This dissertation presents a new method to solve topology optimization problems by
replacing finite element analysis with meshless techniques that are specially tailored

to solve these problems efficiently.

1.1 Goals

We wish to develop a new analysis technique to determine approximate solutions to
the elasticity equations associated with topology optimization problems. We seek an

analysis technique that

1. performs analysis at rates that are independent of resolution

and
2. may be easily incorporated into a topology optimization environment.

To meet these goals, we turn to meshless, wavelet-based, Galerkin techniques.

1.2 Approach

Our approach is based on meshless methods for analysis that accurately predicts the
elastic response of highly detailed mechanical components. Discretization is devel-
oped from an image of the component embedded into a simple fictitious domain: a
square in two dimensions and a cube in three dimensions. The fictitious domain is
used to simplify the analysis, embedding the potentially complex package space in
a simple domain. In two dimensions, this process corresponds to a pixel level dis-
cretization of the component, while in three dimensions the discretization is similar
to a voxel discretization (the three—dimensional equivalent to pixels.) Given the im-
age discretization, wavelet bases defined at the resolution of the component image
are introduced in a traditional Galerkin scheme. This approach leads to a discrete
system of equations, similar to the stiffness matrix and load vector in finite element
methods. Due to the potentially high resolution of the image discretization, an iter-
ative preconditioned conjugate gradient solver is incorporated to solve the resulting
linear equations. A special preconditioner that results in convergence rates that are

insensitive to the resolution of the image discretization is used. This approach is

easily integrated with traditional topology optimization techniques. Like finite ele-
ments, each pixel or voxel in the image discretization is associated with a single shape
variable. Hence, this approach requires minimal modifications to standard topology

optimization algorithms.

1.3 Layout of Dissertation

The remainder of this dissertation is as follows. Chapter 2 discusses various formu-
lations and solution strategies for the topology optimization problem and current
limitations in these methods. Chapter 3 discusses wavelets and their use in numerical
analysis. Chapter 4 presents a new method of solving tOpology optimization problems
using a fictitious domain approach with wavelet-Galerkin solution techniques. Various
examples are solved using the approach developed in Chapter 4 and are discussed in
Chapter 5. These ideas are expanded for three-dimensional problems and presented

in Chapter 6. Conclusions about this new method are provided in Chapter 7.

Chapter 2

Topology Optimization

2.1 Background

Structural optimization has been addressed using three types of problem formula-
tions: sizing problems, shape problems, and generalized topology or layout problems.
Sizing problems focus on a structure with fixed connectivity. Section properties (e.g.,
thickness) are varied to Optimize the structure’s performance. A typical sizing prob-
lem consists of a truss with a fixed number of members and a fixed connectivity,
where the cross sectional area of the truss members are varied to maximize the struc-
ture’s stiffness. An example of this type of problem is shown in Figure 2.1. Shape
optimization problems are problems where boundaries may change shape, but no new
boundaries are introduced. Figure 2.2 illustrates an example of a shape optimization
problem consisting of a loaded plate with a hole. Minimization of the (maximum)
stress may be affected by changing the size and location of the hole, but no new holes

are introduced. Generalized topology optimization problems, introduced by BendsOe

 

 

 

 

 

 

 

 

\\\\\\\\\\\
\\\\\\\ \ \

(a) Original truss members (b) Optimal truss members

Figure 2.1: Example of a sizing optimization problem

and Kikuchi [3], allow the section properties and the layout of the structure to vary

during the optimization process. Initially, only the

1. design domain

2. boundary conditions

3. loading conditions

4. amount of material allowed in the final design

are specified. The design domain, or package space, refers to the region where
the design is created. The result of generalized topology optimization is the
optimal distribution of material within the design domain. Usually, optimal

corresponds to a structure that minimizes compliance or maximizes stiffness.

 

Plate

 

\\\\\\\

 

 

(a) Original shape

 

Plate

 

\\\\\\\

 

 

(b) Optimal shape

Figure 2.2: Example of a fixed-topology Optimization problem

 

 

 

 

/ Design
? Domain
/ co
(a) Problem definition (b) Optimal topology

Figure 2.3: Example of a generalized topology optimization problem

This definition will be used throughout this dissertation. Note that topology opti-

 

mization produces a shape and a connectivity. The main benefit of generalized topol-
ogy optimization is that the general shape and connectivity of the optimum design are
not specified a priori but are determined from the Optimization. Figure 2.3 illustrates
an example of a generalized topology optimization problem. Generalized topology
Optimization problems are often referred to synonymously as topology optimization or
layout optimization.

The remainder of this chapter will discuss various topology optimization formula-

tions and solution methods.

2.2 Formulations of Topology Optimization

This section discusses various mathematical formulations of the topology Optimization
problem. Again, topology optimization is used to determine the optimal distribution

of a fixed amount of material in a design domain for a given set of boundary and

 

 

 

 

 

/

/

/ Design

/ Domain

/

/

/ (o

(a) Problem statement (b) Optimal topology

Figure 2.4: Example of a topology optimization problem using a structure character-
istic function x(x)

loading conditions. This description of the material may be characterized in several
different ways. For instance, let w, C w be a set that describes the layout of the

structure. We may then define the structure indicator function as

I if x E w,
X(X) = (2.1)
0 if x E w \ w,
x(x) is the structure indicator function that characterizes the layout of the design.

Figure 2.4 illustrates the topology Optimization problem statement, optimal topology,

and structure characteristic function. With this notation, the topology Optimization

problem may be stated as

Find x(x), x E w that
minimizes f truly

I (2.2)
subject to fw x(:r)dw g meas(w) x V

Equilibrium equations (e.g., plane elasticity)

where w is the design domain, t represents the applied loading (i.e., tractions), 'y C 8w
is the portion of Ba) where the loads are applied, u is the displacement field, and V is
the maximum percentage of the design domain occupied by material. The equilibrium

equations correspond to the field equations of elastostatics:

1. Equilibrium equations:

 

 

 

0350' + fi = 0 (2.3)
2. Strain-displacement relationship:
1
(Sq = 5 (Ugj + UL.) (2.4)
3. Stress-strain relationship:
0}} = DijkmEkm (2.5)

The equilibrium equations relate the state u and the design x: the elasticity tensor

Dijkm in the stress-strain relationship depends on the current design.

10

Solving the topology Optimization problem as stated in (2.2) is difficult. In the
literature, there have been three main solution strategies suggested: macroscopic
topology optimization, topology Optimization using homogenization, and topology
optimization using a fictitious material model. The strategy behind these approaches
relies on replacing x(x) with a discrete approximation of x(x), denoted X A(x). It is

typical to approximate the structure indicator function with an expansion such as

X(X) z xirx) = int-(x) (2.6)

i=1

where {q,-(x) }§:’1‘ is suitable a basis. Using this approximation, a discrete problem is

introduced, where the goal is to find the coefficients, p,, as shown in problem (2.7)

Find p 6 ER" that
minimizes f tudy

" (2.7)
subject to 1.. EL. p.q.(x)dw s unease») x V

Equilibrium

Traditionally, the discrete version of the problem is based on a finite element
discretization of the design domain. The finite element discretization is used for two
purposes: to provide a rigorous approximation of the displacement field u and to
introduce an orthogonal basis, {q,(x)}f:'f (mesh size n), for the approximation to

the structure indicator function, XA(x). The space occupied by each finite element,

e, corresponds to the support of the basis function, qe(x). The “indicator” basis

11

 

 

 

 

 

 

Design

 

Domain

 

 

 

 

 

 

 

\\\\\\\

 

 

 

 

(b) Problem finite ele- (c) Approximate optimal

(a) Problem statement . . .
ment discretization topology

Figure 2.5: Example of forming an optimal design when approximating x(x)
functions are defined as

1 if location x is contained in element 6
qe(X) = (2-8)
0 if location x is not contained in element e
Since each basis function is associated with a single finite element, each coefficient,
pe, is often referred to as the effective density for element e. Using this setting,
the Optimal structure is formed by distributing material in various elements in the
design domain. It is also assumed that the elastic tensor in a given element, De, is
constant over the element and uniquely characterized by effective density pe of the
element. This process is shown graphically in Figure 2.5. We now examine attempts
to solve various forms of the topology optimization problem using this type of discrete

approximation.

12

2.2.1 Macroscopic Topology Optimization

Macroscopic topology optimization is noted as the first attempt to solve an approxi-
mate form of the topology optimization problem stated in (2.2). Again, the objective
is to find the stiffest structure formed by distributing a fixed amount of material in
the design domain subject to the specified boundary conditions. The material is dis-
tributed in binary fashion, designating each element in the design domain as either
full of material or empty. This results in an integer programming problem where
the goal is to determine the unknown variables {p3}: that define the stiffest struc-
ture, i.e., the optimal material distribution in the domain w . Formally, this may be

expressed as

Find p e {0,1}n that

minimizes fTu
(2.9)

SUbjCCt to 22:1 ”8 (pa — V) S 0
K(p)u — f = 0

where

n = number of elements

p = (pump-spa)

pe = density of element e, 1 or 0
126 = element volume

K = stiffness matrix (FEM)

V = prescribed percent volume

13

u = nodal displacement vector (FEM)
f 2 force vector (FEM)

This type of formulation may be interpreted as an approach that attempts to deter-
mine the optimum mixture of two materials: solid (p = 1) and void (p = 0). Typically,
the solid material is associated with an isotrOpic material. Here, the mixture of the
two materials is occuring at the macroscopic level determined by the scale of the fi-
nite element discretization. This type of problem is not well-posed as shown by Kohn
and Strang [18, 19, 20] and Murat and Tartar [26], who showed that the sequence of
designs that minimize compliance does not converge in the space of macroscopically
mixed materials. The problem must be reformulated to include mixtures that occur

at a small, microsc0pic scale, i.e., a composite mixture of solid and void material.

2.2.2 Topology Optimization Using Homogenization

A new approximation to the topology Optimization problem is introduced to overcome
the problems associated with the macroscopic topology Optimization. Rather than
determining the mixture of two materials at the macroscopic level, the mixture is al-
lowed to occur at a small scale, or microscopic level. This leads to problem formulation
using material with a microstructure, i.e., a material with micrOSCOpic perforations
of different sizes controlled by the effective density p, which may vary continuously
from 0 to 1. As in macroscopic topology optimization, the case p = 0 corresponds to
a void material and the case p = 1 corresponds to a solid, isotropic material. In this

formulation, however, the intermediate values, 0 < p < 1, correspond to a composite

14

 

(a) Existance of microscopic material (b) Micro—hole unit cell

Figure 2.6: Two-dimensional microcell in macroscopic materials

material where a mixture of solid and void takes places at the microscopic level. This
relaxation complicates the relation between the effective density p and the tensor D
of effective material properties. This relation is typically determined through the use
of a homogenization method. In a homogenization method, a mixture of material
described at the microscopic level is used to determine effective material properties
at the macroscopic level. Because of the geometric complexity of the microstructure,
multiple variables are needed to describe a given material mixture. For example,
when using the well known micro-hole microscopic description [3] in two dimensions,
a set of variables a, b, and 0 are required to describe the material. For a given cell, a
and b denote cell size parameters while 0 represents the cell orientation, as shown in
Figure 2.6. The effective density and effective material properties at the macroscopic

level are determined based on the cell parameters. When a finite element discretiza-

tion is introduced, it is assumed that a single set of cell parameters, ac, be and 0..
are used to describe the material within each finite element, but may vary from ele-
ment to element. This leads to element-wise constant effective densities and material
properties based on ae, be and 08. In two dimensions, the effective density in a given

element may be related to the element’s cell Size parameters by

pe : (1 — ae)(1 _ be) (210)

The effective material properties are obtained using the theory of homogenization
of materials with a periodic microstructure [3]. Since this analysis may be time
consuming, the effective material properties for varying cell parameters are computed
and stored in tabular form. Subsequent evaluations of the effective material properties
are obtained by interpolating from the tabular values.

Remark: The micro-hole model may be expanded to three dimensions as Shown by
Suzuki and Kikuchi [30, 31]. This leads to three cell variables a, b, and c to describe
the size of the perforation in the microcell and three orientation variables 11), (b, and 0.

Formally, topology Optimization based on a homogenization strategy (in two di-

16

mensions) is expressed in (2.11)

Find 2 E 323" that

minimizes fTu

subject to 2:1 ue(pe — V) S 0
0<a,-Sl fori=1,n (2-11)
0<b551 fori=1,n

—7r/2 < 0,- S 7r/2 for i = 1,n

K(z)u — f = 0

where

n = number of elements

2 = (a1,b1,01, . ..,a,,,bn,0n)

p. = (1 - ae)(1 - be)

ve = element volume

K = stiffness matrix (FEM)

V = prescribed percent volume

11 = nodal displacement vector (FEM)

f = force vector (FEM)

Unlike the macroscopic topology optimization problem, homogenization-based tOpol-
ogy optimization is well-posed.

17

TOpology optimization using homogenization is widely used to solve design prob-
lems, but it does possess a few drawbacks. Since the effective densities are allowed to
vary, the optimal design may contain solutions that have microscale mixtures when-
ever 0 < pe < 1. Such areas are impossible to manufacture and must be replaced by

some reasonable, practical engineering approximation.

2.2.3 Topology Optimization Using A Fictitious Material

Model

Another approximation to the topology Optimization problem statement (2.2) is to
formulate the problem using a fictitious material model [25]. Similar to topology
Optimization using homogenization, the effective densities p may vary between 0 and 1.
However, rather than introducing a microstructure to determine the effective material

properties, the following relationship is used:

fjlcl = pngjkh P>1 (2-12)

where ijk, corresponds to the effective material tensor for element e and Dg-k, corre-
sponds to an isotropic material tensor. Here, intermediate densities are penalized by
the exponent p. The effect of p is illustrated in Figure 2.7. When p = 1, the material
strength and the amount of material used are directly proportional. However, as the
value of p increases, the strength of the material at a given value of p becomes in-
creasingly weak, hence penalizing the intermediate value Of p. As in other approaches

discussed, this approach assumes that the effective material properties are constant

18

 

Material
Strength

 

 

 

 

0 0.2 0.4 0.6 0.8 1
Densityp

Figure 2.7: Effect of material exponent p

in each element.
Formally, topology optimization using a fictitious material model (in two dimen-

sions) is expressed in (2.13).

Find p E 32“ that

minimizes fT u

subject to 22:1 ue(pe — V) S 0 (2-13)
0<p,-_<_l fori=1,n

K(p)u - f = 0
where
n = number of elements

p = (plap21-°°Ipn)

19

v,- = element volume

K = element stiffness matrix
V = prescribed percent volume
11 = nodal displacement vector

f = force vector

Topology Optimization using a fictitious material model is another popular alternative
to solving the original topology Optimization problem. The selection of p is usually
made based on experience and on the type of problem being solved, and it typically
ranges from p = 2 to p = 3. The approach is simple and computationally efficient:
there is only one design variable per element and the relation between p and the
elastic tensor D is straightforward. However, in this case too, unwanted intermediate

densities are Often present in Optimal designs.

2.3 Solution Methods

In all formulations of the topology Optimization problem the number of design vari-
ables becomes increasingly large as the number Of elements used to discretize the
design domain grows. To overcome the difficulties associated with a potentially very
large number of the design variables, the optimization problem is solved using an
iterative solution technique, based on the Lagrangian optimality conditions [3]. Such
algorithms require that the structure’s performance be evaluated at least once per
iteration, i.e., via finite element analysis. But, as shown in Figure 2.8, in the finite

20

 

Design

 

Domain

 

 

 

 

 

\\\\\\\\\

(b) Optimal topology: (c) Optimal topology:

(a) Problem statement coarse appoximation refined appoximation

Figure 2.8: Effect of discretization level on optimal shape resolution

element based approach, the finite element discretization provides a basis for both the
analysis and the shape of the structure. Thus, a detailed resolution of the geometry
is accompanied by a proportional increase in the size of the finite element model.
As the problem size grows, the solution of the finite element equations becomes the
dominating computational factor in the entire solution process, governing the com-
putational effort required to solve the problem. This may make solving the problem
computationally unrealistic. This is especially true in three-dimensional problems

where a fine mesh resolution is required for adequate shape representation.

2.4 Summary

Alternative solution strategies to solve the topology optimization problem were dis-
cussed. In all cases, the analysis of a structure’s performance for a given design is

21

noted as the dominant factor computationally. Finite element methods are typically
used for this analysis; however, as the problem size grows, traditional finite element
methods may require excessive computational resources. In the following chapters we
suggest replacing traditional finite element methods with alternative solution tech-

niques, i.e., meshless methods based on Wavelet-Galerkin techniques.

22

Chapter 3

Wavelet Approximations in

Numerical Analysis

3.1 Introduction

This chapter discusses some basic concepts related to wavelet analysis as it pertains to
this dissertation. By necessity, the discussion is brief. The interested reader is referred
to references [4] and [35]. Wavelets are functions that satisfy Special properties and
are Often used in signal analysis for their ability to capture localized phenomena in
both space and frequency. Since the wavelets are localized in space, they are typically
non-zero only on a short interval (compactly supported). Since wavelets are localized
in frequency, their non-zero portion contains frequency information that varies from
wavelet to wavelet. Because of their localization properties in both space and fre-
quency, a wavelet decomposition of a signal provides information about high and low

frequency components and where they occur. In many applications, this may have

23

definite advantages over standard Fourier analysis, which identifies frequency infor-
mation but no spatial information. In addition, a variety of very efficient algorithms
for wavelets have been developed, making wavelet decompositions competitive even
with very computationally efficient methods such as Fast Fourier 'II‘ansforms. Over
the past few years, wavelets have become a very popular research topic. For example,
wavelets have been applied to general signal processing, image processing, solving
ordinary differential equations, and solving partial differential equations [4].

The wide range of applications of wavelets in analysis has lead to the development
of many different types of wavelets that are Specially tailored to effectively solve the
problem at hand. The different types of wavelets are Often referred to as wavelet
families. Usually, wavelet families are based on two special wavelets: the father
wavelet and the mother wavelet. The father wavelet is called a scaling function and
the mother wavelet is called a wavelet function, or simply a wavelet. The father and
mother wavelets are used to develop a set Of functions to form a wavelet system in
a given wavelet family. Each wavelet system contains functions that vary in space
and frequency to produce bases that may be used to describe data at varying levels
of detail. Two characteristics of father and mother wavelets are traditionally used
to evaluate their ability to capture local phenomena: support and smoothness. The
support measures the length of interval in which a wavelet is non-zero, while the
smoothness measures how well a wavelet expansion can represent polynomials. In
general, the support of a wavelet and its smoothness are closely related. As the
support of the wavelet becomes shorter, the wavelet becomes less smooth. Likewise,

as the support of the wavelet becomes wider, the wavelet becomes smoother. This

24

 

 

     

 

 

 

 

 

 

-0.5
I -1 I
4.5 J J I l l I I 4.5 I I I l I I I
-4 -3 -2 -l 0 1 2 3 4 -4 -3 -2 -l 0 1 2 3 4
x x
(a) Daubechies D6 wavelet (b) Daubechies D12 wavelet

Figure 3.1: Sample Daubechies wavelets

feature is illustrated in Figure 3.1, where two wavelets from the Daubechies family
are displayed. Note that a wavelet from the Daubechies D6 wavelet has much shorter

support and is less smooth than a wavelet from the Daubechies D12 wavelet.

3.2 Multiresolution Analysis

Wavelets are usually introduced in one of two ways: via the continuous wavelet trans-
form or in conjunction with multiresolution analysis [32]. Here, wavelets will be pre-
sented in a multiresolution analysis context. Multiresolution analysis is the ability to
represent functions at varying level of detail. A formal definition of multiresolution
analysis using wavelets was introduced by Daubechies [13] and Mallat [24]. A mul-

tiresolution analysis is achieved by generating hierarchical bases where basis functions

25

consist of scaled translates of a single function, ¢(r), the scaling function. A formal
description Of a multiresolution analysis may be stated as follows:
Given that there exists a function 43(27) 6 V0, such that {¢(a: — k)}lceZ is a Riesz

basis for V0, we wish to find a sequence Of embedded subspaces V,- such that

(P1) 0c---cV_1cVocV1C---cL2(éR)

(P2) fljesz=0

(P3) f($)€Vj®f(2$)€VJ-+1 (3.1)
(P4) f($)€Vo®f($+1)€Vo

(P5) Ujez V]- is dense in L2(§R)

where L2(§R) is the space of all square integrable functions.

The property (P1) in 3.1 indicates a relationship among the embedded subspaces.
Note that increasing the level of resolution, denoted by j, produces subspaces Vj
that more closely approximate L2(§R). Property (P2) states that although spaces
V, are embedded, they are disjoint. Property (P3) and (P4) are often referred to
as the scaling and translation properties, respectively. These properties Show that
the dilations of a given function in V} can be represented in another space Vj: and
translates of a given function in V,- remain in V,. Property (P5) states that increasing
the level or resolution will provide sufficiently rich approximations of functions in
L2(§R).

More restrictive forms of multiresolution analysis may be developed based on the

prOperties of ¢(a:). For example, an orthogonal multiresolution analysis occurs when

26

{¢(2: — k)}kez is an orthonormal basis for V0. Wavelets associated with this type of
multiresolution analysis are called orthogonal wavelets. When {¢(:c — k)}kez is a non-
orthogonal Riesz basis for V0, non-orthogonal wavelets are generated. For the work
presented here, it is assumed that both analyses, orthogonal and non-orthogonal,
are constructed by translates and dilates of a single function ¢(a:). This leads to a
shift-invariant wavelet system constructed with shift-invariant wavelets.

Remark: The term wavelets is often used to identify scaling functions and/or

 

wavelet functions. A wavelet system may contain scaling functions and wavelet func-

tions.

3.3 Defining Wavelets

Shift-invariant wavelets are defined recursively by a scaling equation. The scaling
equation is also called the refinement equation or the dilation equation. The scaling

function is recursively defined by the scaling equation

¢(rr) = f: a.¢(2:c -— k) (3.2)

Ic=-oo

The wavelet function is defined by

112(20) =. i bra/1&2: — k) (3.3)

k=-oo

27

The relationship between the coefficients oh and bk is

bk = (—1)"_lak (3.4)

The scaling equation coefficients ak are also referred to as filter coefficients because
filter Operations are often associated with wavelet analysis. A nice feature of wavelets
is that a finite and usually small number of filter coefficients are needed to define them.
The finite number of filter coefficients used to define a wavelet control its support
and smoothness. The wavelets shown Figure 3.1 demonstrate the effect of varying
the number of non-zero coefficients. The Daubechies D6 wavelet is constructed with
six non-zero coefficients while the Daubechies D12 wavelet is constructed with twelve
non-zero coefficients. In general, a larger number of coefficients leads to longer and

smoother scaling functions and wavelets.

3.4 Defining Wavelet Bases

Wavelets may be used to make up a complete set of basis functions. The basis
functions are defined by translations and dilations of the scaling function and/or
wavelet function. Translations and dilations of the scaling and wavelet functions are
of the following form:

¢,,,(x) = 2j/2¢(2j:c - k) (3.5)

and

1,151,017) = 217%(21'1: — k) (3.6)

28

The first index j, or level, controls the dilation of the function, while the second
index k controls the translation of the function. Figure 3.2 displays a few examples
of translates and dilates for various wavelets. Note how changes in k translate the
function, preserving its shape, while changes in 3' scale the function. Translates and

dilates of the scaling functions and wavelets can be used to define spaces of functions,

V,- and W,-, as follows:

V,- = span {¢,k(x) : k E Z} (3.7)

and
W,- 2 spam {115-fix) : k E Z} (3.8)

Note that

V,- C L2(§R) and W,- C L2(§R) (3.9)

and
L2(§R) = V0 EB éWj] (3.10)

j=0

 

The wavelet spaces V,- and W,- also posses another property based on the embedded

structure

V1+1 = V,- EB WJ’ (3-11)

This relationship helps explain how the multiresolution analysis works. If a function is
approximated in the space V,- and a finer approximation is desired, it may be obtained
by adding to the approximation a function in the space W,-. Since scaling functions

are used as the basis functions for the subspaces V, and wavelets are used as the

29

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2.5 I I I I I I I 2 5 I I I I I I l
2 - - 2 - -
| I
1.5 '- | - 1.5 P l -
1 - - 1 - -
| |
0.5 - - 0.5 - d
o :. a o : A .

'05 _ -‘ -05 I- V —I
-1 I— .4 -1 I— I -I
-1.5 h I .. -1 5 r I _
_2 l L I L I I I _2 l I I I I L I

-2 0 2 4 6 8 10 12 -2 0 2 4 6 8 10 12
x x
(a) Daubechies D6, Ibo,o(a:) (b) Daubechies D6, ¢o,Io($)
2-5 l I I I I I I 2 5 l I r r I I I
2 - - 2 -‘
| l
1.5 - l - 1.5 l -
l - l - 1 F | -
0.5 - d 05 - -
l I
o h o A l——
-0.5 - -* -05 - I -
-1 I- l -I -1 I- I .-
-1.5 - | -+ -1.5 l- l n
_2 l L I I I I I _2 l I L I I I L
-2 0 2 4 6 8 10 12 -2 0 2 4 6 8 10 12
x x
(c) Daubechies D6, $1,-2(a:) (d) Daubechies D6, $1,13(a:)

Figure 3.2: Effects of dilation and translation of wavelets

30

basis functions for the subspaces W,-, it follows that wavelets approximations provide

a method of easily increasing the resolution of a given function approximation.

3.5 Examples of Orthogonal Wavelets

This section presents examples of two orthogonal wavelets: the Haar wavelet and
the Daubechies wavelet. These wavelets are discussed because they will be used

extensively in upcoming chapters.

3.5.1 Haar Wavelets

Haar wavelets are the oldest and simplest orthogonal wavelets. Haar wavelets have
scaling equation coefficients a = {1, 1}. The Haar father wavelet and mother wavelet

are defined as

¢haar (113) = X[0,1) (2C) (3.12)

and

1P”“"'@) = X[0,1/2)($) — X[1/2,1)($) (3.13)

where x(a:) is an indicator function. An example of both the compactly supported
Haar scaling function and wavelet are Shown in Figure 3.3. For a given level, j,
translates of Haar scaling functions are orthogonal and translates of Haar wavelet

functions are orthogonal. Also, any Haar scaling function (i.e., any level or transla-

31

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

j l I T _I I I r
1 - I - 1 - ] J
I
05 - a 0.5 L l -
|
o — — o — —
|
-0.5 - I -' -0.5 '- -
1 | l
l J l J l I l l l I II]
-1 -0.5 0 0.5 1 1.5 2 -1 -0.5 0 0.5 1 1.5 2
X X
(a) Haar scaling function (b) Haar wavelet function
Figure 3.3: Example Haar scaling function and wavelet
tion) is orthogonal to any Haar wavelet function, i.e.
(A1) 13°00 ¢jk($)¢fl($)d$ = 6k_1 j, 16,1 6 Z
(A2) 13°... anon-ave = 6H j, w e z (3-14)

(A3) If?» ¢jk($)1/)jrkl($)d$ = 0 j, k,j’, k, E Z

Note that properties (A1) and (A2) are guaranteed by the compact support of the
Haar scaling functions and wavelet functions. However, property (A3), a much more
difficult condition to satisfy, is guaranteed from the relationship between ¢”““'(a:),
1b’““"(:r) and the translation and dilation formulae (3.5) and (3.6).

The Haar scaling function and wavelet are both discontinuous and therefore may
be applied in a limited number of analyses. Next, we examine an orthogonal wavelet

that is continuous and has continuous derivatives, the Daubechies D6 wavelet.

32

 

 

 

 

 

 

 

 

 

(a) Daubechies D6 scaling function (b) Daubechies D6 wavelet function

Figure 3.4: Example Daubechies D6 scaling function and wavelet

3.5.2 Daubechies D6 Wavelets

Modern use of wavelets is largely due to the work of Daubechies [13]. In this work,
wavelets were constructed based on the orthogonal multiresolution analysis properties.
Daubechies created a wavelet family with infinite wavelet systems characterized by
the number of non-zero coefficients in the scaling equation. Of interest here are the
Daubechies wavelets constructed with six non-zero coefficients (n = 6). Examples of
the compactly supported Daubechies D6 wavelets are shown in Figure 3.4.

A special property of the Daubechies wavelets is that they satisfy a vanishing

moment condition based on the number of filter coefficients used to define the wavelet.

33

Namely, Daubechies wavelets with n filter coefficients satisfy
/ x‘¢(x)dx = 0, l = 0,1,...,n/2 (3.15)

Note that all Daubechies scaling and wavelet functions satisfy the properties listed
in (3.14). Also note that unlike the Haar wavelet system, there is no explicit expression
for the Daubechies D6 scaling functions or wavelet functions. This is true for all
Daubechies wavelets except for the case when n = 2, where the Daubechies wavelet

system reduces to the Haar wavelet system.

3.6 Examples of Non-Orthogonal Wavelets

Non-Orthogonal wavelets are generated in conjunction with a non-orthogonal mul-
tiresolution analysis. An example of non-orthogonal wavelets are spline wavelets
which are based on Cardinal B-Splines [9] and [33]. Bosplines of order m may be

expressed as

Nm($) = :pkNm(2-T _ k) (3'16)
k=0
where
m
pk = 2—m-I-l (3.17)
It

Like Daubechies wavelets, the spline wavelet family as an infinite number of spline
wavelet systems based on the number of non-zero filter coefficients in the scaling

equation. Unlike Daubechies wavelets, one may use (3.16) and (3.17) to develop

34

explicit expressions for spline father functions and mother functions of order m

Me) = Nmm = ipuvmrzx — k) (3.18)
k=0
and
3111—2
212“") (x) = Z qum(2x — k) (3.19)
k=0
where,
m m
qr = (-1)"2""‘ Z N2m(k + 1 - l) (3.20)
l=0 l

Two spline wavelets will be discussed here: the linear spline wavelets (2nd order

B-Splines) and the quadratic spline wavelets (3rd order B-Spline).

3.6.1 Linear Spline Wavelets

Linear spline wavelets, denoted as S3 wavelets, are based on linear B-splines (order
m = 2) and are constructed with the use Of three, (m + 1), non-zero filter coefficients

a = {%, 1, -;—}. The father and mother wavelets for linear spline wavelets are defined

as
a: 0 _<_ a: < 1
¢53($) = N2($l = 2 — a: 1 S a: < 2 (3-21)
0 otherwise
and

1933(2) = 11—2N2(23:) — §N2(23: — 1) + %N2(255 “ 2)
(3.22)

—%N2(2a; — 3) + fliN2(2a: — 4)

35

 

 

0.8
0.6
0.4
0.2

 

‘0.2 - —(
-0.4 _ -l

 

 

 

 

 

 

'0.6 l l I l l

 

(a) Linear spline scaling function (b) Linear spline wavelet function

Figure 3.5: Example linear spline scaling function and wavelet

An example of both the compactly supported linear spline scaling function and
wavelet are shown in Figure 3.5. Unlike Daubechies D6 scaling functions, linear
Spline scaling functions satisfy (15(2) 2 0 Va: E 33. This behavior coupled with the
support length restricts the possibility of using translations of 42(23) as orthogonal basis
functions. Also note that linear spline scaling functions and wavelets are symmetric
about the midpoint of their support. In many analysis problems, the symmetry of

the linear spline wavelets help in producing more accurate results.

3.6.2 Quadratic Spline Wavelets

Quadratic spline wavelets, denoted S4 wavelets, are based on quadratic B-splines
(order m = 3) and are constructed with the use of four, (m + 1), non-zero filter

coefficients a = {31, %, g, i}. The father and mother wavelets for quadratic spline

36

wavelets are defined as

 

$172 0 S a: < 1

S4 i-(w-%)2 IS$<2
‘1’ (3’) = N3“) = i (3.23)

%($ — 3)2 2 S a: < 3

0 otherwise

L
and

$54”) = TsoN3(2$)- r2890N3(2:z: — 1) + 1—3—3N3(2x — 2)

—:——2:Na<2x - 3>+ :—::N3<2x — 4)— :—;gzv.<2x - s) (3.24)
+4280—2N3(2$— 6)— 4—80—1_N3(2$_ 7)

An example of both the quadratic Spline scaling function and wavelet are shown in Fig-
ure 3.6. Comparing the linear and quadratic spline scaling functions and wavelets, one
notices the wider support and smoother behavior for the quadratic system. Similar to
the linear spline wavelets, quadratic spline scaling functions satisfy (15(3) 2 0 Va: 6 3?

and are symmetric about the midpoint of their support.

3.7 Wavelet Representations and Algorithms

This section briefly discusses how wavelets may be used to approximate functions and
algorithms associated with these approximations. Up to now, all wavelets and wavelet
bases were defined along the real line and infinite series were used for approximat-

ing functions. In practical applications, however, one usually investigates bounded

37

 

 

 

 

 

 

 

 

 

 

(a) Quadratic spline scaling function (b) Quadratic spline wavelet function

Figure 3.6: Example quadratic Spline scaling function and wavelet

domains with finite approximations. Typically, the number of terms in discrete ap-
proximations using wavelets is governed by the level of analysis, j. The number
of wavelets used in this approach (both scaling functions and wavelet functions) is
N = 23'. It is also assumed that the function is periodic. For some situations, i.e.,
image compression and signal analysis, the periodicity assumption does not interfere
with the analysis. However, in other situations, i.e., solution of ordinary and partial
differential equations, the periodicity may cause some difficulties. For now, we will
examine how wavelets are used to approximate periodic functions and discuss the
associated algorithms.

Assume that the function f (x) E L2(§R) is periodic with period L. There are two
methods for approximating the function f (2:) using wavelets: fixed-scaled expansions

and multi-scale expansions. In fixed-scaled expansions, a function is approximated

38

using translations of scaling function at a fixed level, J, i.e.

2J—1

f($) z 2 CJk¢Jk($) (3.25)
k=0

When orthogonal scaling functions are used in the approximation, the coefficients CJk

are found by

CJk = All f($)¢1k($)d$ (3.26)

When non-orthogonal scaling functions are used in the approximation, the coefficients

ck are found by solving the system of equations:

2J—1

L L
2 / CJk¢Jk($)¢Ji($)d$ = / f($)¢Ji($)d-’B for i = 0, 1, - - -,2J - 1 (3-27)
k=0 0 0

Although the system of equations (3.27) seem very tedious and large for modest values
of J, it may be solved very efficiently using Fourier expansions and the computation-
ally efficient Fourier transforms. In multi-scale expansions, a function is approximated

using a base scaling function and translations of wavelets at various levels, i.e.

J 21—1

“33) z Coo¢oo($) + Z Z djkI/ij($) (3.28)

j=0 k=0

When orthogonal scaling functions and wavelets are used in the approximation, the

39

coefficients 000 and djk are found by

000 = foLf($)¢oo($ld$
djk = foLf(x)1,l),-k(x)dx

(3.29)

When non-orthogonal scaling functions are used in the approximation, the coefficients

Coo and dfl, are found by solving the following system of equations

123-1

/0 when“: :3 / d,.w,-.(x) )¢oo(x)dx = /0 f(z)¢oo(z)dx (3.30)

j=0k= 0
and f0rj'=0,1,,_,,Jand k’=0,1,...,2j’-1

J2j—l

[L COO¢00 ($)¢jlkl($) )dx + Z 2 _/0L djch/ijcv (low-ilk! ($)d$ = ‘/(;L f($)’1pjlkl($)dx
j=0 k=0
(3.31)
Note that multi-scale approximations of a given function may contain fixed-scaled

approximations at higher levels than shown. For example, a common multi-scale

approximation has the form

21‘1—1 2i— 1

x) z : cj1k¢j1k($ )+ z Z djk¢jk(x) (332)
k=0

J=Jil== 0

The fixed-scale coefficients, {eykfijj‘l are often called the scaling function co-
efficients for level J. Similarly, the portion of the multi-scale coefficients {d,-k},,:2J 1
is called the wavelet coefficients for level j. It is also typical to identify a fixed-

Scale approximation of a given function simply by its coefficients {ch}::3J‘1. This is

40

also true for multi-scale with approximates in the form {c,-,k}§:(2,j’_1 and {d,-qc : j’ =

j1,J— 1, k =0,2J" — 1}

In most wavelet applications, a high level fixed-scale approximation is used for the
initial representation for the data.‘ Then, a multi-scale representation is generated
using a very efficient wavelet algorithm, the forward wavelet transform. Typically,
the multi-scale representation is then analyzed and possibly modified. Whenever
it is necessary, the modified multi-scale representation may be converted back to a
higher level fixed-scale approximation using the inverse wavelet transform. Both the
forward and inverse wavelet transform are applied, in stages, one level at a time.
For example, given a level j fixed-scale representation of a function, {cjk}’,::gj"l,
application of one stage of the forward wavelet transform would produce the scaling

function coefficients for level j — 1, {c(,_1)k}fi:gj_l‘l, and the wavelet coefficients for

level 3' — 1, {dU_1)k}’,::§j-l“l. A similar pattern occurs when applying the inverse
wavelet transform. However, in the inverse transform application, the level j fixed-

scale coefficients are generated from the level j - 1 scaling function and wavelet

coefficients.

3.8 Wavelets In Multiple Dimensions

Up to now, wavelets have been discussed for one-dimensional examples. For the prob—
lems we wish address here, wavelets in multiple dimensions are required. Here, only
fixed-scale expansions using scaling functions in multiple dimensions are discussed.

For a discussion of wavelets in multiple dimensions see Daubechies [13]. One simple

41

method for extending wavelets to multi-dimensions is with the use of tensor products.

A scaling function in two-dimensions may be defined as

4531:1617, ll) = ¢jk($)¢jr(y) (3-33)

A level j fixed-scaled approximation requires an array of 23' x 21' scaling function co-
efficients. This type of representation is often compared to an image with dimensions
2i x 23' in pixels. In the case of Haar fixed-scale expansions, the support of each scal-
ing function exactly represents the area occupied by a pixel. Hence, the fixed-scale
coefficients may directly related to corresponding pixel intensity.

Three-dimensional scaling functions may be expressed in a similar fashion, namely

¢jkzm(x, y, Z) = ¢jk($)¢jt(y)¢jm(2) (334)

Of special interest are Haar fixed-scale expansions. Here, the support of each scaling
function exactly represents the volume occupied by a voxel, the three-dimensional
equivalent of pixel. Similar to pixels, voxels can be used to represent three-dimensional
images and the fixed-scale coefficients may directly related to corresponding voxel

intensity.

42

3.9 Wavelet Applications

3.9.1 Function approximations

We first examine how various wavelets are used to approximate functions. Figure 3.7
shows how the Haar basis may be used to approximate a function at varying levels
Of detail. The multiresolution analysis is performed in the following steps. A level
3' = 8 fixed-scaled representation of the original functions is obtained. This requires
28 = 256 data points that represent the coefficients of 256 scaling functions. Then the
forward wavelet transform is applied to develop representations at lower levels. The
lower level representations shown in Figure 3.7 are level j = 2, 22 = 4 data points;
level j = 4, 24 = 16 data points; level j = 6, 26 = 64 data points. Note that as
the level of approximation increases, the approximations become more accurate. In
many cases, the Haar basis is not sufficiently smooth to obtain reasonable approxi-
mations. We now examine the same multiresolution analysis when using Daubechies
D6 wavelets. As shown in Figure 3.8, we see that the multiresolution analysis of the
given function using Daubechies D6 wavelets provides close approximations at rela-
tively low level. Unlike the Haar basis, the Daubechies based approximation achieves

near exact representation at level 6, requiring only 26 == 64 coefficients.

3.9.2 Image Compression

Wavelets may also be used for image compression. Here a two-dimensional image
of a mechanical component is photographed using a digital camera. The original
image is Shown in Figure 3.9. This image was then analyzed using two-dimensional

43

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.1 I I l I 0.06 | 1 I l
0.04 - -
0.05
0.02 - -
0 0 I—- — —— — -— — — _—
-0.02 L -
‘0.05 _004 ,_ _,
-0.06 - ..
-0.1
-0.08 - -
-o.15 -0.1 ' ' ' '
0 0.2 0.4 0.6 0.8 1
x
(b) Level 2 representation
0-1 I I I I 0.1 I I I I
0.05 0.05
0 0
-0.05 -0.05
-0.1 -0.1
-0.15 -0.15
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 l
x x
(O) Level 4 representation ((1) Level 6 representation

Figure 3.7: Multiresolution using haar wavelets

44

 

 

0-1 I I I I 0.1 I I I I

0.05

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0,1 , , , , 0.1 . I . I
0.05 0.05
0 0
-0.05 -0.05
-0.1 -0.1
-0.15 -0.15

0 0.2 0.4 0.6 0.8 l 0 0.2 0.4 0.6 0.8 l
x x
(c) Level 4 representation (d) Level 6 representation

Figure 3.8: Multiresolution using Daubechies D6 wavelets

45

 

Figure 3.9: Original component image

wavelet Daubechies D4 wavelets (Daubechies wavelet constructed with four non-zero
coefficients). A full two-dimensional multi-scale approximation was found. Rom
this representation, only the largest (absolute value) 5% of the scaling function and
wavelet coefficients were stored. The inverse wavelet transform was performed on
the compressed data. The resulting image is shown in Figure 3.10. Comparing
the reconstructed figure to the original figure shows the power of the multi-scale
representations. In this case, only a small fraction of the wavelet coefficients were
used to construct the image and yet differences between the two images are difficult
to detect. This example shows the amazing data compression that may be obtained

using wavelet based approximations.

46

 

Figure 3.10: Component image created from 5% of the wavelet coefficients

3.9.3 Wavelet-Based PDE Solvers

Solving partial differential equations with a Wavelet-Galerkin scheme has received
considerable attention in recent years. The main idea of the Wavelet-Galerkin ap-
proach for solving partial differential equations is to use wavelet subspaces, i.e., fixed
scale spaces V,- or multi-scale spaces Vo 699:0 W,-:, as the finite-dimensional subspaces
in traditional Galerkin methods. Most efforts in this area fall under one of two cate-

gories based on the type of approximation spaces used in the Galerkin technique:
1. Refinable shift-invariant wavelet approximation spaces
2. Explicit geometric wavelet approximation spaces

The two approaches arise primarily from using different strategies for satisfying

boundary conditions and approximating the domain of the problem. The graph in

47

 

Satisfying Boundary Conditions with

 

 

 

 

 

Wavelet-Galerkin Methods
Refinable Refinable
Geometric Shift-Invariant
Wavelet Spaces Wavelet Spaces

 

 

 

 

 

 

 

 

 

 

1D Operator n-D Biorthogonal Wavelet Element Capacitance Penalty Lagrange
Wavelets Wavelets Method Method Method Multipliers

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3.11: Wavelet based approaches for solving ordinary and partial differential
equations

Figure 3.11 illustrates how different approaches stem off the type of wavelets used
to satisfy the boundary conditions. Each block in the bottom row of Figure 3.11
corresponds to a method for solving ordinary and / or partial differential equations. A

description of the various methods and references to appropriate literature follows.

Refinable Geometric Wavelet Spaces

This class of Wavelet-Galerkin solvers is based on the development of wavelets that
satisfy boundary conditions by construction. Given a wavelet basis where each ba-
sis function satisfies the homogeneous form of the boundary conditions, all functions
constructed with that basis must also satisfy the homogeneous form of the boundary

conditions. First attempts of this type were used to solve ordinary differential equa-

48

tions with the use of 1D operator wavelets [36] and [32]. Here, special biorthogonal
wavelets were integrated to create operator wavelets leading to an operator wavelet
basis. In some cases, the matrices associated with the discrete system of equations
are diagonal when represented in the operator wavelet basis. These ideas quickly
progressed into multiple dimensions to solve partial differential equations. Although
these methods did not produce diagonal system matrices, the condition number of
these matrices with appropriate preconditioning are independent of the level of re-
finement. Examples of this type of approach are wavelets based on stationary sub-
division [10], composite wavelet bases [12], and divergence free wavelets to solve the
Stokes equations [34]. Recently, there have been developments in the interaction of
biorthogonal wavelets and Spectral Element Methods [6] and [7]. This work appears
very promising offering solutions of partial differential equations with fairly general

domains / boundaries.

Refinable Shift-Invariant Wavelet Spaces

This class of Wavelet-Galerkin solvers is based on developing efficient solution strate-
gies by exploiting the structure of shift-invariant wavelet spaces. Some methods in
this area were based on solving periodic problems and applying boundary conditions
with the‘use of a capacitance matrix method [27] and [17]. These approaches cen-
tered on the special structure of the system matrix for periodic problems. Other
approaches focused on the use of fictitious domains where the application of bound-
ary conditions is accomplished with the use of penalty methods [28], [29] and [16].

These methods lead to efficient solution of problems with arbitrary domains with

49

the use of multi-scale preconditioning. Lastly, one approach incorporates boundary
conditions by appending them by means of Lagrange multipliers [21]. This type of
problem formulation leads to saddle point problems where a class of multi-scale pre-
conditioners provides condition numbers of the system matrix that are independent

of resolution.

3.10 Summary

Examples of Haar, Daubechies, and Spline wavelets were introduced and defined.
The properties of multiresolution analysis were presented and the use of wavelets as
basis functions was discussed. Various applications that use wavelet techniques were
reviewed with a special interest in methods for solving partial differential equations.
The next chapter discusses a new meshless method to solve the associated elasticity
equation arising in a topology optimization setting. The new method heavily relies

on the use of wavelet-based techniques.

50

Chapter 4

Fictitious DOmain Approach Using
Wavelet-Galerkin Techniques For
Solving Partial Differential

Equations

4.1 Background

This chapter introduces a new method for solving partial differential equations using
meshless Wavelet-Galerkin techniques. Here, we wish to find approximate solutions
of elasticity problems associated with large-scale topology optimization problems.
Traditionally, these problems are solved using a finite element method. However,

for large-scale problems, finite element methods possess two drawbacks. First, large—

51

scale problems typically require tedious meshing of the domain of interest. Second,
determining solutions for large-scale problems may be too computationally expensive.
The purpose Of this chapter is to introduce a wavelet-based solution scheme that does
not possess these drawbacks. Here, we suggest the use of a fictitious domain approach

with wavelet-based subspaces used in a Galerkin solution scheme.

4.2 Solving Two-Dimensional Elasticity Problems
Using a Fictitious Domain Approach and

Wavelet Bases

This section discusses a new approach for solving elasticity problems using a fictitious
domain approach and wavelet bases. The presentation focuses on two—dimensional
elasticity problems, but the ideas may be extended to three dimensions. The elas-
ticity problems considered are those that typically arise in the solution of topology
Optimization problems. The approach uses a fictitious domain to recast a given elas-
ticity problem with an arbitrary domain into a periodic problem with a simple, reg-
ular domain. Then, a discrete version of the periodic problem is developed using
a Wavelet-Galerkin technique. The discrete problem is then solved using an effi-
cient preconditioned conjugate gradient solver. The following sections discuss this

approach.

52

 

 

 

V

 

\\\\\\\\\\\\\\\\\

 

 

Figure 4.1: An example topology optimization problem statement

4.2.1 Problem Statement

The topology optimization problem was discussed in Chapter 2. In topology opti-
mization, the goal is to determine the optimal material distribution within a design
domain as given a fixed amount Of material. The optimal structure is determined
from the material distribution, characterized by the design variable 2. The effective
density and effective material at a point (at, y) in to depend on this design variable,
i.e., p : p(z;a:, y) and D : D(z;x, y). Note that the notations (ray) and (231,252) are
used interchanagably.

Figure 4.1 illustrates a typical topology Optimization problem. In the figure,
to represents the design domain or package space, I‘u C 6w is the section of the
boundary where the displacement constraints are applied and I“ C 3w is the section

of the boundary where the loads (tractions) are applied. Note that I’“ U I“ = Ba)

53

and I‘“ r) I“ = (0. The elasticity problem associated with the current design 2 may be

expressed as

Find uw E Kw that
(4.1)
minimizes II w(uw)— — %/ Dwe( (uw) ()uwdw —./1" tude‘
where
Kw = {um E H1(w) : uw, = 0 on I‘" C 6w} (4.2)

is the set of all kinematically admissible solutions. Note that uw = (uw,,uw,) and
t = (t1, t2). In (4.1), Dw(a:, y) is the tensor of elastic properties.

To solve problem (4.1) using a fictitious domain approach, we first embed the
package Space w inside a simplified domain 52 = [0, d] x [0, d], as shown in Figure 4.2.
The domain I), called the fictitious domain, acts as the domain for a new Q—periodic

elasticity problem stated as

Find u 6 Kg that
(4.3)
minimizes Hn(u =—/n Dne( (u)e( u)df2— / fudfl
where
Kn = {11; 6 V9 I Bu,- = 0 on F“ g 6w} (4.4)
and
V9 = {u,- E H1(D) : u,- is Q-periodic} (4.5)

where B is a linear operator on V9, u = (u1,u2), f = (f1,f2) and Dn(a:,y) is the

54

 

 

 

 

 

 

 

é __>

? I“—>

; L-—>
u/ ;
I‘/

/

/

/ a)

/

/

/

/

/

/ n

 

 

 

Figure 4.2: Package space embedded into a fictitious domain

material tensor over the entire fictitious domain (2. Because of the periodic setting of
problem (4.3), the material tensor D9 and forcing function f must be Q—periodic. To
recast the original elasticity problem (4.1) into the fictitious domain problem (4.3),

one must:

1. Define an appropriate material tensor D9 over the entire domain {2

2. Approximate the applied tractions t,- by body forces f;

3. Define Kn in a way such that restriction of u to w satisfies the original boundary

conditions on uw

With reasonable choices of the material tensor Dn, the body force f, and the con-
straint set K9, the solution of the fictitious domain problem (4.3) approaches the

solution of the original problem (4.1) in the domain of interest, i.e., u approaches uw

55

in w. For a more complete discussion of the fictitious domain approach, see the work
of Glowinski et a1 [29] and Bakhvalov and Knyzev [2].
In our approach, the choices for material tensor D9, the body force f, and the

constraint set K Q are as follows:

1. Choosing the material tensor D9. The material properties within the original

 

package space are assigned values of the material tensor 0,, for the current
design z. In the added fictitious domain (2 \ w, the material properties are set
to those associated with a very weak material when compared to the material

in w, where Dweak < D“. We define D9 as

Dweak z, y E Q \ w
139(2): 3’) = (46)

Dw($ay) $.21 E w

2. Approximating tractions with body forces. A body force f defined over 9 is

 

used to approximate the tractions applied on I“. We choose a body force f

such that the work of f in 9 equals the work of the traction t on I“, i.e.

fnfudn = Attuwdw (4.7)

The selection of f is closely coupled with the discretization used to approximate
the domains w and D, so a method for determining f will be discussed in

upcoming sections that address discretization.

3. Constraining the periodic displacement field. To incorporate the essential

 

56

boundary conditions into the fictitious domain problem, we use a Lagrange

multiplier approach. The Lagrangian functional associated with (4.3) - (4.5) is
l(u,)I) = II(u) + [F ATudI‘ (4.8)

where A is the Lagrange multiplier defined on the boundary 1‘". T is an operator
that maps u onto 8w, and is needed because we must map the two-dimensional
displacement field u onto the one-dimensional boundary of w. Stationary con-

ditions of (4.8) lead to an alternative form of the Fictitious Domain Problem

[ane(u)e(v)dQ — [fifvdfl + fr“ AawTvdI‘ = 0

pawTUdF = 0

(4.9)

[‘15

for all v 6 V9 and arbitrary functions A3,, and paw defined on I‘“.

The Fictitious Domain Problem has some advantages over the original elasticity
problem due to the simplified geometry of the problem domain. A typical method
for finding approximate solutions to either of these problems is by using a Galerkin
scheme. In the original elasticity problem, the basis functions are defined over w (a
possibly complex domain); however, in the fictitious domain approach, the periodic
basis functions are defined over the simplified domain (2. Using the simplified periodic
domain to our advantage, we choose shift-invariant wavelet bases to describe the

Galerkin subspaces.

57

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Q

 

Figure 4.3: Level 4 discretization of an example topology Optimization problem

4.2.2 Wavelet Discretization

Fixed-scale wavelet approximations at level j are introduced to solve the Fictitious
Domain Problem (4.9). To aide in the development of the approximations for the ma-
terial tensor, displacement field, and boundary integrals, we introduce an alternative
description of the domain S2. The domain I) is first treated as a square of dimensions
Nx N where N = 2". We then designate each unit square X(k, l) = [k, k+1) x [l, 1+1)
in D as a pixel. This effectively yields a level j, N x N raster description of the do-
main (2 as shown in Figure 4.3. We now introduce the necessary approximations to

develop a discretized form of the Fictitious Domain Problem.

58

Material tensor discretization

The material tensor in two-dimensional elasticity may be expressed in the following

form:

r -

d“(x,y) (11205.31) d’3(:v,y)

D(x,y) = d12(a:,y) d”(a:,y) d23(x,y)

 

(113(3), 31) 612303.11) 6133(27, y)
E . . (4.10)

3111105,?!) E1112(:r,y) 13111203,?!)

" 13112203,?!) E2222(:c,y) E2212(:c,y)

 

 

_ E1112(ar,y) 13221201132!) 19121201831).

The discretization of the material tensor is performed with following wavelet expan-

sion:
21' —1 2i-I

dmn($,y) = Z Z dfifafidx, y) (4.11)

k=0 1:0

where fl? are the expansion coefficients for material tensor entry d'”” and at?“ are
Haar scaling functions at level j. In two dimensions, a fixed-scale Haar wavelet rep-
resentation of an Object is equivalent to an image of the object: each coefficient in
the Haar expansion is related to an associated pixel intensity in the image. There-
fore, the expansion of the material tensor may be directly related to an image of the
current material distribution, i.e., an image of the domain D. Figure 4.4 illustrates
how a shape is more accurately represented by increasing the level of Haar approxi-
mation. When using Haar expansions to represent the material distribution, a high

level resolution may be needed for detailed shape description.

59

 

(a) Level 5 (b) Level 6
(c) Level 7 (d) Level 8

Figure 4.4: Example of using Haar approximations of a circle at various levels

60

Displacement field discretization

The discretization of the displacement field is of the following form:

21-12i-1

“(55:31) = Z Z ujkl¢ykz($,y) (4.12)

k=0 (:0

where u,“ are the expansion coefficients and <15ka are scaling functions that possess
the necessary smoothness. Possible choices for 3‘“ are the Daubechies D6 scaling
function, the quadratic Spline S4 scaling function, or the linear spline S3 scaling
function; all Of which were discussed in Chapter 3. The weight functions v(a:, y) are

approximated in a similar fashion using the expansion

21—121—1

V($,y) = k: 2 ijl¢ykz($ay) (4.13)

:0 1:0

Boundary integral discretization

Here we discuss a method of approximating the boundary integral along 6w with an
integral over the domain 9. From measure theory [15], the integral of any smooth
function 433,, along the boundary 8w may be approximated by an integral involving

a suitable extension of 4530,, ¢(a:, y), and the boundary measure Haw“, i.e.

(a... «an = /n ¢(x, ymawudn (4.14)

61

The boundary measure ||6w|| is expressed as

Haw“ = -wa ° 11 (4-15)

where xw is the indicator function for w and n is the outward normal vector along
6w. In the Fictitious Domain Problem, the function (36,, along 8w is extended to a Q-
periodic function ¢(.z, y). AS in the material discretization, we use a raster description

of the domain to approximate w.

(I) = U(k,()€1wX(k, l) (4.16)

X (k, l) represents the pixel located at (k, l) in the image and 1,, is a set of all pixels
contained in the domain w. Given this representation of w, a discrete pixel-based
representation of the boundary is created and denoted by Aw. The approximate
boundary is formed by collecting the edge pixels of w. An example of such an Oper-
ation is illustrated in Figure 4.5, where a one pixel wide discretized boundary Aw is
shown. Using a discretized representation of the boundary Aw, one could approxi-
mate the boundary measure by

paw if (x, y) 6 Aw
||3w||(x,y) *3 (4.17)

0 otherwise

where the constant paw is adjusted to guarantee equivalent perimeter measurements

62

 

Figure 4.5: Example of edge pixel formation of Aw

of 6w

fawn“: [a ||6deQ (4.18)

Using (4.14), the boundary integral may now be approximated in the following fash-
ion:

/,, an up... Amado (4.19)

Boundary integrals over I‘" and I” may now be approximated by integrals over the

domain 0. For example

fr“ (158de z pp. [AN ¢(:I:,y)df2 (4.20)
and
A, 453de 2 pp. [Art ¢(z, y)dfl (4.21)

where constants pp“ and pr! are set so Pa; = pru + 171‘! and Aw = AF“ U Al“

The above boundary approximations are based on the assumption that the do-

63

main w is of finite perimeter and that the boundary 8w is piecewise smooth. In
most topology Optimization problems, the package space has a simple form so these
approximations are sufficient. If more complex design domains are investigated, the
boundary measure ||8w|| may be approximated using wavelet expansions. See Wells
and Zhou [16] and Glowinski et al [29] for details.

The integrals over the boundary I“ involving traction t are approximated by in-
tegrals over the domain 9 involving a body force f. The body force is represented

using the pixel-wise constant basis expansion

'1‘ , A1“
f(a:,y)= f“ l (x ”E (4.22)

0 otherwise

where f“ is associated with the value of f (3:, y) over the pixel X (k, l). The values
f“ must be adjusted to provide a suitable approximation for the applied traction t.

This is done by forcing f and t to be statically equivalent, meaning

frttude‘z [0 mm (4.23)

We may express these integrals as

t de‘ = f t de‘ 4.24

At u ‘ k2, rtnX(k,z) u ( )
and

dc: do 425

in 21.04» < >

64

Assuming that t is piecewise constant over each pixel boundary, we find that

w _ = 4.2
2: (tklA-‘tnX(k,l)u dF fkl./X(k,l) udfl) 0 ( 6)

(N)

This means that for large-scale problems, the approximation f“ = tkl is sufficient for
representing the applied traction as a body force along the boundary.
It is convenient to express the pixel-wise constant representation of the body force

as a Haar basis expansion of the form

21-12i—1

“3:31) = Z 2: fjkl glgrr($’y) (427)

k=0 (:0

where Haar expansion coefficients fjkl are determined from f“.

Discretized form of the Fictitious Domain Problem

We now use the previously discussed discretizations to develop a discretized form of
the Fictitious Domain Problem stated in (4.9). This is achieved by first converting
all the boundary integrals in (4.9) to domain integrals. For example, through (4.21),

we use the approximation

AawTvdI‘ 8 pr. Avdfl (4.28)
1‘“ AF“

65

where A is the Q-periodic extension of A3,). A similar conversion is made with the

constraint integral, namely

/ pawTudI‘zppu / )2de (4.29)
F" AF“

where p is the Q-periodic extension of How-

Remark: Since we are only interested in /\ along the discretized version of the
boundary Aw, it is convenient to approximate A with a level j fixed-scale Haar ex-
pansion. In this expansion, only the terms corresponding to Haar functions along the
discretized boundary Aw are potentially non-zero.

The approximation of boundary integrals by domain integrals leads to the follow-
ing version of the Fictitious Domain Problem:

[Dn€(u)€(v)dQ-/fvdfl+ppu/ Ade = O for allv 6 V9
9 9 AT“ (4.30)

pl'w AI‘ pudfl = 0 for all p E 62’“

where Q1, is the space of pixel-wise constant functions along the discretized boundary

Aw.

4.2.3 Constructing the Discrete System of Equations

With the introduction of the approximations for the material distribution, displace-

ment functions, boundary conditions and loading conditions into the Fictitious Do-

66

main Problem 4.30), we obtain the discretized problem

A BT U F
= (4.31)
B 0 A 0
where, in two dimensions, A, U, A and F are of the form
A11 A12 111 /\1 f1
A = , U = /\ = and F = (4.32)
A12 A22 112 My f2

The vectors 11,- contain the wavelet expansion coefficients of the displacement field,
and the vectors A,- contain the wavelet expansion coefficients (again only along the
boundary Aw) of the multipliers. Next, we discuss a method for constructing and

solving this system of equations.

Constructing stiffness matrices associated with a heterogeneous material

distribution

Here we examine the construction of the stiffness matrix A when a heterogeneous
material distribution exists in Q. Expressions for entries in the stiffness matrix A are

found from the approximations of d", u and v. The stiffness matrix is of dimension

2N 2 x 2N2 (N = 21) and has entries of the form

_ 11 (1.0.1.0) 33 (0.1.0.1) 3 (1.0.0.1) 13 (0.1.1.0)
(All)(kl)(mn) _ Z dij‘qulmn + dquijqklmn + dJl'quJ'qulmn + dqucjmklmn (4'33)
m

_ 13 (1.0.1.0) 23 (0.1.0.1) 33 (1.0.0.1) 12 (0.1.1.0)
(A12)(kl)(mn) — Z djmcquklmn + dim J'qulmn + dim J'qulmn + dijjmklmn (4'34)
M

67

_ 13 (1,0,1,0) 23 (0.1.0.1) 12 (1.0.0.1) 33 (0.1.1.0)
(A21)(kl)(mn) " Z dijquklmn + dqucquklmn + dquijqklmn + djmcquklmn (4'35)
m

_ 33 (1.0.1.0) 22 (0.1.0.1) 23 (1.0.0.1) 23 (0.1.1.0)
(A22)(kl)(mn) _ Z dquCJ'qulmn + dquCJ'qulmn + dquijqklmn + djmcquklmn (4'36)
M

where
$3435.73 = (/n :.<z)D°¢:.(x)D<¢;+m(4)49.) x (4.37)
(/n 24.4w" 35(3)” tow.)
and Di denotes the derivative operator.
Terms Cgfligfg in (4.33) through (4.37) are products of connection coefficients and

depend only on the choice of basis functions used to approximate d", u, and v. Con-

nection coefficients for each type of wavelet may be computed once for all problems

 

and stored as data. The entries in A are computed by a material weighted sum of
connection coefficients, eliminating the need for the numerical integration associated
with tradition finite element approaches. Due to the compact support of wavelets
and their derivatives, the number of non-zero terms in the summation is only de-
pendent on the type of wavelet used. In addition, there exist accurate and efficient
methods for computing the connection coefficients, see Latto et al [23] and Dahmen
and Micchelli [11] for details.

As in the case of finite element stiffness matrices, stiffness matrices associated
with wavelet approximations, A, are symmetric, sparse and handed. Examples of the
structure of two-dimensional stiffness matrices using Daubechies D6, quadratic spline
S4 and linear spline S3 approximations for u and v are shown in Figures 4.6, 4.7,
and 4.8 respectively.

Remark: The stiffness matrix A is positive semi-definitive because only periodic

68

 

Figure 4.6: Structure of A when using Daubechies D6 fixed scale expansions

69

 

 

 

 

 

 

 

 

 

Figure 4.7: Structure of A when using quadratic spline S4 fixed scale expansions

70

 

 

 

 

 

Figure 4.8: Structure of A when using linear spline S3 fixed scale expansions

71

boundary conditions are imposed. The matrix A has two zero eigenvalues with eigen-

vectors

p = (4.38)

These eigenvectors represent rigid body modes.

Remark: The overall banded structure is caused by the periodicity in Q and the
selection of periodic wavelet basis functions for u and v. The number of possible non-
zero entries in the stiffness matrix is directly related to the support of the wavelet
used to approximate the displacement field. From Figures 4.6, 4.7, and 4.8 we notice
that the number of non-zero entries decreases as the support of the wavelet decreases.
Therefore, to promote sparse stiffness matrices, wavelets with the shortest support

possible should be used to approximate the displacement field.

Constructing stiffness matrices associated with homogeneous material dis-

tributions

The stiffness matrix associated with a homogeneous material distribution in Q has
some special properties. The homogeneous material stiffness matrix, Ah, is parti-
tioned in a similar manner as the stiffness matrix, A, namely
h h
A11 A12

A), = (4.39)
A53 A32

72

Entries in A), are of the form

Jklmn jkImn jkImn jkImn

(Alll)(k()(mn): d110(1 010)+ (1330(0’1’0 1) + d13C(10 01) + d13C(0110) (4.40)

(A3,)(HW): (113011; H0 1 0’ + c1230”; 1 ° 1) + d33C(-l}’ H0 0 I) + d12C(-°” ,.1 10’ (4.41)

jkl mu 1]: Imn 3k Imn jlclmn

(A1211)(H)(mn) _d13C(1010) + d23C(0101) + d12C(1 00 1) + d330(01 1 0) (4.42)

JkImn jkImn jkImn jkImn

jkImn 3'”: mu 1]: Imn 1k Imn

(A3,) (Wm =d33C(-1” H” 1 0’ + (1220(0): 1 ° 1’ + c1300}: ° ° 1) + d23C(-°;’ ..1 1 0’ (4.43)

(0.19.90)

where terms C Hm" are again connection coefficients, but of the form

333:3: (/D°¢-.(x)n<¢;..(x)dn.) ((1034:.(y)m¢;+.(y)da) (4.44)

Examples of the structure of stiffness matrices A), associated with various wavelet
approximations are shown in Figures 4.9, 4.10, and 4.11.

Now we examine some special properties of the stiffness matrices associated with
a homogeneous material distribution. Each of the sub-matrices of Ah, A?“ A’l‘z, A31,
and A32, is block circulant. To understand the structure of a block circulant matrix,
the structure of a circulant matrix is examined. A circulant matrix may be completely
defined by a single row. Every other row in the matrix may be obtained by a right
circulant shift of the given row. For example, the 2"" + 1 row of a circulant matrix
may be obtained by a right circulant shift of one unit of the 2"" row. A block circulant
matrix is a matrix composed of M 2, M X M circulant matrices. The 2"” + 1 block

row is obtained by a right circulant shift of one unit block of the 2"" block row. Block

73

 

 

 

 

 

Figure 4.9: Structure of A), when using Daubechies D6 fixed scale expansions

74

 

 

 

 

 

 

 

 

 

Figure 4.10: Structure of A), when using quadratic spline S4 fixed scale expansions

75

 

 

 

 

 

Figure 4.11: Structure of A). when using linear spline S3 fixed scale expansions

76

circulant matrices can be stored very efficiently (only one row is needed) and inverted
very easily using Fourier transforms [8]. This is because block circulant matrices are
diagonal in the Fourier basis. Recalling that each sub-matrix of A), is block circulant,
A), is tridiagonal when represented in the Fourier space.

Remark: As in the case with A, A), has two zero eigenvalues associated the rigid
body modes in (4.38).

Remark: An interesting trend is noted when comparing the structure of stifl'ness
matrices associated with heterogeneous material distributions with the structure of
the stiffness matrices associated with homogeneous material distributions. When the
Daubechies D6 approximations are used for the displacement field, A). contains far
fewer entries than A. This is due to the orthogonality of the Daubechies D6 wavelets.

Due to the constant material properties in the homogeneous setting, the inner prod-

(1 ,0, 1 ,0)
CJkImn

not between non-identical Daubechies wavelets vanish, i.e., the terms and

o,1,o,1 . 0 .
Cfikmn ) are zero unless I = n and k = m respectively. However, 1n a heterogeneous

material distribution with rapidly varying material properties, the inner product be-
tween non—identical wavelets usually vanishes only if the supports of the wavelets
do not intersect. This is not the case with the spline wavelets since they are non-

orthogonal basis functions. This is apparent in the similarities seen when comparing

matrices A and A), corresponding for the spline approximations for the displacement

field.

77

Constructing the constraint matrix

The constraint matrix B describes the essential boundary conditions. Since the bound-
ary integrals are replaced with domain integrals, entries of B are determined from
examining the displacement field along the discretized boundary AI‘". Once again,
connection coefficients are used to determine entries in B

_ hour 11
B " AI‘“ ¢jkI jmndQ (445)

Note that the Haar wavelets (fig-TI" are associated with the approximation of the La-
grange multipliers /\ along the Aw. Since 8w is approximated by one pixel width

boundary, the constraints in Bu = 0 force the average displacement over each pixel

along the approximate boundary AP“ to be zero.

Constructing the load vector

Based on the term [a fvdfl in (4.30), the load vector becomes

21-121—1
F(mn) = 2 z fjkICjklmn (4.46)
k=0 (:0
where
0,-1.1"... = £1,393? ymndo (4.47)

Recall the coefficients fjk, are associated with the Haar approximation of the body

load f expressed in (4.27).

78

Summary

Using wavelet approximations, a discretized version of the material distribution, dis-
placement field, boundary conditions, and loading conditions were introduced. These
approximations were used to construct a discrete form of the Fictitious Domain Prob-
lem. We now discuss an efficient method to solve the system of equations associated

with the discrete Fictitious Domain Problem.

4.2.4 Alternative Fictitious Domain Problem Formulation

This section discusses an approach that results in an alternative form of the discrete
problem that facilitates its solution. The discretized system of equations associated

with the Fictitious Domain Problem is

ABT U F

(4.48)
B 0 A 0
For purposes of this discussion, it is convenient to view these equations as arising

from the problem

Find U E 3321‘” that
minimizes %UTAU — UTF (4-49)

subject to BU = 0

79

Now, we define the following mean value zero vectors:

U = U — pa
(4.50)
F=F-pfl

where p are the rigid body modes in (4.38). The vectors U, F, a and 6 are of the

form

- 111 .. f1 C11 [51
U: , F: , a: and 6: (4.51)

112 f 2 a2 ,32

where a1, a2, 61 and 62 are set to guarantee that 211, 112, f1 and f2 have mean value

zero respectively. Using these mean value zero vectors, we rewrite (4.49) as

Find U E 322”: and a1, a2 6 R that
minimizes %(U + pa)TA(U + pa) -— UTF -- aTfi (4-52)

subject to B(U + pa) = 0

Since functions 21,- and f,- have mean value zero, pTU = 0 and pTF = 0. Also,
since p contains the eigenvectors of A associated with zero value eigenvalues, Ap = 0.

This reduces (4.52) to the problem

Find U E §RzN2 and a1,a2 E R that
minimizes %UTAU — UTF — aTfi (4-53)

subject to BU + ga = 0
where ga E Bp. Note that problems (4.49), (4.52) and (4.53) are equivalent prob-

80

lems.

We now add the following penalty term to the objective function of (4.53)
1 “ T “ 1 T
P = -2-I~tU QU + 56(BU) (BU) (4.54)

where Q = ppT and K. and e > 0 are scalars. One may use a mean value representation

for the constraint matrix B, i.e.

B = B + BQ (455)
With this relationship, ones finds
BU = BU + Bpa = BU +ga (4.56)
This results in the problem
Find U E 322”? and a1, a2 6 R that
minimizes fiUTAU — UTF — aTfl + P (4-57)

subject to BU + ga = 0

which is equivalent to (4.53) since at a solution of (4.53), QU = 0 and BU = 0 and

the penalty term P vanishes. The stationary conditions of the Lagrangian associated

81

with problem (4.57) yield the following discrete system of equations.

 

A + nQ + EBTB eBTg BT U F
eBTg eng gT a = )6 (4-58)
B g 0 A 0

 

where A are Lagrange multipliers. Examining (4.58) more closely, one observes that

gTA = ,6. Using this relationship in (4.58), we find

(A + 50 + 613"}? — BTg(ng)"ng§) 9T - I"3Tg(ng)"gT U
13 - g(ng)“gT13 793T - %g(ng)“gT A
13‘ - 9Tg(ng)‘lfi

'ng - %g(sTg)"fl
(4.59)

where 'y is an arbitrary scalar.
Remark: System (4.59) is equivalent to system (4.31) with arbitrary scalars n and
7 provided 6 > 0 . However, as upcoming sections will show, this system provides

more flexibility in the implementation of the solution technique.

4.2.5 Preconditioned Conjugate Gradient Solution Tech-
nique

We investigate the solution of system (4.59) using the Preconditioned Conjugate

Gradient (PCG) algorithm of Bramble and Pasciak [5]. The algorithm is designed to

82

solve problems of the form

'11:

A BT U
=- (4.60)

B—C A

Q:

where A is positive definite and C is positive semi-definite. Bramble and Pasciak

show that a preconditioner of the form

P = (4.61)

leads to system of equations with condition numbers approaching unity for an appro-

priate choice of A0. A “good” A0 satisfies the condition

0 < ((A — A,))U, U) _<_ n(A_U, U) for all U 94 0, U e 322”“ (4.62)

for some constant 17. In general, A0 will be “close” to satisfying (4.62) if
1. A, > 0
2. A — A, 2 0
3. BAglfiT > 0
In addition, for the preconditioner to be useful in computations we require:
4. The spectrum of A0 closely approximates the spectrum of A

5. A0 is easily invertible

83

A preconditioner involving A). is a good candidate, as A), and A have similar structure
(see figures 4.6 through 4.8 and figures 4.9 through 4.11) and is block-circulant and
thus easy to invert. However, the spectrum of A), differs from the spectrum of A

because of the terms in A that result from the boundary conditions. 80, we choose

A0 = A), + nQ + eBTB — BTg(ng)”1gTB (4.63)

where A), is constructed with a weak isotropic homogeneous material distribution
over the entire domain 9. This choice of A0 satisfies the conditions (1-5) above and
results in an efficient preconditioner.

The introduction of the preconditioner P leads to the following system of equa-

tions:
0 ..
M =F (4.64)
A
where
A;1 0 A ET
M :
BA? —1 E —C
(4.65)
AglA Agll'BT
13434—13 mysuc
and
~ AglF
F = (4.66)

EAglF — G

84

Note that M is not symmetric, but it is positive definite and symmetric with respect

to the weighted inner product

, = ((A — Ao)w, y)§R2 + (:r, z)g;2 (4.67)

Bramble and Pasciak proved that system (4.64) has a condition number that is in-
dependent of resolution. Therefore, system (4.64) could be solved using a PCG al-
gorithm based on the weighted inner product defined in (4.67) at convergence rates
independent of scale (but depend on the constant 1)).

Determining constants KS, 6 and 7

The constants It, 6 and 7 are selected as follows:
1. Select 6 and 7 so the matrix C = —7ggT+%g(ng)'1gT is positive semi-definite.

2. Given 6 and 7, select It so the matrices A = A+nQ+£BTB —cBTg(ng)'1gTB

and A0 = A). + «Q + cBTB — eBTg(ng)"lgTB are positive definite.

With proper selection of It, 6 and 7 we may now investigate how to use the primary

preconditioning matrix A0.

85

Inverting A0

The matrix A0 must be inverted to apply the preconditioner P. This is done with a

Sherman-Morrison-Woodbury formula [1]
(A + UVT)-l = .4-1 — A-1U(I + 1/TA-1U)-1VTA-l (4.68)

We express A, as

A0 2 Ah + ch + eBTB — BTg(ng)'1gTB
= (Ah + 6Q) + lA31"(<51 - g(ng)"gT) 1’»
(4.69)
= (Ah + 6Q) + (BTLXLTB)

= A+UVT

Remark: Using (4.68) and (4.69) require the construction, storage, and inverse of

a square matrix of dimensions equal to the total number of constraints.

4.3 Post Processing Wavelet-Based Solutions

This section discusses the post processing of the fixed-scaled wavelet expansion of the
displacement field to determine the strain and stress fields. The wavelet expansion is

used to determine the average displacement, strain, and stress over each pixel in the

86

domain w. The average displacement over each pixel X (a, b) is given by

1
-. _ __ . — 1 Q
(11,)(ab) ( ,b) / (a,b) u,(:1:,y)d$2 2_ 711/11 x( a, b) )u- (21:, y)d

2J-12J1

= 2’ )3 2114...] 31231449114. 1140 (4.70)
Ic=0 I: O
2J—12J-l C(00)

= 2j Z 2 “1°qu jabkl
k=0 I: 0

where A(a, b) = 2‘” is the area of pixel X (a, b) and x(a, b) is an indicator function
corresponding to the area occupied by the pixel X (a, b). The terms 0‘35} are again

connection coefficients of the form

05-31532 = ([4 arrow“: (x)dn)><

 

(4.71)
([945 ¢?é‘“"(y)Dfi¢ ;1(y)dfly)
The strain field is computed based on the small strain relationships
1 611.5 6113'
.. : _ _ 4.72
5” 2 (an, + 62,-) ( )
The average strain field for a given pixel X (a, b) is given by
2.21—121—1 (10)
(€11)(ab) = 2 J 2 Z uljkICjabkm
Ic_=0 I=0
.21- —121'— 1 C(0 1)
(€22)(ab) = 22] Z 2 1121110 jabkm (4'73)
1:0 1: 0
- 2j-12j—1 (01) (1,0)
(€12)(ab) : 221—1 2: Z (uljklcjabkm +u2jkICjabkm)
k=0 I=0

87

where connection coefficients 0:32,) are defined in (4.71).

Remark: Expressions for the average stress field are based on the average strain

field and follow the stress/strain relationship: 6 = DE.

4.4 Summary

This chapter discussed a fictitious domain approach using wavelet-based Galerkin
techniques. A discrete Fictitious Domain Problem statement was derived. A PCG
solver is used to solve this problem using a preconditioner based on the homogeneous
stifl’ness matrix Ah. The next two chapters provide examples of topology optimization
problems in two and three dimensions that are solved using versions of PCG methods

discussed here.

88

Chapter 5

Wavelet-Based Topology

Optimization in Two Dimensions

5.1 Introduction

The previous chapter introduced wavelet-based analysis using fictitious domains. This
chapter presents examples of two-dimensional problems solved using this method.
We first address an elastostatic analysis problem to determine how effective various
wavelets are in solving the elasticity problem. The wavelet-based solution method
is integrated with standard topology Optimization, leading to wavelet-based topol-
ogy Optimization. Various topology optimization problems are solved to confirm the

effectiveness of wavelet-based methods.

89

V

 

Young’s Modulus: El

Poisson’s Ratio: \11

 

 

 

 

mam“ .

 

X

 

A
v

Figure 5.1: Example elastostatic problem definition

5.2 Elastostatic Analysis Using Wavelet-Based
Analysis

This section examines solutions of an elastostatic problem using wavelet-based analy-
sis with fictitious domains. We wish to examine how the choice of (1) the wavelet used
to approximate the displacement field and (2) the strength of the fictitious material
affect the accuracy and performance of the PCG solver. An example problem, shown
in Figure 5.1, is solved using wavelet-based and finite element analyses, using similar
domain discretization, boundary and loading conditions.

First, a wavelet-based discretization is developed for the domain w in Figure 5.1.
The domain in this example can be represented exactly using square pixels. In the

example, the domain w is discretized in raster form, using 30 rows by 48 columns of

90

 

Figure 5.2: Example problem placed in a periodic fictitious domain

square pixels. The fictitious domain discretization is created by placing the raster
representation of w into a fictitious domain 0, which is discretized using 64 by 64
square pixels. The pixels in w are associated with an isotropic material with Young’s
Modulus E1 and Poisson’s ratio 111. The pixels in the fictitious domain, 0 \ w, are
associated with an isotropic material with Young’s Modulus E2, where E2 << E1, and
Poisson’s ratio V2 = V1. The discretization of the domains w and Q are displayed
in Figure 5.2. In the wavelet-based approach, constraints and loading conditions are
enforced on an appropriate subset of the boundary pixels Aw, where Aw is a discrete

approximation of 6w. In this discretization, Aw is made up of the outer rows and

91

Loaded
Pixels

Constrained
Pixels

 

Am:-

Figure 5.3: Wavelet-based discretization

columns of pixels in w. The fixed displacement constraint on the left side of the
structure is enforced by requiring that the average displacement of each pixel along
the left edge of Aw is zero. The point load is approximated by placing a body force
over a small number of pixels neighboring the point. The body force is distributed
over two pixels located at the center (vertically) of Aw along the right edge. The
extent of Aw and the implementation of the boundary and loading conditions for the
wavelet-based approach is shown in Figure 5.3.

The finite element discretization begins with a raster description of the domain

92

 

Edge of
Constrained
Nodes

Point

 

Figure 5.4: Finite element discretization

w. The finite element mesh is created by treating each pixel as an isoparametric, four
noded square finite element. In this case, the boundary 6w is represented exactly. The
material properties of each element correspond to an isotropic material with Young’s
Modulus E1 and Poisson’s ratio 111. The fixed displacement boundary conditions along
the left edge of w are applied by constraining the degrees of freedom of appropriate
nodes. For this problem, all degrees of freedom of the 31 nodes located along the left
edge of the domain are constrained. A point load is placed on the node centered on
the right edge of w. This is a typical finite element implementation for this problem

and is shown in Figure 5.4.

5.2.1 Accuracy Of The Wavelet-Based Fictitious Domain So-
lution

Here, we examine the accuracy of the wavelet-based fictitious domain approach when

compared to the finite element discretization. The objectives are to (1) determine

93

Top Edge Pixels

Right
Edge
Pixels

 

30

Figure 5.5: The areas of comparison for the wavelet-based and finite element dis-
cretizations

the effectiveness of various wavelets used to approximate the displacement field and
(2) find an appropriate magnitude for the fictitious material’s Young’s Modulus, E2.
Evaluations are made by comparing the wavelet-based and finite element stress fields
along the top and the right side boundaries of w. In the wavelet-based approach, the
average stress field over each pixel is computed based on the average strain field as
discussed in Chapter 4, Section 4.3. The finite element stress field is computed at the
center of each element, which is standard in finite element analysis.

The stress fields are compared along the right edge of pixels/elements and along
the top edge of pixel/ elements as shown in Figure 5.5. The traction along the bound-

ary may be expressed as

tlfi) 011 012 n1
tan) 012 022 712

94

where t is the traction and f1 = [711, n2]T is the normal vector to the boundary. Along

the right edge, we have

1 011
f1 = , t = (5.2)
0 012
and along the top edge, we have
0 012
f1 = , t = (5.3)
1 022

On a traction free boundary, t should vanish. It is important for the wavelet scheme to
capture this feature accurately. A non-zero stress vector on a traction free boundary
indicates that the fictitious material is applying an appreciable load on w.

We now examine the effect of the strength of the fictitious material when using
various wavelets to approximate the displacement field. For each wavelet, Daubechies
D6, quadratic S4 spline and linear S3 spline, fictitious materials of varying strength
will be used. The material in the domain w is isotropic with E1 = 1.0 and
121 = 0.3. The fictitious material is also isotropic with varying Young’s Moduli
E2 6 {10'1, 10‘2,10‘3, 104} and V2 = V1.

We first compare the wavelet-based results using Daubechies D6 wavelets with
the results found with the finite element analysis. Figure 5.6 displays the norm of
the traction along the top and right edge of w. Figure 5.6(a) displays the norm of
the traction along the t0p row of pixels/elements within w. Notice the differences in
stresses at the first pixel/element along the top edge when comparing wavelet-based

results for various strengths of E2. As the fictitious material becomes weaker, the

95

 

E PEN} 9-
6 = 10- + -
'1) E: = 10-3 U

5 - E2 = 10-2 X ._
4 I: _

”tn 5

 

 

 

 

 

 

24
Element / Pixel Along Top Edge

(a) Norm of the traction along top edge of w

 

 

 

18 I c-. l
16 _ FEM '9'- _
= 13*: +
_ 2 = 1 _ D
12 1- E2 = 10-1 A .1
10 - '—
lltll
— 0 0 -1
6 *- _.
P 9
4F OQAAQO
24% .. {,3 QA AQ 3 - z:
t f. 4 A 4-52.4345; -"- A A ‘4' *7 A A ‘ a 1'

0 1.‘ . . ., J
1 10 20 30
Element / Pixel Along Right Edge

 

(b) Norm of the traction along right edge of w

Figure 5.6: Norm of the traction along the top and right edge for D6 wavelet-based
and finite element solution techniques

96

wavelet—based stress approaches the finite element stress value, but still remains much
lower in magnitude. This type of behavior is present in all wavelet-based solutions
presented here, but this occurs only in a subset of Aw, which for large-scale problems,
is very small compared to the problem domain w.

Figure 5.6(b) displays the norm of the traction along the right edge of pix-
els/elements within w. The figure shows that neither method models the stress free
condition along the right edge of w. The effects of the non-zero stress near the load
are present over a majority of the boundary. This type of response is believed to
be primarily due to the coarse discretization used in this example. However, as the
strength of the fictitious material is reduced, the wavelet-based results better approxi-
mate the finite element results. Note that the Dacbechies D6 wavelets do not produce
symmetric results along the right edge of w. This loss of symmetry is related to the
non-symmetric nature of the Daubechies D6 wavelet. When using a finer discretiza-
tion for the domain w, the traction free boundary along the right edge away from
the load is more accurately modeled, but is still not symmetric. This is shown for a
192 x 120 pixel discretization of w in Figure 5.7.

We now examine the effectiveness of the spline S4 and S3 wavelets in modeling
the stress free boundaries. In general, the spline S4 and 83 results model the stress
condition similar to the Daubechies D6 wavelet, but slight differences are present.
Figure 5.8 illustrates a comparison of the norm of the traction along the right edge of
w for the Daubechies D6, spline S4, and spline S3 wavelets at a fixed weak material,
E2 = 10". From Figure 5.8 we see that unlike the Daubechies D6, both the spline

S3 and S4 wavelets produce symmetric results.

97

 

 

 

 

20 1 1
13 ~ -
16 - -
14 - -
12 - ~
Iltll 10 - r
8 _ -
6 .— —1
4 u- q
2 — —1
0 _/
1 40 80 120

Element / Pixel Along Right Edge

Figure 5.7: Norm of the traction along the right edge for a 192 x 120 pixel discretiza-
tion and D6 wavelet-based solution

We now investigate the performance of these different wavelets in solving the
elastostatic problem where traditional finite element analysis is replaced by wavelet-

based analysis.

5.2.2 Performance Of The Wavelet-Based Fictitious Domain

Solution Method

In this section, we discuss the performance of the solution method described in Chap-
ter 4 when solving the elastostatic problem depicted in Figure 5.1 for the three
wavelets discussed in the previous section. Recall that the solution scheme is centered
on the use of a PCG solver. The performance of the solution method is determined
by the number of iterations required for the PCG solver to meet a specified error

tolerance. The error is measured by a normalized norm of the residual, i.e., er-

98

 

 

 

 

 

 

l 24 48
Element / Pixel Along Top Edge

(a) Norm of the traction along top edge of w

 

lltll

w
l

 

 

 

0 {7 ‘ ‘3’ 3:. 1: 1: ‘ - I ' _ _
l 10 20 30
Element / Pixel Along Right Edge

 

(b) Norm of the traction along right edge of w

Figure 5.8: Norm of the [traction along the top and right edge for various wavelets
and E2 = 1.0—4

99

 

 

 

1 I I I I I I E
(4 D6 49— :
~ ,, S4 -l— ‘
" “ ~ 4* S3 8— 1
0.1 .,: ‘3
“at i
h°¢°4 . ‘
[11,11] “ I " ' ° 9 g l
[F,F] 0.01 \s “ E
\‘
u *e’e
\\ . ’4 9 “ d
0.001 n V s‘ 4 -.
i
.4
0.0001 7 l I l l l l
0 5 10 15 20 25 30 35
PCG Iteration

Figure 5.9: The performance of the PCG solver using various wavelet approximations

ror = [R, R] / [F, F], where [~, -] is defined in (4.67). Figure 5.9 displays magnitude of
error during the PCG iteration history of the problem solved with E2 = 10‘4 using
the Daubechies D6 wavelet, spline S4 wavelet, and spline S3 wavelet approximations.
For this example, the initial guess was U = 0 over 9.

The figure shows that better performance is achieved using the S3 wavelets. This
behavior is typical and was observed in other problems with different geometries and
boundary conditions. It is speculated that 83 wavelets lead to better preconditioners.

It is interesting to compare the performance of the PCG solver in the wavelet-based
approach with traditional iterative (using a diagonal preconditioner) finite element
methods. The domain w depicted in Figure 5.1 is discretized at three resolutions
as described in Table 5.1, where w DOF corresponds to the number of degrees of
freedom associated with the displacement field in w and Total DOF corresponds to

the number of degrees of freedom associated with the displacement field over the

100

fictitious domain. The PCG performance for the wavelet—based method using linear

Table 5.1: Discretization levels for analysis of PCG performance

 

 

 

 

 

 

 

Level w Discretization w DOF Q Discretization Total DOF
6 30 x 48 2,880 64 x 64 8,192
7 60 x 96 11,520 128 x 128 32,768
8 120 x 192 46,080 256 x 256 131,072

 

 

 

 

 

 

 

83 spline wavelets at the three discretizations is shown in Figure 5.10(a). While the
PCG performance for the traditional iterative finite element methods for the same
discretization levels is shown in Figure 5.10(b). (As is the previous examples, the
fictitious domain approach is not used in the finite element solution technique.) From

Figure 5.10 we note a few important trends.

1. The PCG performance of the wavelet-based approach does not decline as the
resolution level increases. The wavelet-based approach exhibits level insensitive

convergence properties.
2. Convergence using wavelet based schemes is typically monotonic.

3. As is typical with iterative finite element approaches, significant computational
effort is spent in a large portion of the iterations without significant reductions

in the error.

4. The PCG performance of the traditional iterative finite element decays as the

resolution level increases.

This comparison of the PCG performance for the two solution methods shows one of
the potential benefits of the wavelet-based solution method over tradition iterative
finite element methods: monotonic convergence that is independent of resolution.

101

 

1% I F l l I r

 

Level 6 G-
.‘:\,_- .. .. five] 7 g:
aru§g=gg V8 8
0.1 \x
l3\s:\
[11,11] \-:':
[F.F] 0'01 ‘\‘,
e a
0.001 9 F
El
(A0001 7 ' ' 1 l I 1 ,
O 2 4 6 8 10 12 14

PCG Iteration

(a) Performance of the wavelet-based PCG solver using linear S3 Spline wavelets for
various levels of resolution

 

11).“ a

 

 

|R,R]
0.001 1:
Level 6 Level 7 Level 8 ;
00001 I I I I I I
0 100 200 300 400 500 600 700

PCG Iteration

(b) Performance of a finite element PCG solver for various levels of resolution

Figure 5.10: Performance of the PCG solver for wavelet-based methods using 83
spline wavelets and traditional iterative finite element methods

102

5.3 Solving Topology Optimization Problems Us-

ing Wavelet-Based Analysis

We now incorporate the wavelet-based fictitious domain solution strategy and topol-

ogy optimization using the following algorithm:

1. Embed the package space w into a sufficiently large square fictitious domain (2,

filling Q \ w with a suitable weak fictitious material.

2. Generate a raster description of the fictitious domain of size (N x N) where

N=2J.

3. Define the design variables based on the raster description of the package space

w. The design variables depend upon the material model used in the topology

optimization algorithm. For example, if homogenization is used with a micro-

hole material model, one uses the expansions for the cell parameters a, b and 0

of the form

9(2, y) =

b(:vy)=

0(1),?» =

21-121—1

Z Z ajk1¢???'($,y)

k=0 I=0

(5.4)

2J—12J— l

A: Z bJ-n <15???“ m y)

k==0l0

(5.5)

21-121—1

Z Z 93°1c1¢2131ar($, y)

k=0 I=0

(5.6)

If a simplified (“density approach”) material model is used, then an expansion

for p is used, i.e.

p(x, y):

2J-12J— 1

Z 2 ijl ?£f’(x,y)

k:0l=0

(5.7)

103

In both approaches, the coefficients of the design variables for all $1,??'(:1:,y)
whose support is outside of the package space w are fixed and are set to values
that correspond to the weak fictitious material properties. Define the material

tensor d’"" corresponding to the design variables.

Following this algorithm, one may solve topology optimization problems using

wavelet-based analysis. However, a few questions must still be addressed:

1. Which wavelet, Daubechies D6, spline S4 or spline S3, should be used to ap-

proximate the displacement field?
2. What is a suitable strength for the fictitious material?

3. What characterizes the size of a sufficiently large square fictitious domain?

5.4 Adjusting Wavelet-Based Analysis For Topol-
ogy Optimization

This section discusses how the selection of the (1) wavelet for the displacement field
approximation, (2) strength of the fictitious material, and (3) extent of the ficti-
tious material affect the results obtained using wavelet-based analysis with topology

optimization.

104

5.4.1 Effect Of Wavelet Used For Approximating The Dis-
placement Field: Checkerboard Patterns

From the elastostatic analysis problem, we saw the various advantages and disadvan-
tages in using the Daubechies D6 wavelet, spline S4 wavelet, or spline S3 wavelet for
analysis. This investigation determines how these wavelets perform in a topology op-
timization setting. The problem domain of interest is the same as used in elastostatic
analysis and is displayed in Figure 5.1. The domain w is selected for the package space
with a volume constraint V = 0.35 (meaning that 35% of the domain may be filled
with material in the optimal design.) For the problems addressed here, the maximum
strength of the design material allowed corresponds to an isotropic material with
Young’s Modulus E = 1.0 and Poisson’s ratio V = 0.3. Using a raster description of
the problem, a wavelet-based discretization is generated based on a fictitious domain
of size 256 x 256, with a weak material E2 = 10" and V2 = 0.3. The design domain
w is discretized at a resolution of 120 x 192 pixels. The point load is approximated
with a body force distributed over two pixels centered along the right edge of the
discretized boundary Aw. Figure 5.11 displays the optimal topologies found using
the wavelets of interest.

The three results displayed in Figure 5.11 are similar to one another and resem-
ble the optimal topology reported in the literature [22]. However, after examining
each result closely, one notices the differences between the three topologies and finds
how they deviate from the results reported in the literature. The results based on

Daubechies wavelets are slightly non-symmetric and contain regions of checkerboard-

105

 

(a) Daubechies D6 (b) Spline S4

 

(c) Spline S3

Figure 5.11: Optimal topologies using various wavelet approximations for the dis-
placement field

106

 

E=0.0005

y ‘ I 20‘) Y

 

 

X X
(a) Material (Haar) and displacement (S4) (b) Material (Haar) one level lower
at the same level than displacement (S4)

Figure 5.12: Material mixture in layer form for two approximations

like material distribution. These checkerboard patterns resemble those found when
using isoparametric, four noded finite elements to solve topology optimization prob-
lems [14]. The results based on the spline S4 wavelet approximation are symmetric
but exhibit an unusual sub-optimal layered material distribution. The results based
on the spline S3 wavelet approximation are symmetric but contain regions that exhibit
sub-optimal checkerboards.

The checkerboards and layered material patterns result from numerical phe-
nomenon that arise in the interaction between the displacement approximation and
material approximation. This phenomenon may be explained using wavelets to de-
termine the effective properties of a two material mixture that is mixed at scales
determined by the numerical approximation of the displacement and material fields.

When the mixture occurs at the same scale (material and displacement at the same

107

level), the effective properties of a patch as shown in Figure 5.12(a) are

 

 

0.50 0.15 0
D= 0.15 0.50 0 (5-8)
0 0 0.175

The effective material properties displayed in (5.8) show that the layered material
is artificially stiff in the layered direction. This artificial stiffness results from the
fact that, in this case, the material distribution varies faster than the support of the
displacement field. This artificial stiffness may be eliminated if one decreases the level
of approximation for the material distribution as displayed in Figure 5.12(b). With
a level 3' — 1 approximation of the material distribution and a level j approximation

of the displacement field, the effective material properties of the layer material are

0.0013 0.0004 0

D= 0.0004 0.46 0 (5-9)

 

 

0 0 0.00047

which better reflect the expected properties of the layered material. This finding
leads to a straight forward method to eliminate the checkerboard and layered material
patterns found in the optimal topologies reported earlier. Namely, approximate the
material distribution (and the domain) one level lower than the approximation level

of the displacement field. This leads to the following design variable approximations

108

when using homogenization:

20-11-120-1)—1

000,31): 2: Z a(j—1)k1¢[’ffl)n($,y) (5-10)

k=0 l=0

2(j-l)_12(j-1)_1

144.1): 2 Z 116-61:41:35,..(44) (5.11)

k=0 I=0

When using the fictitious material model, this leads to the following design variable

approximation:
20-11-120-1L1
P($,y)= Z Z 10(j-1)k1¢(ffi)1¢1($,y) (5-12)
k=0 I=0

From experience, we have found the most success approximating the microcell orien-
tation 0 at the same level as displacement. Using an approximation of the material
distribution one level lower than that of the displacement field leads to the optimal
topologies shown in Figure 5.13. Although adjusting the relative discretization levels
corrected the checkerboard patterns found in the Daubechies D6 results, the results

displaced in Figure 5.13 are still slightly non-symmetric.

Conclusions

To avoid the presence of checkerboards and layered material in the Optimal topology,

the material field must be approximated one level lower than the displacement field.

5.4.2 Effect Of The Strength Of The Weak Material

This section examines the effect Of the fictitious material strength on the Optimal

topology. Again we use the design domain shown in Figure 5.1 with a volume con-

109

  

V
A

(a) Daubechies D6 (b) Spline S4

)0

(c) Spline S3

    

 

Figure 5.13: Optimal topologies found when approximating the material field one
level lower than the displacement field

110

  
 

(b) E2 = 10—3E1

(c) E2 = 10-4131

Figure 5.14: The effect Of the strength of the fictitious material on the optimal topol-
ogy

straint of V = 0.35. Here, we examine solutions for fictitious materials of various
strengths: E2 = 10‘2E1, E2 = 10‘3E1 and E2 = 10‘4E1. The Optimal topologies
using the spline S4 wavelet level j = 8 displacement approximation and Haar wavelet
level j = 7 material approximation are shown in Figure 5.14. From Figure 5.14 we
see that the solution to topology optimization problems are definitely sensitive to the
strength of the fictitious material. As the strength Of the fictitious material increases,
the periodicity assumption allows the loads to be transmitted through the fictitious

material, to the structure in nei hborin “ eriodic domains.” Non-zero tractions will
8 g P

111

appear on what should be traction free boundaries, causing erroneous results.

Conclusions

The strength Of the fictitious domain is selected no greater than E2 = 10‘4E1 to

guarantee negligible tractions on traction free boundaries.

5.4.3 Effect Of The Extent Of The Fictitious Domain

One Of the drawbacks of a fictitious domain approach is that it increases the number
of degrees of freedom needed to solve the problem by adding the fictitious domain
surrounding the domain w. Here, we investigate the effect Of varying the extent Of
the fictitious material surrounding the design domain. We examine various design
domain sizes for the bracket problem depicted in Figure 5.15. The problem is solved
at displacement discretization level 3' = 7 (material discretization at level j : 6) using
a spline S4 wavelet expansion. The problem is solved, for the design domains sized
96 x 96, 112 x 112, 120 x 120 and 124 x 124, each centered in a 128 x 128 domain
(Note the number of fictitious domain pixels surrounding the design is no less than
half of the support of the wavelet used to approximate the displacement field.) The

results for these problems with volume constraint V = 0.35 are shown in Figure 5.16.

Conclusions

The results show that the extent of the fictitious material surrounding the design

domain has little effect, if any, on the quality Of the solution Obtained.

112

 

31>:

Design

Domain

a} v

P

 

 

 

 

Figure 5.15: Bracket problem definition

5.5 Performance Of Wavelet-Based Topology Op-
timization

The performance Of the wavelet-based algorithm will be evaluated based on the num-
ber Of PCG iterations required to solve the elasticity problem once during each it-
eration of the topology optimization sequence. In the remaining examples in this
chapter, the discretization for the material is assumed to be one level lower than the
displacement discretization. Also, in all examples, the number Of fictitious domain
pixels surrounding the design domain is no less than the half of the support of the

function used to approximate the displacement field.

113

 

 

 

 

 

 

 

(a) 96 x 96 (b) 112 x 112

 

 

 

 

 

 

 

(c) 120 x 120 (d) 124 x 124

Figure 5.16: Optimal topologies for varying amounts Of fictitious material

114

 

(a) Level 6 (b) Level 7

(c) Level 8

Figure 5.17: Solutions of the 8x5 design domain problem for varying levels Of resolu—
tion

5.5.1 Effects Of Varying Resolution

The wavelet-based optimal topology results for the design domain depicted in Fig-
ure 5.1 for varying level Of resolution are examined. The design domain is embedded
into fictitious domains Of varying size corresponding to level j = 6 (64 x 64), j = 7
(128 x 128) and j = 8 (256 x 256). The degrees of freedom associated with these
discretizations are shown in Table 5.1 on page 101. All problems are solved with a
volume constraint of V = 0.35. Figure 5.17 displays the Optimal topologies for these

three examples. The performance of the wavelet-based solver is also examined for

115

 

1'00 D I I I I I I I

 

 

 

 

90 .13 Level 6 O -
+45 Level 7 +

80 - Level 8 El -

70 _ -
O

60 - —1

PCG 50 _ O _
Iterations @

4O [-

30$

20 -

10 -

0 I I I I I I
0 5 10 15 20 25 30 35 40

Topology Optimization Iteration

Figure 5.18: Number of PCG iterations required for solution of the 8x5 design domain
at varying levels Of resolution

these problems. Figure 5.18 displays how many PCG iterations are necessary during
each iteration Of the topology Optimization. Here we see that the performance Of the

PCG solver is relatively independent Of resolution.

 

5.5.2 Effects Of Varying The Extent Of The Fictitious Ma-

terial

In the example presented in section 5.4.3, we saw that varying the extent Of the
fictitious material had little effect, if any, on the optimal tOpologies found (See Fig-
ure 5.16). Now we would like to see how the extent of the fictitious domains affect
the performance of the PCG solver. Figure 5.19 displays the PCG performance for
bracket design domain sizes: 96 x 96, 112 x 112, 120 x 120 and 124 x 124. Again,

each is centered in a 128 x 128 domain. From the figure we see that as the design size

116

 

 

 

 

100 T 1 l l I
90 _+ 96 x 96 0
112 x 112 +
80 ,_ 120 X 120 CI _
O 124 x 124 x
70 ->< -
POO 6("El ‘
Iterat1ons 50 _ O _
wx *EX
E? E] x x -
fl +QQ é x¥
” 0 090390 3 x D ‘
20 F- og 68 X9 ""
m - . ” 1&99424949994911
15 20 25 30

0 5 10
Topology Optimization Iteration

Figure 5.19: Performance of the PCG solver for various design domain sizes

increases, the performance Of the PCG solver declines only for the first few iterations
of the topology Optimization sequence.
Conclusions

The examples provided show that the PCG solver converges at rates that are rela-
tively insensitive to the level of resolution and to the extent of the fictitious material

surrounding the design domain.

5.6 Wavelet-Based Topology Optimization Exam-

ples

This section presents more examples of problems solved using wavelet-based tOpology

Optimization.

117

 

Design

Domain 1
V

Figure 5.20: Problem statement for truss example 1

 

 

 

 

P

5.6.1 Single Load Case Problems

This section discusses three single load case problems solved using wavelet-based

topology Optimization.

Truss example 1

Figure 5.20 displays the problem definition for this example. The problem is solved
with a volume constraint of V = 0.35. Spline S4 wavelets are used to discretize a
design domain embedded into a fictitious domain at a level j = 8 resolution. The
discretization results in 131,072 displacement degrees of freedom with approximately
50,000 describing the displacement field inside the design domain. Figure 5.21 displays
the Optimal topology, performance of the PCG solver and the Optimization objective

function (compliance) history. From the figure, we see that the PCG solver requires

118

 

 

 

 

 

 

 

 

 

 

 

70 I I I I I I
PCG
Iterations
0 5 10 15 20 25 30
Topology Optimization Iteration
(b) Performance of PCG solver
350 I I I I I -]
300
250 J
- 200 4
Compliance 150 _ _
100 - -
50 1. _
0 l I 1 1 1
0 5 10 15 20 25 30

Topology Optimization Iteration

(c) Compliance history

Figure 5.21: Results of topology optimization for truss example 1

119

the highest number of iterations solving the elasticity problem associated with first
and second iterations (The initial PCG solver iteration is the 0‘” iteration.) The
performance Of the PCG solver improves after those iterations. This type of behavior

is typical and is present in all the examples solved with the wavelet-based solution

strategy.

120

 

Design

Domain

:4 a"

A ‘
‘ r

2

 

 

 

 

Figure 5.22: Problem statement for truss example 2

Truss example 2

This example is based on the same design domain as the previous example but has
different loading conditions as shown in Figure 5.22. The results and performance for

this problem as shown in Figure 5.23.

121

PCG
Iterations

350
300
250
Compliance 200
100

50

0

 

 

I

 

 

10 15 20
Topology Optimization Iteration

(b) Performance Of PCG solver

 

 

150 -

 

 

I l l l 1

IIIIII

 

 

10 15 20 25
Topology Optimization Iteration

(c) Compliance history

Figure 5.23: Results of topology optimization for truss example 2

122

 

Design

Domain

 

 

 

 

 

V

A? A}

Figure 5.24: Problem statement for the L domain example example

 

L domain example

This example illustrates the performance Of wavelet-based topology optimization in
a problem with a slightly more complex geometry. The problem definition is shown
in Figure 5.24. The problem is solved with a volume constraint Of V = 0.35 at a
level Of j = 8. This results in approximately 32,000 degrees Of freedom associated
with the displacement field in the design domain. The results and performance Of the
wavelet-based topology optimization are displayed in Figure 5.25. In this example, the
more complex design domain leads to slightly less efficient PCG solver performance.

However, the number of iterations required for the PCG solver follows a similar trend

123

 

(a) Optimal Topology

 

140 1 l I I T I

120
100

PCG 80
Iterations 60

r- -
40 .l
20 T WWW—r
0 1 7 7 1 1 1 7 1 1
0 5 10 15 20 25 30
Topology Optimization Iteration

 

 

 

 

 

 

 

 

 

 

 

(b) Performance of PCG solver

 

Compliance

IIII
LIIIII

 

 

 

0 I L I I J

0 5 10 15 20 25 30
Topology Optimization Iteration

 

(c) Compliance history

Figure 5.25: Results of topology Optimization for the L domain example

124

as seen in the previous examples.

5.6.2 Multiple Load Case Problems

This section discusses three multiple load case problems solved using wavelet-based

topology Optimization.

Five bar truss example

The five bar truss problem statement is shown in Figure 5.26. The two loads are
applied independently (multiple load case) and the Objective function is the average
Of the compliance associated with the individual loads. The problem is solved with
a volume constraint of V = 0.35 at a level Of j = 8. The results and performance
of the wavelet-based topology Optimization are displayed in Figure 5.27. Only one
PCG iteration history is provided since the PCG solver solves multiple load cases in

parallel.

125

 

Design
9 Domain Q

 

 

 

 

P
2

Figure 5.26: Problem statement for the five bar example

126

 

(a) Optimal Topology

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

60 I I f I I I
50 - "‘ d
40 - _l -
Iterations
20 1- -
10 f -
0 1 1 7 7 1 1 7 1 1
0 5 10 15 20 25 30
Topology Optimization Iteration
(b) Performance Of PCG solver
200 l I F r l
150 -
Mean
Compliance 100 - -
50 - .4
O L I I I I
0 5 10 15 20 25 30

Topology Optimization Iteration

(c) Compliance history

Figure 5.27: Results Of topology Optimization for the five bar example

127

A
V

 

Design

Domain

 

 

 

P P P

Figure 5.28: Problem statement for the bridge example

Bridge example

This example has three separate load cases and the problem statement is as shown
in Figure 5.28. Again, the problem is solved with a volume constraint Of V = 0.35
at a level of j = 8. The results and performance Of the wavelet-based topology
Optimization are displayed in Figure 5.29. Unlike previous examples, here there is a
relatively high number of pixels with intermediate design material, i.e., with effective
density 0 < p < 1). Again, we note the efficient performance of the PCG solver. As in

the previous example, the PCG solver solves all three loading conditions in parallel.

128

 

PCG

 

 

 

 

 

 

 

 

 

Iterations
I
0 5 10 15 2O 25 30
Topology Optimization Iteration
(b) Performance of PCG solver
300 I I I I I
250 -
Mean :23 :
Compliance
100 - -
50 r- d
0 I I I I I
0 5 10 15 20 25 30

Topology Optimization Iteration

(c) Compliance history

Figure 5.29: Results of topology Optimization for the bridge example

129

 

 

 

 

 

[l
é Design
/ Domain 5
/
é
/ 11

Figure 5.30: Problem statement for the dual load 8x5 design domain example

Dual load 8x5 design domain example

The last multiple load case example is a two load case example based on the 8x5
design domain studied previously. The problem statement for this example is shown
in Figure 5.30. The problem is solved with a volume constraint of V = 0.35 at a level
of j = 8. The results and performance of the wavelet-based topOlOgy Optimization

are displayed in Figure 5.31.

130

 

 

 

 

 

 

 

 

 

 

 

120 I I I I I I
PCG
Iterations
I I I
0 10 15 20 25 30
Topology Optimization Iteration
(b) Performance of PCG solver
600 1 r T 1 1
500 -
400 _
Compliance 300 - -
200 - -
100 - -
0 I I l I I
0 5 10 15 20 25 30

Topology Optimization Iteration

(c) Compliance history

Figure 5.31: Results of topology Optimization for the 8x5 domain dual load example

131

5.7 Summary

This chapter presented examples that illustrate the effectiveness Of topology optimiza-
tion using wavelet-based analysis. From various investigations, appropriate values for
the strength Of the fictitious material were determined. Numerous examples show the
resolution insensitive convergence properties Of the PCG solver used in the fictitious

domain approach.

132

Chapter 6

Wavelet-Based Topology

Optimization in Three Dimensions

6.1 Introduction

This chapter discusses wavelet-based topology optimization in three dimensions. An
overview of the discretization process in three dimensions is followed by the discretized
equations associated with the three-dimensional problem. Examples solved using this

method are discussed.

6.2 Fictitious Domain Problem Statement

Three-dimensional topology Optimization problems are stated in a similar fashion as
problems in two dimensions. It is assumed that a design domain or package space

w E R3 is specified with associated boundary and loading conditions. Again, we use

133

a fictitious domain approach. The domain of interest w is embedded into a cube and
the added region, Q \ w is filled with a weak material. As before, an approximate
solution for the elasticity problem in the domain w is extracted from an Q-periodic
solution. As the examples in the two dimensions showed, if the fictitious material is
sufficiently weak and the boundary and loading conditions are represented accurately,
the solution of the periodic problem in Q is a good approximation of the solution of
the original problem in w. We now discuss the discretization and the solution of the

discretized fictitious domain problem in three dimensions.

6.3 Problem Discretization In Three Dimensions

We begin by discretizing the fictitious domain (2. Here, the set (2 is discretized using
small cubes or voxels, the three—dimensional equivalent of pixels in two dimensions.
Voxels are indexed in raster form where X (k, l, m) corresponds to the voxel occupying
the unit cube [k, k + 1) x [l, l + 1) x [m, m+ 1), as shown in Figure 6.1. This results in
an N x N x N voxel discretization of 0. Again we select N = 2J for a fixed positive
integer j. The individual voxels are equivalent to three—dimensional Haar scaling
functions at level j and are used to approximate w and the discretized boundary

Of w, Aw.

134

 

 

 

 

 

 

 

 

 

 

 

/ / / / /
/
(0.0.4) (4,0 4) ///
//
/ /
/
z y (01°10) (4,0,0) /

 

 

 

 

 

 

X

Figure 6.1: Example voxel discretization

6.3.1 Approximation Of The Material Distribution And Dis-

placement Field

We first develop an approximation for the material distribution in the domain 9. In

three dimensions the material tensor D has the form

d1 1 C112 6113 d14 d15 d16
(112 (122 d23 d24 d25 d26
d13 d23 d33 d34 d35 d36
D(=v, y, z) = (51)
d” d24 d34 d44 (145 6146

d15 (125 d35 d45 d55 d56

 

 

d16 d26 d36 d46 d56 d66

135

The tensor D may also expressed as

D(rc, y. z) =

E1113

E1112

E1122

E2222
E2233
E2223
E2213

E2212

E3312

(6.2)

 

 

where d” (and EW) are functions Of the spatial coordinates 3:, y and 2. Assuming
that the material prOperties over each voxel are constant, one may relate each voxel
X (k, l, m) to a three-dimensional Haar scaling function ¢jk1m(x, y, z) in the following

representation Of the material tensor:

21—121— 121- 1

=2 2: Z diiiln 3’37 2: 9,2)

Ic=0 I: 0 m=0

d"( (:1: y,,z (6.3)

In the most general case, twenty-one expansions are required to approximate the
material tensor D. Depending on the level of resolution, storing the material tensor
expansions may require a large amount of computer memory. To help reduce these

requirements, we use a simplified material model. When using a simplified model, the

 

material tensor at any location is represented by a scaled version of fixed isotropic

136

material tensor D0, where D0 has the form

 

 

 

1 — V0 11° V0 0 0 0
V0 1 — 1x0 11° 0 0 0
0 E0 V0 11° 1 — V0 0 0 0 ( )
D = 6.4
1 + 11° 1 - 2 0
( ) ( V ) 0 0 0 132"" 0 0
0 0 0 0 132”” 0
0 0 0 0 0 113‘”
where E0 and 11° are the Young’s Modulus and Poisson’s ratio, respectively.
The material distribution in Q is expressed now as
1303.31.21) = 10(3, 9, Z)”D° (6-5)

with exponent p > 1. This reduces the representation of the material distribution to

a single expansion. The expansion of the design variable p is

21—1 21-121-1

M15131, Z) = Z: Z Z pjklm ?£?r($,ya Z) (66)

k=0 l=0 m=0
The discretization Of the displacement field is Of the form

21—1 21-1 21—1

“(mil/12) : Z Z: : ujklm¢;k[m($,y, Z) (67)

k=0 I=0 m=0

where uJ-kzm are the expansion coefficients and 34“", are some sufficiently smooth

137

wavelets. Based on the results Observed in two-dimensional topology Optimization,
linear spline S3 scaling functions are used to approximate the displacement field. In a

similar fashion, the weight functions v(:1:, y, z) are approximated using the expansion

21-121—121—1

V(.’B, 31,2) = Z Z Z vjklm¢jklm(x1yaz) (68)

k=0 I=0 m=0

Remark: The boundary integrals associated with the boundary and loading con-
ditions must be converted tO volume integrals. This procedure is a simple exten-
sion of the two-dimensional derivation. The essential boundary conditions associated
with Bu = 0 lead to enforcing that average displacement is zero over each voxel
individually. The natural boundary conditions associated with applied tractions are
converted into body forces applied over voxels in the one voxel thick boundary ap-

proximation Aw.

6.3.2 The Discretized Equations

Using the approximations Of the material distribution, displacement functions, bound-

ary conditions and the loading conditions, we develop the discretized problem

ABT U F

B 0 A 0

138

In three-dimensional problems, it is convenient to express the discretized equations

in the partitioned form

 

 

A11 A12 A13 111
A : A21 A22 A23 1 U = 112
_ A31 A32 A33 _ “3

/\1 f1

)‘2 and F = f2

A3 f3
(6.10)

where the u,- are vectors containing the wavelet expansion coefficients Of the displace-

ment field and A,- are vectors containing the wavelet expansion coefficients Of the

multipliers (only along the discretized boundary Of w).

Constructing stiffness matrix A

Entries in the stiffness matrix A are again a material weighted sums of connection

coefficients. Using the approximations for D, u and v we find that

1,0,o,1,0,0)
(A11)(abc)(klm) Z dlzlwqrcipqrabcum ‘1'
p99,,-
__ 2 11010103130)
(A12)(abc)(klm) _ Z dqurcqurabcklm +
new

(1,0,0,0,0,1)

(A13)(abc)(klm) = Z dqurcqurabcklm ‘1' d

P1917

1,0,0,0,1,o
(A21)(abc)(klm) = 2d“ ijqrabcklm) +d1'2 C

J
P1917

1,0,0,1,0,0
(A22)(abc)(klm) = 2 dis ipqrabcum) +61!”

p.q.r
_ 23 (0.1.0.0.0.1)
(A23)(abc)(klm) — Z dquerpqrabcklm +
par

139

(0,1,0,0,l,0) 55 (0,0,1,0,0,1)
qurabcklm + d C

Jpqr qurabcklm

0. .0. .0.0
(.5314...) (6.11)

55 C(oioal 111010)
'pqr qurabcklm

(0,1,0,1,0,0)
'pqrabcklm

jgfigflfif,’ + 53,060,501» (6.12)

qurabcklm

(0.0.1.0.110)
pqrabclclm

_ 55 (1,,,,,00001) 13 (0,,01,1,,00)
(A31)(abc)(klm) — 2 dj pququrabclclm + dquerpqrabclclm

p’q’r
0W10001 0W01010
(A32)(abc)(klm) : qur ijqrabclclm)+ dquerpqrabclclm) (613)

pm.

55 (1,,,,,+00100) 4 (0,,,,,10010) 33 (0,,,,,01001)

(A33)(ab6)(k1m) 2 dj pququrabcIclm+ J'pquJ'pqrabcklm + deququrabclclm
pIq’r

where

0.1353319 = (/a Tom: ,.(4>D<4:.(x)49.) x
([4 413115054 4,.,(y)D"4 ”(4)442 )x (6.14)
(I... 4513:1401 144094 1....(z)dn)
and Di denotes the derivative operator.
_Ile_m_ar_k: The structure Of the stiffness matrix A is similar to that found in two
dimensions: it is symmetric and sparse.

Remark: The stiffness matrix A is only semi-positive definitive and has three zero

eigenvalues with eigenvectors that correspond tO the rigid body modes

 

 

T
(1 1 0 0 0 0)

p = 0 0 1 1 0 0 (6.15)
\0 0 00 1 1)

Constructing stiffness matrix A),

Here, we examine the stiffness matrix associated with the case where the domain 9

contains an isotrOpic, homogeneous material. The homogeneous stiffness matrix, Ah,

140

is partitioned in a similar manner as the heterogeneous stiffness matrix, A, namely

)1 h h
A11 A12 A13

Ah 2 A33 A32 A33

 

 

h h h
A‘31 A32 A33

Entries in A), are of the form

where

(A’fi)(abc)(klm)
(A112) (abc)(klm)

(A113) (abc)(klm)

(A31) (abc)(klm)
(A512 )(abc)(klm)

(A33 )(abc)(klm)

(A91) (abc)(klm)
(Ag2)(a5c)(k1m)

(A33 )(abc)(klm)

(a’flI’Y’C’n’o)
Cjabclclm

1H001,1,0,0 OWIOOIO mW01001
dllcjabck ) + d660(al>c’kzflni ) + d55C jabcklm )
1W00010 0W10100
412C}.....: ’+ 4601.....4 ’

1,,001,0,,01 0,,,,,01100
(1136405)“ ) + d550<abeklm )

6 (1,0,0,o,1,0) 12 (0,1,o,1,0,0)
‘16 Cjabeklm +d Cjabcklm

d560(100100)+d220(010010)+d44C(001001)

3 (0,1,0,0,0,1) 44 (0,0,1,0,1,0)
d2 Cjabcklm + d Cjabcklm

55 (11010101011) 13 (01011111010)
d Cjabcklm + d Cjabcklm

44 (0.1.0.0.0.1) 3 (0,0,1,0,1,0)
d Cjabeklm +512 Cjabcklm

1,001,160 o,1,0,,,010 0,0,1,o,0,1
(155021” ) + (1440(036“ (m ) + d330]amzm )

(f: 004;,(x)o< ;,(x)dn,) x

(f 131’; We" ,(y)do,)x

(/ 01 new" 1.,.z.()49)

(6.16)

(6.17)

(6.18)

(6.19)

(6.20)

Recall that in the two—dimensional case, the sub-matrices Of the homogeneous

141

matrix were block circulant. In three dimensions, the partitions of the homogeneous
stiffness matrix Ail are two—block circulant matrices. This means that each matrix At,
is block circulant and its block matrices are themselves block circulant. Fortunately,
two-block circulant matrices are also diagonal when represented in the Fourier space
and a simple three-dimensional Fourier transform can diagonalize the sub-matrices of
Ah.

Remark: The homogeneous stiffness matrix Ah is only semi-positive definite and

has three zero eigenvalues associated with the rigid body modes shown in (6.15).

Constructing the constraint matrix

The constraint matrix is associated with constraining the average displacement over

a given voxel. As in the two-dimensional case, we have

B<abc).(klm) = [Au 5‘33! E’kzmdfl (6-21)

Constructing the load vector

After the tractions are approximated by body forces, the load vector is constructed

based on the Haar discretization of the body force (f,)jk,m.

21—121-12j-1

00156):ch Z :0 (f') jklm Cjklmabc (6.22)

=Ol=0 =0

where

¢haar u
Cjklmabc :/¢ jklm jabcdg-2 (6.23)

142

6.4 Post Processing Three-Dimensional Solutions

The solution to (6.9) is a fixed-scaled wavelet expansion of the displacement field. The
wavelet expansion is post processed to determine the average displacement, strain,
and stress over each voxel in the domain w. The average displacement over each voxel
X (a, b, c) is given by

1

(”0““) V(a, b, c) X(a,b,c)

1
ui(x1y: Z)dQ = 5:31-- A X(a) b: C)’U5(.’L', ya Z)dfl

21—12i—121—1

= 23j/2 Z Z Z lIlium/I ?$($1 y, z)¢_t1"klm($a y: Z)dQ (6'24)

k=0 1:0 m=0 9

2i-1 2L1 2i-1

23"” Z Z Z mjk.m0,‘-3L‘lf3m

k=0 l=0 m=0

where V(a, b, c) = 2‘3”2 is the volume of voxel X (a, b, c) and x(a, b, c) is an indicator
function corresponding to the volume occupied by the voxel X (a, b, c). The terms

01(3ng are again connection coefficients of the form

0135172,, = (fa: ¢92M<z>D° yummy) x
(/n ¢;:°'(y)Dfi¢;-‘.(y)dny) x (6.25)
(/n 9:“'(z)DV¢;-‘m(z)dnz)

The strain field is computed based on the following strain / displacement relationship.

 

 

_ 1 Bu, an
Eij — 5 (6133' + 6112;) (6.26)

143

The average strain field for a given voxel X (a, b, c) is given by

21— 121— 121— 1

_ 1,,00
(€11)(a11c) = 23172 2 Z 2 “111.11,. 03(11me

k=01=0 m=0
21—121—121—1

_ - 0,10
(€22)(abc) = 231/2 2 Z: Z “211.111.03‘a51ck3n1

k=0 (:0 m=0
2j—12j—l2j—l

_ 0,0,1
(€33)(abc) : 23j/2 Z Z Z u3jk1mC§abck2nl

k: 0 1:0 m==0
21— 121— 121— l (6'27)

(€23)(abc) : 2j :0 g :0 (“23111111 Cégbgklbu'b u3jk1m0;2b::r)ru)
_ m_

21— 121— 121—1

(513)(abc) = 2j Z Z Z (uljkzmciobgklbu+u31'kthJg100) )

k_=01=0m=0
21— 121-121—1

(5.12)(abc) : 2i 2 2: Zn

k=0 1:0 m=0

C(01110) (11010)
“1,1111. Cjabckml + 112111ij11111;ka

where connection coefficients 0;;ka are defined in (6.25).
Remark: Expressions for the average stress field are based on the average strain

field and follow the stress/strain relationship: 6 2 DE.

6.5 Performance Of The Fictitious Domain Ap-

proach

This section examines the performance of the fictitious domain approach to solve elas-
ticity problems associated with three-dimensional topology optimization. A series of
‘ elastostatic analysis problems were solved to determine which values of the fictitious
material provide reasonable approximation of the solution in the domain of interest
w. From these investigations, the fictitious domain corresponds to an isotropic ma-

terial with E2 = 10’5E1 and V2 2 121, where the strongest possible design material

144

 

 

Design

Domain

 

 

 

 

 

 

 

A
V

Figure 6.2: Example problem statement

corresponds to an isotropic material associated with E1 and V2.

In this investigation, we use a simple extension of the PCG solver used in two
dimensions to solve the discrete system of equations (6.9). The performance of the
PCG solver is again determined by the number of iterations required for the PCG
solver to converge. Here we are interested in how the performance PCG solver changes

for (1) different levels of resolution and (2) varying amounts of fictitious material.

6.5.1 Effect Of Varying Resolution

A topology optimization problem with design domain as displayed in Figure 6.2
and volume constraint V = 0.25 is solved with wavelet-based fictitious domain dis-
cretizations associated with three resolution levels: j = 4, j = 5 and j = 6. The

resulting degrees of freedom in these discretizations are displayed in Table 6.1 (The

145

domain dimensions in Table 6.1 correspond to the discretizations of the design do-

main in voxels.) The simplified material model with p = 2.5 is used in the topology

Table 6.1: Discretization levels for analysis of PCG performance

 

 

 

 

 

 

Level Domain Dimensions Design Voxels Total Voxels Total DOF
4 12 x 6 x 6 432 4,096 12,288
5 24 x 12 x 12 3,456 32,768 98,304
6 48 x 24 x 24 27,648 262,144 786,432

 

 

 

 

 

 

 

 

 

optimization problem for each of the discretizations listed in Table 6.1. The optimal
topologies found for these discretizations are displayed in Figure 6.3. We note from
the optimal topologies that the results appear to be converging to the same structure
represented at different resolutions. Recall that the material description and displace-
ment field were approximated at the same discretization level. Although there exist
checkerboards in the level j = 4 solution, no checkerboards exists in the level j = 5
and level j = 6 solutions. Unlike the checkerboards observed in the two-dimensional
problem when the material distribution and displacement field were approximated at
the same level in two-dimensional problems, the checkerboards in the level 4 solution
are due to the extremely coarse discretization of the design domain used. As the
discretization is refined (level j = 5 and j = 6), the checkerboards disappear. The
absence of checkerboards is beneficial because we no longer have to approximate the
material distribution at a lower level than the displacement field, which results in a
loss in the detail of the shape description.

The performance of the PCG solver during the topology optimization sequence is
shown in Figure 6.4. Although the total number of degrees of freedom varies greatly

in the three models, the performance of the PCG solver is relatively consistent among

146

 

(a) Level 1' :4 (b) Level 1' = 5

 

(c) Level j = 6

Figure 6.3: Optimal topologies found for various discretizations

147

80

 

To+ r I r l I
Level4 0
70 ' +CI Leve15 + ‘
60 Cl Level6 Cl

50- SE _

<>
PCG 0 80

 

 

 

Iterations 4 — D + 0 d
30 - DB

20 L 8&89 D 0 ++ L

10 - [3+ 0 +$ o 13$ 33:] -

++ + 0 1:10 1:1

0 1 1 1 1 + 199%

0 5 10 15 20 25 30

Topology Optimization Iteration

Figure 6.4: PCG convergence properties for various discretizations

them. This suggests level insensitive convergence in the solution method.

Another way to compare the effect of discretization level on convergence is to
examine the PCG residual during the initial iteration of the topology optimization
sequence. Here the analysis is performed on the same initial shape. The residual
for the three models during the initial topology optimization iteration is graphed in
Figure 6.5. (For this investigation, the PCG convergence tolerance for the residual
[R, R] / [F, F] is 0.05.) Although the three models have very different numbers of de-
grees of freedom the PCG performance for the three discretizations is nearly identical.

This again verifies that the PCG algorithms insensitive to the discretization level.

148

 

8 Level 4 o
Level 5 +
@ Level 6 Cl

I‘ll

0.1

E0
:u

'5
.3

0.01

 

 

0.001 ‘ ' '
0 5 10 15 20
PCG Iteration

 

Figure 6.5: PCG convergence for the initial topology Optimization iteration

6.5.2 Effect Of Varying The Amount Of Fictitious Material

Here, we examine the effect of varying the amount of fictitious material on the per-
formance of the PCG solver. We examine various discretizations the of cube design
domain shown in Figure 6.6. Cubes of various sizes are embedded into a fictitious
domain at discretization level j = 5. We consider four discretizations of the design
domain with voxel dimensions of 30 x 30 x 30, 28 x 28 x 28, 24 x 24 x 24 and 16 x 16 x 16,
each centered in a 32 x 32 x 32 fictitious domain. The Optimal topologies found at
these discretizations are shown in Figure 6.7. From the results, we note that the solu-
tions are converging to similar structures represented in varying detail. The amount
of fictitious material does not appear to affect the optimal topology determined. We
next examine the efficiency of PCG solver for these problems. The PCG solver itera-

tion history during the topology optimization sequence for each problem is displayed

149

 

 

Design 1

Domain

 

 

 

 

 

 

 

P

Figure 6.6: Problem statement for fictitious material investigation

in Figure 6.8. Unlike the results seen in two dimensions, there appears to be a corre-
lation between the amount of fictitious domain and the PCG solver convergence: as

the amount of fictitious domain decreases, the number of PCG iterations increases.

6.6 A More Memory Efficient Solution Technique

In moving from two-dimensional problems to three-dimensional problems, the com-
putational resources required to solve a level j problem become approximately 2j
times larger; in two dimensions, the level j problem size is 2N 2 x 2N 2, while in three
dimensions the level j problem size is 3N 3 x 3N 3 (again N = 23'). The problem size
may significantly increase the computational time and memory requirements. The

wavelet-based solution technique helps to reduce the computational time required to

150

 

33

 

(a) Domain 30 x 30 x 30 (b) Domain 28 x 28 x 28
(c) Domain 24 x 24 x 24 (d) Domain 16 x 16 x 16

Figure 6.7: Results of topology optimization for various cube discretizations

151

 

 

 

 

120 I I I l I T
00 30x30x30 O
100 _ 28x28x28 + _
++ 24x24x24 D
C] 16x16x16 X
80 - 9x -
PCG
Iterations 60 _ éfiooo<>0¢0 +00 ‘
40 ' ES ++ mes -
o sé D 9: ee one
a £8 M08 “ago
20 fi 5x5 ano‘
X x X X +++*+m++-§-
0 . 1 . XXXXQXX x x
0 5 10 15 20 25 30

Topology Optimization Iteration

Figure 6.8: PCG solver performance for various cube discretizations

solve these problems iteratively, but memory requirements are still significant with
increasing problem size. In our approach, a substantial portion of the main mem-
ory requirements are related to the data vectors associated with the PCG algorithm.
The use of the fictitious domain greatly increases memory requirements: adding the
fictitious domain adds more degrees of freedom (more unknowns) to our system of
equations. Depending on the shape of the design domain w, the number of unknowns
in the discrete system of equations associated with degrees of freedom in the ficti-
tious domain could easily dominate the size of the problem. This trend was observed
in two-dimensional examples and is even more apparent in three dimensions. To
overcome this problem, a reduced system of equations is solved. The reduction is
achieved by setting a significant portion of displacement coefficients associated with

the fictitious domain to zero. This is implemented in the following manner:

152

1. Displacement coefficients associated with scaling functions whose support is

entirely within the fictitious domain, 9 \ w, are set to zero.

2. Appropriate terms in the unknown vector U are eliminated, reducing the system

of equations.

This approach reduces to the total number of degrees of freedom to be approximately
equal to the number of degrees of freedom associated with the domain of interest w.
However, the preconditioner must still be applied on vectors of dimension N 3 because

inversion of the preconditioner requires Fourier analysis on a vector of size N 3.

6.7 Three-Dimensional Topology Optimization

Examples

The examples in this section use the memory saving techniques outlined in the pre-
vious section. A simplified material model with material exponent p = 2.5 is used.
The isotropic reference material tensor D0 is built with E0 = 1.0 and 11° = 0.3. The
fictitious material is set to an isotropic material with E = 10‘5 and u = 0.3. Also, the
material discretization using fixed-scale Haar expansions and the displacement dis-
cretization using linear spline S3 fixed-scale expansions are represented at the same
level. Finally, the PCG convergence tolerance for the residual [R, R] / [F, F] is set to

10’3

153

 

 

D/D ‘ 1.

Design 5

Domain

1
f

 

 

 

 

 

 

 

V

8

Figure 6.9: Problem statement for truss example 1

6.7 .1 Single Load Case Examples
Truss example 1

Figure 6.9 illustrates the setting of truss example 1. The problem is solved at dis-
cretization level j = 6 with a volume constraint of V = 0.35. This discretization
level results in a design domain with dimensions (in voxels) of 48 x 30 x 12 leading to
17,280 design voxels. The optimal topology and performance of the PCG solver are
shown in Figure 6.10

The solution method generates a highly detailed symmetric optimal topology
without any checkerboards. We also note that the PCG solver follows the typical
performance characteristics of the problems presented here, requiring the most PCG
iterations in the first few topology optimization iterations and requiring less as the

154

 

(a) Optimal Topology

 

 

 

 

 

 

 

 

 

 

140
120
100
PCG 80
Iterations 60
40
20
0 l I l
0 10 15 20 25 30
Topology Optimization Iteration
(b) Performance of PCG solver
3:9”? l— I I I I I
1800 - -l
1500 ~ -
Compliance 1200 F -
900 - -
600 - -
303 r-l l l I l l -
0 25 30

10 15 20
Topology Optimization Iteration
(c) Compliance history
Figure 6.10: Results of topology optimization for truss example 1

.155

 

 

l>/> 1

Design

Domain

v2 2

L
‘

 

 

 

 

 

 

 

V

P

Figure 6.11: Problem statement for truss example 2

optimization sequence continues.

Truss Example 2

This problem has the same design domain discretization as the last example, but with
different loading conditions. The problem statement is illustrated by Figure 6.11.
This problem is solved with a volume constraint of V = 0.35. The optimal topology
and performance of the PCG solver are shown in Figure 6.12 As seen in Figure 6.12,
the wavelet-based topology optimization again generates highly detailed resolutions

with very efficient PCG solver performance over the optimization sequence.

156

 

(a) Optimal Topology

 

E
4

   

PCG 80
Iterations 60
40

  

IIIII
IIIIII

 

 

 

 

 

 

 

 

20
0
0 10 15 20 25 30
Topology Optimization Iteration
(b) Performance of PCG solver
3000 I l I I l l
2500 - -
2000 - ..
Compliance 1500 - -
1000 - -
500 - —
0 I I I l I l
0 5 10 15 20 25 30

Topology Optimization Iteration
(c) Compliance history
Figure 6.12: Results of topology optimization for truss example 2

157

6.7 .2 Multiple Load Case Examples

This section discusses two multiple load case examples. Recall that in multiple load
case examples, the loads are applied separately and the PCG solver determines the

displacement field associated with each load case in parallel.

X—truss example

Figure 6.13 displays the problem for the X—truss example. The problem is solved for
three different volume constraints: V = 0.20, V = 0.15 and V = 0.10. A discretization
level j = 6 is used and represents the design domain with voxel dimensions 58x 26 x 26,
which generates 43,904 design voxels. The optimal topologies associated with the
three volume constraints are shown in Figure 6.14. As expected, the Optimal
topologies become “slimmer” as the volume constraint decreases. The performance
Of the PCG solver for these problems is shown in Figure 6.15. For each volume
studied, the performance of the PCG solver is as in the previous examples. However,
when comparing the PCG performances for the different volumes tO one another we
notice that as the volume percentage decreases, the number of iterations in the PCG
solver increases. This suggests that the fictitious domain approach may not be well

suited for topology Optimization problems with very small volume constraints.

Bridge example

The last multiple load case problem is shown in Figure 6.16, where an inner rect-
angular prism is set as non—design material. This is achieved by denoting the inner

rectangular prism as part Of the fictitious domain. This problem is also discretized at

158

 

 

 

 

 

 

 

 

 

Figure 6.13: Problem statement for the X-truss example

159

 

(a) Optimal topology for V = 0.20

 

(b) Optimal topology for V = 0.15

 

(c) Optimal topology for V = 0.10

Figure 6.14: Results Of topology optimization of the X-truss example for various
volume constraints

160

160
140
120

PCG 100
Iterations 80

40
20

160
140
120

PCG
Iterations

ss§

20

160

140

120

PCG 100
Iterations 80

40
20

 

 

 

 

 

 

 

 

 

0 5 10 15 20 25

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

30
Topology Optimization Iteration
(a) Performance of PCG solver for V = 0.15
r I I I I I
p-
1 I 7 7 l 1 1
0 5 10 15 20 25 30
Topology Optimization Iteration
(b) Performance of PCG solver for V = 0.15
I I I I I I
F r I
'- 'l
IT I I I I I I I I
0 5 10 15 20 25 30

Topology Optimization Iteration

(c) Performance of PCG solver for V = 0.15

Figure 6.15: PCG performance for various volume constraints for the X-truss example

161

Design Non-Design

Domain Domain

 

 

 

 

 

 

 

._ __. _____ ?____..______

 

 

 

 

 

 

 

P P P

Figure 6.16: Problem statement for the bridge example

level j = 6. At this discretization, the large domain as voxel dimensions 56 x 28 x 28
and the inner non-design domain as voxel dimensions 56 x 12 x 12. This configuration
leads to 35,840 design voxels. Figure 6.17 displays the results of the topology opti-
mization for this example. Again, we see highly detailed results with efficient solution

histories.

6.8 Summary

The fictitious domain approach using wavelet-Galerkin solution techniques was ex-

tended to three dimensions and used to solve three-dimensional topology optimization

162

 

(a) Optimal Topology

 

G)
O
r1

PCG igt
Iterations 30

IIllIII

 

 

0 5 10 15 20 25 30
Topology Optimization Iteration

 

O

(b) Performance of PCG solver

 

1400 l I l I I

1200

1000

Mean 800
Compliance 600
400

200

0

TIIIII
IIIIII

 

 

 

I I l J

5 10 15 20 25 30
Topology Optimization Iteration

 

o-

(c) Compliance history

Figure 6.17: Results Of the topology optimization for the multiple load case bridge
example

163

problems. The method produces highly detailed shapes with level insensitive conver-

gence rates.

164

Chapter 7

Conclusions

7 .1 Summary

Wavelet-based techniques were formulated and used to solve two- and three-
dimensional topology optimization problems. Below we list the significant conclusions

in the areas of PCG performance, checkerboard patterns, and Optimal topologies.
o PCG convergence

1. The wavelet-based fictitious domain approach combined with the PCG

solver exhibit level insensitive convergence rates.

2. The convergence rates Observed with PCG solver are substantial improve-

ment over iterative finite elements methods with diagonal preconditioners.

3. The PCG convergence rates appear to be dependent upon the wavelet used
to discretize the displacement field. Wavelets with shorter support result
in better convergence rates.

165

o Checkerboard patterns

1. When same level approximations are used for material and displacement
approximations, checkerboard patterns were Observed in two-dimensional
problems with Daubechies D6 and linear 83 spline displacement field ap-
proximations. These results are similar to traditional tOpology Optimiza-

tion results found with linear, four noded quadrilateral finite elements.

2. Under the same conditions, quadratic S4 Spline displacement field approx-

imations produce results that possess a “layered” material distribution.

3. The checkerboard and “layered” material distributions are caused by nu-
merical phenomena and vanish when the Haar material approximation is

performed one level lower than the displacement field approximation.

4. In three dimensions, the checkerboard-like phenomena were not observed

even when using the same level material and displacement approximations.
0 Optimal topologies

1. In all relevant examples, the wavelet-based Optimal topologies correlated
well with results found using traditional topology Optimization based on

finite element discretizations.

2. The extent Of the fictitious domain showed little effect on the Optimal
topologies, while the strength Of the weak fictitious material has great

importance.

166

3. The Optimal topologies did not exhibit any significant resolution dependent
results, as is Often seen in traditional topology Optimization using finite

elements.

7 .2 Areas Of Future Work

Based on the results Observed, the following areas should be investigated more

throughly.

o Quantify convergence behavior: Although the convergence properties Of the
wavelet-based solution technique produce accurate results very efficiently, a

more complete investigation Of the convergence characteristics is needed.

0 Multi-scale approximations: All the work presented here is based on fixed-scale
expansions for the material distribution and displacement field. Multi-scale
expansions Of these variables and with appropriate compression techniques may

greatly improve the method.

0 Eigenvalue analysis: Determine how the wavelet-based methods may benefit
eigenvalue analyses and incorporate with topology Optimization problems that

have frequency based Objective functions and constraints.

167

BIBLIOGRAPHY

168

Bibliography

[1] O. Axelson. Iterative Solution Methods. Cambridge University Press, 1996.

[2] N. S. Bakhvalov and A. V. Knyazev. Fictitious domain methods and compu-
tations of homogenized properties of composites with a periodic structure Of
essentially different components. In G. I. Marchuk, editor, Numerical Methods
and Applications, pages 221-266. CRC Press, 1994.

[3] M. Bendsoe and N. Kikuchi. Generating Optimal topologies in structural design
using a homogenization method. Computational Methods in Applied Mechanics
and Engineering, 71:197—224, 1988.

[4] J. Benedetto and J. Frazier, editors. Wavelets: Mathematics and Applications.
CRC Press, 1994.

[5] J. H. Bramble and J. E. Pasciak. A preconditioning technique for indefinite
systems resulting from mixed approximations for elliptic problems. Math. Comp.,
50:1—17, 1988.

[6] A. Tabacco C. Canuto and K. Urban. The wavelet element method, part i:
construction and analysis. Preprint NO. 1038, Istituto del Analisi Numerica del
C.N.R. Pavia, to appear in ACHA, 1997.

[7] A. Tabacco C. Canuto and K. Urban. The wavelet element method, part ii:
realization and additional features in 2d and 3d. Preprint NO. 1052, Istituto del
Analisi Numerica del C.N.R. Pavia, to appear in ACHA, 1997.

[8] C. Chao. A remark on symmetric circulant matrices. Linear Algebra and its
Applications, 103:133—148, 1988.

[9] C. K. Chui. Wavelets: A Mathematical Tool for Signal Analysis. SIAM, 1997.

[10] S. Dahlke and I. Weinreich. Wave]et-galerkin-methods: an adapted biorthogonal
wavelet basis. Technical Report A-91-25, Freie Univ. Berlin, 1991.

[11] W. Dahmen and C. A. Micchelli. Using the refinement equation for evaluating
integrals Of wavelets. Siam J. Numer. Anal, 30:507—537, 1993.

[12] W. Dahmen and R. Schneider. Composite wavelet bases for Operator equations.
Perprint TU Chemnitz, 1996.

169

[13] I. Daubechies. Orthonormal bases of compactly supported wavelets. Comp. Pure
Appl. Math., 41:909-996, 1988.

[14] A. Diaz and O. Sigmund. Checkerboard patterns in layout optimization. Struc-
tural Optimization, 10:40—45, 1995.

[15] H. Federer. Geometric Measure Theory. Sprin—Verlag, 1969.

[16] R. 0. Wells Jr. and X. Zhou. Wavelet solutions for the dirichlet problem. Tech-
nical Report Computational Mathematics Laboratory TR92-02, Rice University,
1992.

[17] S. Qian K. Amaratunga, J. Williams and J. Weiss. Wavelet-galerkin solutions
for one-dimensional partial differential equations. Int. J. Numer. Meth. Eng,
37:2703—2716, 1994.

[18] R. V. Kohn and G. Strang. Optimal design and relaxation of variational prob-
lems, 1. Communications on Pure and Applied Mathematics, 39(1):113—137, 1986.

[19] R. V. Kohn and G. Strang. Optimal design and relaxation of variational prob—
lems, II. Communications on Pure and Applied Mathematics, 39(2):139—182,
1986.

[20] R. V. Kohn and G. Strang. Optimal design and relaxation Of variational prob-
lems, III. Communications on Pure and Applied Mathematics, 39(3):353—377,
1986.

[21] A. Kunoth. Multilevel preconditioning - appending boundary conditions by
largrange multipliers. Advances in Computational Mathematics, 4:145—170, 1995.

[22] A. R. Diaz M. P. Bendsoe and N. Kikuchi. Topology and generalized layout
Optimization of elastic structures. In M. P. Bendsoe and C. A. Mota Soares,
editors, Topology Design of Structures, pages 159—205. DordrechtzKluwer, 1992.

[23] Y. Maday, editor. The evaluation of connection coefi‘icients of compactly sup-
ported wavelets, Princeton University, 1992.

[24] S. Mallat. Multiresolution approximation and wavelet orthonormal bases of l2.
Trans. Amer. Math. Soc., 315:69—88, 1989.

[25] H. P. Mlejnik and R. Schirrmacher. An engineering approach to optimal material
distribution and shape finding. Comp. Meth. Appl. Mech. Engng, 106:1—26, 1993.

[26] F. Murat and L. Tartar. Calcul des variations et homogeneisation. Coll de al
Dir. des Etudes et Recherches de Electricite de France, pages 319—370, 1985.

[27] S. Qian and J. Weiss. Wavelets and the numerical solutionof partial differential
equations. J. Comput. Phys, 106:155—175, 1993.

170

[28] R. O. Wells Jr. R. Glowinski, A. Rieder and X. Zhou. A wavelet multi-grid pre-
conditioner for dirichlet boundary value problems in general domains. Technical
Report Computational Mathematics Laboratory TR93-06, Rice University, 1993.

[29] R. 0. Wells Jr. R. Glowinski, T. W. Pan and X. Zhou. Wavelet and finite el-
ement solutions for the neumann problem using fictitious domains. Technical
Report Computational Mathematics Laboratory Technical Report, Rice Univer-
sity, 1994.

[30] K. Suzuki and N. Kikuchi. Shape and topology Optimization for generalized
layout problems using the homogenization method. Comp. Meth. Appl. Mechs.
Engng., 93:291—318, 1991.

[31] K. Suzuki and N. Kikuchi. Shape and topology Optimization in three-dimensional
solid-structures. Technical report, University of Michigan, USA, 1991.

[32] W. Sweldens. Construction and Applications of Wavelets in Numerical Analysis.
PhD thesis, Katholieke Universiteit Levuen, 1994.

[33] M. Ueda and S. Lodha. Wavelets: An elementary introduction and examples.
Technical Report UCSU-CRL 9417, University Of California, Santa Cruz, 1995.

[34] A. Kunoth W. Dahmen and K. Urban. A wavelet-galerkin method for the stokes-
equation. Computin, 56:259—302, 1996.

[35] M. V. Wickrhauser. Adapted Wavelet Analysis form Theory to Software. A K
Peters, 1994.

[36] J. C. Xu and W. C. Shann. Galerkin-wavelet methods for two-point boundary
value problems. Numer. Math, 63: 123—142, 1992.

171

nICHIGAN smIE UNIV. LIBRARIES
[Ill[Ill]lllllllllllllllllll[ll][lllllllllllllll
31293017743968