PLACE IN RETURN BOX to remove this checkout from your record.
To AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 cJCIRC/DatoDue.p65-p.15

COMPUTATIONAL DESIGN OF MECHANICAL
STRUCTURES IN ELASTICITY USING
MULTI-RESOLUTION ANALYSIS

By

S udarsanam Che llappa

A DISSERTATION
Submitted to
Michigan State University
in partial fiilfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Department of Mechanical Engineering

2003

 

C OMPL’TATI

A form
elastic system
various resol
homogenizatit
resolution an:
elasticity OPE
information
mathces that
models of la
in‘v’OMng CO
intensive CO,
is applied u
heterogmem
mmPIES in ‘

presented.

ABSTRACT

COMPUTATIONAL DESIGN OF MECHANICAL STRUCTURES IN ELASTICITY
USING MULTI-RESOLUTION ANALYSIS
By

Sudarswram Chellappa

A formal methodology for reducing the size of models used in the analysis of
elastic systems is presented. This involves an explicit representation of the model at
various resolutions and is accomplished using a projection generated numerical
homogenization procedure. The fiamework for this analysis is derived from the multi-
resolution analysis associated with the construction of wavelet bases. This is applied to
elasticity operators to average fine scale properties and behavior while limiting loss of
information. In discretized form, the method produces equivalent smaller stiffness
matrices that can be used as building blocks (super-elements) to construct reduced
models of larger systems. The principal application envisioned is in design problems
involving complex structural systems, such as in crash-worthiness design, where very
intensive computations demand computational efficiency. This model reduction scheme
is applied to the problem of layout optimization of structures involving multi-scale
heterogeneities of sizes that may be comparable to the size of the structure. Numerical
examples in which the heterogeneities are in the form of perforations of various sizes are

presented.

T 0 my parents

iii

l woult
Your guidar
Dr. Andre E
time spent i
like to thanl

Lastly.

motivation '

ACKNOWLEDGEMENTS

I would like to express my Sincere gratitude to my advisor Dr. Alejandro Diaz.
Your guidance and support has been invaluable. I thank the members of my committee:
Dr. Andre Benard, Dr. Michael Frazier and Dr. Farhang Pourboghrat. Thank you for the
time spent in reviewing this document and providing valuable suggestions. I would also
like to thank Dr. Martin Bendsrae. Thank you for your insight and suggestions.

Lastly, I would like to thank my family. Thank you for all your support and

motivation without which this endeavor would not have been possible.

iv

LIST 1
LIST (
1 Intro

2 Wave

3.4.:

3.4 3
35 Comp

 

TABLE OF CONTENTS

LIST OF TABLES viii
LIST OF FIGURES ix
1 Introduction ........................................................................................................ l
2 Wavelets And Multi-Scalc Representation Of Functions ................................. 7
2.1 Introduction ................................................................................................. 7
2.2 Mum-Resolution Analysis ............................................................................ 8
2.3 Example Using the Haar Basis ...................................................................... 10
2.4 Generalized Orthogonal Scaling Functions and Wavelets .............................. 10
2.5 The Discrete Wavelet Transform .................................................................. 15
2.6 Constructing The Basic Scaling Function and Wavelet ................................. 19
2.7 Function Approximation Using Orthogonal Projection ................................. 23
2.8 Wavelets In Multiple Dimensions ................................................................. 25

3 Model Reduction In Elastostatics ...................................................................... 28
3.1 Scales And Material Properties In Elasticity ................................................. 29
3.2 The Fine-Scale Elasticity Problem ................................................................ 32

3.3 A Multi-Resolution Analysis Of Material Distribution ................................... 34

3 .4 A Multi-Resolution Analysis Of Displacement .............................................. 37
3.4.1 A Periodic Multi-Resolution Reduction Scheme ................................. 37

3.4.2 Transformation To Finite Element F orrn ............................................ 44

3.4.2.1 Method I ............................................................................... 45

3.4.2.2 Method H ............................................................................. 46

3.4.2.3 Method III ............................................................................ 48

3.4.3 Computational Aspects ...................................................................... 50

3.5 Comparison Of Model Reduction Schemes ................................................... 53

DJ

a
.3 .

 

3.5.1 Comparison Between MRA Of Material Distribution And MRA Of

Displacements ................................................................................... 53
3.5.2 Comparison Between MRA Techniques And Classical
Homogenization ................................................................................ 56

3 .6 Computing The Non-Periodic Reduced Stiffness Matrix Of A

Substructure ................................................................................................ 62

3 .7 Computing Fine-Scale Stresses By Augmentation ........................................ 66
3.7 .1 Computing The Periodic Large-Scale Nodal Displacements ............... 68
3.7.2 Conversion To A Wavelet Basis ......................................................... 71
3.7.3 Periodic Multi-Resolution Refinement ................................................ 72

3.8 Numerical Examples .................................................................................... 74
3.8.1 Example 1 ......................................................................................... 74
3.8.2 Example 2 ......................................................................................... 81
3.8.3 Example 3 ......................................................................................... 87

4 Multi-Scale Layout Optimization Of Structures ............................................... 93
4.1 Background: Topology Optimization Of Structures ..................................... 94
4.2 Optimal Design Using Macro-Scale Heterogeneities .................................... 97
4.2.1 Building A Model Using Reduced Substructures ................................ 99
4.2.2 Constructing Libraries Of Perforated Substructures ........................... 101
4.2.3 Sensitivity Analysis ............................................................................ 103

4.3 Compliance Minimization Problems .............................................................. 105
4.4 Optimization Using Fixed Layouts Of Reduced Substructures ...................... 108
4.5 Optimization Using Variable Layouts Of Reduced Substructures .................. 110
4.5.1 Perforated Substructures And Layout Dependency ............................. 110
4.5.2 A Dividing Approach To Optimization With Variable Layouts ............ 113
4.5.3 A Merging Approach To Optimization With Variable Layouts ............ 115

4.6 Examples ..................................................................................................... 119
4.6.1 Example 1 .......................................................................................... 119
4.6.2 Example 2 .......................................................................................... 121
4.6.3 Example 3 .......................................................................................... 122

5 Conc

L"
n—

.U‘
k)

lil’t
b.)

BIBLIO'

 

4.6.4 Example 4 .......................................................................................... 125

4.6.5 Example 5 .......................................................................................... 132

5 Concluding Remarks .......................................................................................... 139
5.1 Summary ...................................................................................................... 139
5.2 Conclusions ................................................................................................. 140
5.3 Areas Of Future Work ................................................................................. 141
BIBLIOGRAPHY .................................................................................................. 142

ME

Me

A]...
a).

List Of Tables

3.1 Main Results For Example 1 .......................................................................... 79
3.2 Main Results For Example 2 .......................................................................... 86
3.3 Main Results For Example 3 .......................................................................... 92

viii

3_1 Ha:
22(3) Ori
22th) Apr
22(0) Apr
2.2(d) Der
2.3 Sin;
2.4 The
2.5 The
2.6 Sing
27 The
28 Son
29 Ortl
210 Am
3 1 A st
3-2 Sun.
33 Basi:
3'4 Com
35 Com
3‘6 Com
3'7 Choir
3'8 Boun

using

List Of Figures

2.1 Haar scaling firnction and wavelet ................................................................. 11
2.2(a) Original function ........................................................................................... 12
2.2(b) Approximation at scale 6 ............................................................................... 12
2.2(c) Approximation at scale 5 ............................................................................... 12
2.2(d) Detail at scale 5 ............................................................................................. 12
2.3 Single stage discrete wavelet decomposition .................................................. 17
2.4 The convolution operation ............................................................................. 17
2.5 The dyadic down sampling operation ............................................................. 18
2.6 Single stage discrete wavelet reconstruction .................................................. 18
2.7 The dyadic up sampling operation ................................................................. 19
2.8 Some commonly used orthogonal scaling functions and wavelets ................... 22
2.9 Orthogonal projections fiom wavelet to finite-element spaces ........................ 23
2.10 A multi-scale representation of two-dimensional functions ............................. 27
3.1 A structure built from substructures of difl‘erent scales ................................... 31
3.2 Structure of L J and HJ_k .......................................................................... 43
3.3 Basis functions in wavelet and finite-element spaces ...................................... 44
3 .4 Conversion from wavelet to finite-element spaces using method I .................. 46
3.5 Conversion from wavelet to finite-element spaces using method II ................. 47
3 .6 Conversion from wavelet to finite-element spaces using method [[1 ............... 49
3.7 Choices of the material distributions for the comparison ................................ 55

3 .8 Bounds on the maximum relative error between the strain energy computed

using material MRA and displacement MRA ................................................. 55
3.9(a) Strain energy comparisons for material distribution (a) .................................. 60
3.9(b) Strain energy comparisons for material distribution (b) .................................. 60
3.9(c) Strain energy comparisons for material distribution (c) .................................. 61
3.9(d) Strain energy comparisons for material distribution ((1) .................................. 61
3.10 Substructure surrounded by a weak fictitious domain .................................... 63
3.11 Constant pre-strains applied to a patch of substructures ................................. 66

3.12

3.13
3.14

3.15
3.16
3.17
3.18

3.19
3.20

3.21
3.22
3.23
3.24
3.25
3.26

3.27
3.28

3.29
3.30
3.31
3.32
3.33
3.34

3.35

Decomposition of a general non-periodic fimction into a coarse, non-periodic

firnction and a periodic fimction ..................................................................... 69
Geometry and boundary conditions for example 1 ......................................... 74
Von Mises stress distribution using the fine-scale model of the structure

in example 1 .................................................................................................. 75
Assembly of substructures for example 1 ....................................................... 76
Von Mises stress distribution in reduced model using material MRA ............. 77
Von Mises stress distribution in reduced model using displacement MRA ...... 77
Von Mises stress in substructure (A) from the reduced model using

displacement MRA ........................................................................................ 78
Von Mises stress in substructure (A) after augmentation ............................... 78
Detail of Von Mises stress near substructure (A) obtained fiom the

fine-scale model ............................................................................................ 79
Geometry and boundary conditions for example 2 ......................................... 81
Von Mises stress distribution from the fine-scale model of the structure ........ 82
Arrangement of substructures for example 2 .................................................. 83
Von Mises stress distribution fi'om reduced model using material MRA ......... 84

Von Mises stress distribution from reduced model using displacement MRA .. 84
Von Mises stress in substructure (A) from reduced model using

displacement MRA ........................................................................................ 85
Von Mises stress in substructure (A) after augmentation ............................... 8S
Detail of the Von Mises stress in a region around substructure (A) fiom the

fine-scale model ............................................................................................ 86
Geometry and boundary conditions for example 3 ......................................... 87
Von Mises stress distribution fiom fine-scale model ....................................... 88
Arrangement of substructures for example 3 .................................................. 88
Von Mises stress distribution fiom reduced model using material MRA ......... 89

Von Mises stress distribution from reduced model using displacement MRA . 89
Von Mises stress in substructure (A) from reduced model using

displacement MRA ........................................................................................ 9O
Von Mises stress in substructure (A) after augmentation ............................... 90

3.36

4.1
4.2
4.3
4.4

4.5
4.6
4.7
4.8
4.9
4.10
4.11
4.12
4.13

4.14

4.15
4.16
4.17
4.18
4.19
4 20
4.2

4 22
4.23
4.24
4.25
4 76

‘5-

'C '1‘). '7‘ W.

’Y)

Ombnnmmi.

3.36

4.1
4.2
4.3
4.4

4.5
4.6
4.7
4.8
4.9
4.10
4.11
4.12
4.13

4.14

4.15
4.16
4.17
4.18
4.19
4.20
4.21
4.22
4.23
4.24
4.25
4.26

Detail of the Von Mises stress in substructure (A) obtained from the

fine-scale model ............................................................................................ 91
A structure with a composite microstructure ................................................. 96
A typical topology optimization problem and solution ................................... 97
A typical structure with macro-scale heterogeneities (perforations) ................ 98
A structure built from an assembly of substructures with circular

perforations ................................................................................................... 99
Assembly of substructures ............................................................................. 100
A non-conforming assembly of discretized substructures ............................... 101
Substructures with circular perforations for various finest levels .................... 103
Substructures with rectangular and circular perforations ................................ 105
Optimization using fixed layout of substructures ............................................ 107
Optimization using variable layouts of substructures ...................................... 107
A design domain and a possible layout of reduced substructures .................... 109
Boundary traction test ................................................................................... 1 10
Comparison of strain energies in substructures with different number of
perforations for various boundary tractions .................................................... l 1 1
Comparison of strain energies in substructures with different number of
perforations for various constant pre-strains .................................................. 112
Illustration of possible evolution of layouts using the dividing approach ......... 113
Illustration of possible evolution of layouts using the merging approach ........ 116
Problem description for example 1 ................................................................. 120
Optimal layout for example 1 ........................................................................ 120
Problem description for example 2 ................................................................. 121
Optimal layout using multiple load cases for example 2 .................................. 121
Optimal layout using single load case for example 2 ....................................... 122
Problem description with uniform layout for example 3 ................................. 122
Problem description with multi-size layout for example 3 ............................... 123
Optimal layout using single-size substructures for example 3 ......................... 123
Layout using multi-size substructure for example 3 ........................................ 124
Problem description for example 4 ................................................................. 125

4428
4429

430

44 4-
b.)

D.)
k)

b)
b)

440

(A

(b

Pn
Ini
lay
Sec
()pi
SITU
Star
corr
Seqt
()pfii
Sch
subst
LayOI

1.15an

 

4.27

4.28
4.29

4.30

4.31
4.32

4.33
4.34

4.35
4.36

4.37

4.38

4.39

4.40

Starting arrangement of substructures for optimization using the dividing

approach for example 4 and the corresponding optimal layout ....................... 126
Sequence of steps in the dividing approach for example 4 .............................. 127
Optimal arrangement of substructures using the dividing approach for

example 4 and the corresponding optimal layout ............................................ 128
Starting arrangement of substructures using the merging approach for

example 4 and the corresponding optimal layout ............................................ 129
Sequence of steps in the merging approach for example 4 .............................. 130

Optimal arrangement of substructures using the merging approach and the
corresponding optimal layout for example 4 .................................................. 131
Problem description for example 5 ................................................................. 132

Initial layout using the merging approach and optimal structure using this

layout for example 5 ...................................................................................... 133
Sequence of layouts obtained using the merging approach for example 5 ....... 134
Optimal arrangement of substructures and the corresponding optimal

structure obtained using the merging approach for example 5 ........................ 135
Starting layout of substructures for the dividing approach and the

corresponding optimal structure for example 5 .............................................. 13 5
Sequence of layouts obtained using the dividing approach for example 5 ....... 136

Optimal arrangement of substructures and the corresponding optimal

structure obtained using the dividing approach with prescribed number of
substructures of size 2L for example 5 ......................................................... 137
Layout of substructures and corresponding optimal structure obtained

using the dividing approach with a perimeter constraint for example 5 ........... 138

xii

Cl

lnt

In the
optimiz
prevent
an autt
CompUll
for exte
cOfltplelt
POSsible
minimUm

ln
44am

is invoked

Chapter 1

Introduction

In the design of complex structural systems or structures with complex behavior to
optimize certain properties such as crashworthinees, the complexity of the problem often
prevents a detailed computational analysis of the structure. A typical computer model of
an automotive structure involves around 106 degrees of freedom requiring days of
computer time to run a single analysis. Under these conditions the current practice calls
for extensive simplifications. The strategies for simplification of the optimal design of
complex structures (or structures with complex behavior) can be divided into two
possible approaches: one that uses a full-scale model but limits the evaluations to a

minimum and the other that uses a reduced model.

In the first approach, the optimization is applied on surrogate models (or response
surfaces), constructed using a strategic (statistical) sampling of a lull-scale model, which

is invoked only sparingly to save effort. The firll-scale model is still used as the principal

source 0'
evaluatior
very effer
[40]). As
to build

attractive ;

Th
full-scale r
noted that
inforrnatior
rePlaces. 1
0f accuraq
Same Ph)'si

the results c

It is
one that a
mode'15) eff.
Satings that
PIGSemS
representatic
Preliminary, i
it ignOres fee

NeVetheleSS

source of information about the system and therefore the computational cost per
evaluation is not changed. Instead, the number of evaluations is reduced. This strategy is
very effective whenever the number of design variables involved is small (see Yang
[40]). As the number of variables increases, the number of firnction evaluations required
to build a reasonable response surface also increases, rendering the approach less

attractive and not of much advantage compared to using the firll-scale model throughout.

The main idea in the second strategy is to reduce computations by replacing the
full-scale model by a reduced model, a model that is less expensive to evaluate. It is
noted that the process of constructing the reduced model may result in the loss of critical
information, causing the reduced model to be significantly less accurate than the model it
replaces. In this strategy the computational cost per evaluation is reduced at the expense
of accuracy. In addition, design variables in the reduced order model may not have the
same physical meaning as the variables in the original problem, making interpretation of

the results diflicult.

It is likely that the most effective strategy for the design of complex structures is
one that combines the two approaches (i.e. , response surface methods and reduced
models) effectively. A response surface methodology would clearly benefit from the
savings that result fi'om a carefirlly crafted, reduced model. With this in mind, this work
presents a formal methodology for model reduction that involves an explicit
representation of the model at various scales. The work presented here is still
preliminary, in the sense that it addresses only the elastic behavior of a structure and thus
it ignores features that are crucial to the full understanding of many complex problems.

Nevertheless, this is a necessary first step and it provides an important understanding of

how a

structu

T
quite sr
primaril
reductic
methods
freedom
degrees
freedom
componer
the dynan
cOITtponen
93d] inten
PfOblem is
PTOpenies r
Coarse]- mor

average the .

The Cc
from fine‘xa

dlflEIEnt Seale‘

 

how a formal procedure can be derived to reduce the complexity of models used in

structural design without losing information that can be relevant in the design problem.

The problem of constructing reduced models of structures has been investigated for
quite some time in the context of vibration and control of structural systems. These
primarily comprise of methods such as pseudo static variants of the classical Guyan’s
reduction method or component-mode based techniques. In the Guyan’s reduction type
methods (see Guyan [22], Friswell et a1 [19], Wilson et al [39]), a set of degrees of
fi'eedom of the system are chosen to be master and another set chosen as the slave
degrees of freedom. The reduction process aims to eliminate the slave degrees of
freedom and express the state equation in terms of the master degrees of freedom. In
component mode based techniques (see Hurty [24], Craig and Bampton [8], Seshu [35]),
the dynamics of a structure are described by selected sets of normal modes of individual
components of the structure, plus a set of static vectors that account for the coupling at
each interface where individual components are connected. In the present work, the
problem is phrased in the language of homogenization and the computation of effective
properties of composites: starting fi'om a structure that is modeled in fine detail, we seek
coarser models of the same structure (i.e., “homogenized” or “effective” structures) that

average the detail without losing the (relevant) fine scale information.

The collection of mathematical methods for extracting the coarse-scale behavior
from fine-scale models is termed homogenization. Typically, such problems are solved
using asymptotic expansion techniques or weak-limits; see Bensoussan et a1 [5]. In these
techniques, there is no accounting for structures with features involving distinctly

different scales. In a recent paper by Pecullan et al [30] the subject of scale effects on the

behavior
tensors <
volume r
numerical
appear in
homogeni;
resolution
linear hom
system of e
original 53,-:
these efl‘ect
could deter
Complete Cl:
[30] applied
eqmvalent tr
method of 1
developed m
more rObUSt
situations tha
techqu to .
homogenizat,

that it Com d l

Dorabantu u‘

behavior of two-dimensional composites is discussed by comparing the apparent stiffness
tensors of two-dimensional elastic composites for various sizes of the representative
volume element. Also, there has been substantial effort to develop methods for
numerical homogenization that facilitate the analysis of problems involving systems that
appear in multiple scales. Brewster and Beylkin [7] outlined a procedure for numerical
homogenization of a system of linear ordinary difl‘erential equations using the multi-
resolution analysis (MRA) associated with the construction of a wavelet basis. Their
linear homogenization procedure consists of an algorithm to produce an effective linear
system of equations whose solutions are the coarse-scale projections of the solution of the
original system of equations, called multi-resolution reduction, and one for augmenting
these efi‘ective equations to produce the homogenized solution, called augmentation. One
could determine the projection of the solution at any intermediate scale and obtain a
complete description of the transition from fine to coarse scale representation. Gilbert
[20] applied the same approach to a system of two ordinary differential equations that is
equivalent to a one dimensional second order elliptic problem and compared the classical
method of homogenization, i.e., the asymptotic expansion method, with this recently
developed multi-resolution technique. It was noted that the MRA scheme is physically
more robust than the classical theory, i.e. it could be applied to many more physical
situations than the classical theory. Dorabantu and Engquist [15] applied the same MRA
technique to a discrete elliptic second-order differential equation. They observed that this
homogenization procedure produced an operator that preserved its divergence form and
that it could be well approximated by a band diagonal matrix. The work by Gilbert and

Dorabantu use the Haar basis (piecewise constants) for the discretization of the

djfi‘eren

difiereni

In
of elastic
using fix:
Glowinski
work, the
elasticity n
could be it
large stitfne
generalized
material in 2

this problem
10 the size 01
is Specifican}

me Principa] ‘
To deve}
83.516015

finescar.

 

differential equations. Beylkin and Coult [6] applied this MRA method to elliptic partial

differential equations and studied the spectral characteristics of the reduced operators.

In linear structural analysis, the associated differential equations are the equations
of elasticity. Techniques for solving the elasticity problem defined on arbitrary domains
using fixed-scale wavelet-Galerkin methods have been investigated for some time (see
Glowinski et a1 [21], Wells et al [37], Diaz [14], DeRose and Diaz [11]). In the present
work, the MRA based numerical homogenization scheme is applied to the equations of
elasticity modeled using a wavelet-Galerkin technique. In this case, the model reduction
could be thought of as a method to produce equivalent smaller stiffness matrices from
large stiffness matrices. An application of this method is then presented in the context of
generalized topology or layout optimization of structures, i.e., the optimal distribution of
material in a given design space subject to prescribed loads. The existing methods for
this problem do not account for the presence of finite scales that may even be comparable
to the size of the design domain. Here, a method that uses a model reduction scheme that
is Specifically tailored to the problem of layout optimization of structures such that finite

scale heterogeneities can be accounted is presented.

The principal goals of this dissertation are:
1. To develop a consistent scheme to compute equivalent reduced models of structural

systems in linear elasticity at various coarse-scales that retain relevant features of

fine-scale models.

accc

size

The rema.
introductic
Chapter 3
framework
Chapter 4
problem 0

heteIOgene

possible dir

2. To develop a method for the layout optimization of structures in elasticity that
accounts for multi-scale heterogeneities of sizes that may even be comparable to the

size of the structure.

The remainder of this dissertation is organized as follows. Chapter 2 gives a brief
introduction to wavelets and the concepts of multi-resolution analysis of functions.
Chapter 3 presents a model reduction scheme using the multi-resolution analysis
framework. Numerical emples that illustrate the proposed scheme are provided.
Chapter 4 discusses the application of the proposed model reduction technique to the
problem of layout optimization of two-dimensional structures in elasticity, in which
heterogeneities of finite scale can be accounted. Finally, some concluding remarks and

possible directions for firture work are presented in Chapter 5.

 

"H
freq
Stu}
“av.
hiera
analc

l Daub

Chapter 2

Wavelets and Multi-Scale

Representation of Functions

2.1 Introduction

“Wavelets are mathematical fimctions that are used to cut up data into diflerent
frequency components and then study each component with a resolution matched to its
scale Ӥ. The concept of multi-resolution analysis is firndamental to the theory of
wavelets. The main idea is the separation of the information to be analyzed
hierarchically into principal and residual parts. In signal processing applications this is
analogous to decomposing a signal into its low fi'equency and high fiequency

components with the knowledge of when they occur. This has definite advantages over

 

’Daubechies [10]

the standai
information
making “*8
Fourier trar
been appliet
solving ordi
emphasis 0
chapter. F t
refer to the

references t1

2.2 Mull

DCfiHIIIOn;

filnctions ~ i.

(i) {0;
(ii) fl}
(iii)g(
(1") st.»

(V) The

the standard Fourier analysis, which identifies fi'equency information but no time
information. A variety of efiicient algorithms using wavelets have been developed
making wavelet transforms on par with computationally efficient methods such as fast
Fourier transforms. Wavelets have become a very popular tool in engineering and have
been applied to a wide range of problems in signal processing, image processing and in
solving ordinary and partial differential equations. A brief introduction to wavelets with
emphasis on applications in computational engineering analysis is presented in this
chapter. For detailed information about the construction and applications of wavelets
refer to the books by Daubechies [10], Frazier [18], Resnikofi' and Wells [31] and the

references therein.

2.2 Multi-Resolution Analysis

Definition: A multi-resolution analysis (MRA) of L2 (IR) - the space of square integrable

firnctions — is a nested sequence of subspaces Vj such that

(i) {0} c ---C V_1cV0 ch C---CL2(IR)

(ii) njVj={o} and LIV—1:3 (1R)

(iii) g(x)€ VI 4:) g(2x)€ Vj+1

(iv) g(x)€V0 4:) g(x—k)€ V0, k E Z

(v) There exists a scaling function <p(x)EV0 such that

{90(x-k), k E Z} is an orthonormal basis for V0

Thus, using the above definition it be can shown that if

V0 2 span{<p(x—k), k E Z}

then

Since V0 C

Thus.

where a). E
the dilation
VJ = span { .

A firnction j

PJf apprOac

Define a "e“!

 

then
V1 = span{<p(2x—k), k E Z}
Since V0 C V1 , any function in V0 can be expressed in terms of the basis function of V1.

Thus,

<p(x)= i ak<p(2x—k) (2.1)

k=—oo

where ak 612(1R) (the space of square summable sequences). This equation is known as
the dilation equation or the scaling relation. Denote: cp J, k (x) = 2‘“ 24;)(2Jx— k) and
V J = span {CpJ’k (x), k E Z}. The dilation parameterJ is known as the scale.

A firnction f E L2 (R) may be approximated by its projection onto the space V0 as
00
P0f= Z co’kcp(x—k) (2.2)
k=—oo
and more generally by its projection onto the space V J as
00
PJf = Z cJ,k (Pu: (x) (2-3)
k=—oo

PJ f approaches f as J —+ co , i.e., higher (finer) the scale of the representation, better
the approximation. Consider now the difference between the subspaces V J__1 and V J.

Define a new subspace WJ_1 such that it is the orthogonal complement of VJ_1 in VJ ,
i.e.,

VJ = VJ—r 69 “0-1

(2.4)
VJ—r i WJ—I

It followst
firnction L'
subspace H
for H}. ll

orthonormal

Then from e

This equatior
a level J -1
since the spar

terms ofthe s

2-3 Exarr

The SimpIeSI

filnction and

 

COUStant funCI

 

It follows that the spaces W j are orthogonal and that G) W] = L2 (IR). Define a wavelet
jEZ

firnction tb(x)EW0 such that {¢(x—k), kEZ} forms an orthonormal basis for the

subspace W0. Then, {¢J,k (x): 21/2¢(2Jx—k), k EZ} forms an orthonormal basis
for WJ. In addition, it follows that {am (x)=2f’2¢(21x—k); j,k€Z} forms an

orthonormal basis for L2 (R). Denote the projection of a firnction f on W J as Q J f .
Then, from equation (2.4) we have

PJf = PJ—1f+QJ—rf (2.5)
This equation means that Q J_l f represents the detail that needs to be added to get fiom

a level J —1 representation of the function to a level J representation. Furthermore,
since the space W0 is contained in the space V1 , the wavelet firnction can be expressed in

terms of the scaling function at the next higher scale,

on): {If Wax—k) (2.6)

k=-oo

2.3 Example Using the Haar Basis

The simplest possible orthogonal wavelet system is generated from the Haar scaling
firnction and wavelet shown in figure 2.1. The Haar scaling function is a piece-wise

constant firnction defined as

1 OSx<1

. (2.7)
0 otherwrse

90(X)=

10

This firnctic

Similarly th

ltcanbese

 

This fimction satisfies the scaling relation (equation (2.1)) with coeficients ao = a1 =1,
90(x) = <p(2x)+ <p(2x-— 1) (2.8)
Similarly the Haar wavelet is defined as

1 05x<05
1/2(x)= —1 0.5$x<1 (2.9)
0 otherwise

It can be seen that the Haar wavelet satisfies equation (2.6) with b0 =1 and bl = -1,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

¢(x)=‘P(2x)-‘P(2x“1) (2.10)
1 s ................ .. -, I 5 .....
1.25» { 1:
1t
0 5
o as L *
0.5 o(
0.15
. -o.s
o’ t
E - l
-o.zs» 1:
’ l
-1 -o s o o s 1 1 5 z -1AA‘QoisidiioA‘iAo‘siwiAiiigi1A5 2

Figure 2.1: Haar scaling firnction and wavelet
Figure 2.2(a) shows an arbitrary firnction in its original form and using equation (2.3) this
firnction is approximated by its projection in the space V6 (26 coeficients) as shown in
figure 2.2(b). Figures 2.2(c) and 2.2(d) show the projection of this firnction in the spaces
V5 and W5 respectively, i.e., the next coarser representation and the detail that has been

removed from the representation at scale—6 to create the coarser representation.

11

€an
0 6'-

04'

0/
412-
41.4-

«0.6-

 

418
0

Figu

0.8%.
0.6!-
0 Cr

02*

412»
414'-
0 6-

08\\_

Figllrl

2.4 c;G
Wa\Iel
1“ applica

comPrOrm

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0:2 (E4 ofe 0:8 1 ' 0 0‘2 01.4 ofe ofa 1
Figure 2.2(a): Original function Figure 2.2(b): Approximation at
scale 6 (64 data points)
1 . 1 r 1 0.5
0.8 ‘ 0.4+
0.6 ‘ 0.3].
0.4 ‘ 0.2,, J
0.2% l 0.1 i
0 t o
43.2) -0.1 i
-O.4 l- r .02 t.
-0 8 -O.3[
418° 0‘1 04-" 01-3 0:8 1 ”“0 03 01.4 oie 0:8 1
Figure 2.2(c): Approximation at F' 2 2 d . Detail t seal 5
scale 5 (32 data points) rgure ' ( )' a e-

2.4 Generalized Orthogonal Scaling Functions and

Wavelets

In applications involving functional analysis there is a need for basis firnctions with other
properties besides orthogonality such as continuity and ditferentiability without

compromising the compactness of the scaling firnctions and wavelets involved. Some

12

important pr

scaling fiJnCt
constant we

on the scaling.

where

In order to u

function is us

Which leads 1.

ECluation (2,1

 

important properties of such scaling firnctions are discussed in this section. A general

scaling firnction <p(x) (of dilation factor 2) is a solution to a dilation equation (2.1). The
constant coefficients ak (filter coefficients) are derived by imposing certain conditions

on the scaling firnction such as orthogonality with respect to integer translations,

fgo(x)<p(x+l)dr=60’1 vrez (2.11)
where
o —1 1:0 (212)
0’1— 0 otherwise '

In order to uniquely define scaling firnctions of a given shape, the area under the scaling
function is usually normalized to unity, i.e.,

fcp(x)dr=l (2.13)

which leads to the following condition on the filter coefficients
()0
Z ak = 2 (2.14)
k=—oo
Equation (2.11) results in the following condition on the filter coefficients
OO
2 akak+21 = 250,, v1 6 Z (2.15)
k=—oo

Using the orthogonality of the wavelet function, 4/2( x) is defined as

00

¢(x)= Z (-—1)" aN_1_k <p(2x—k) (2.16)

k=—oo

l3

where N i5
the interval

.. 1.-
and ((—1) i

The above

coefficient Sj
wavelets, the

be able to e)

the wavelet a

  

”Oments prr
difi‘erentiabili

thEn

 

Thug, equatic__

Coeflicient SV
t‘anishjng

approx"'Tlatior

 

where N is an even integer and the scaling firnction and the wavelet are non-zero over

the interval [0, N — 1] called the support of the functions. The sets of coefficients {ak}
and {(— 1)k aN_1_ k} are said to form a pair of quadrature mirror filters.
The above equations still do not yield a unique set of filter coefficients. In an N

coefficient system, they yield a total of 1:- +1 equations. For the Daubechies family of
wavelets, the other %—1 equations are determined by requiring the scaling firnction to

be able to exactly represent polynomials of order up to %. Using the orthogonality of
the wavelet and the scaling fimction this leads to the following condition on the wavelet,

f£¢(x)=o 1=o,1,...,%-1 (2.17)

i.e., the first % moments of the wavelet must be zero. This is known as the vanishing

moments property of the wavelet and is closely related to the smoothness and

differentiability of (,0 (x) and (,1) (x). The associated constraint on the filter coefficients is
then

00
Z (—1)"akk’=o 1:0,1,...,-12X—1 (2.18)

k=—oo
Thus, equations (2.14), (2.15) and (2.18) uniquely define the filter coefficients for an N
coefficient system. In addition, it is possible to construct wavelet systems by enforcing
vanishing moments properties on the scaling firnctions. This results in better

approximation of the expansion coefizicients by samples (rather than by orthogonal

I4

projections)

of wavelets i

2.5 The

The single 5
means of Ira:
of resolution
for transfom

resolution re

Consider a fi

”d QJf de

where (,3)

 

From eqUalio

Substjmline ( .

 

 

projections) and it causes the scaling firnction to be more symmetric. The Coiflet—family

of wavelets is an example of such systems.

2.5 The Discrete Wavelet Transform

The single stage discrete wavelet transform algorithm of Mallat [28] provides a simple
means of transforming functions from one level of resolution, J , to the next coarser level
of resolution, J —1 , called multi-resolution decomposition or the analysis phase and one
for transforming functions from the coarser level back to the finer level, called multi-
resolution reconstruction or the synthesis phase.

Consider a function f and let R, f denote the projection of f onto the subspace VJ

and Q J f denote the projection of f onto the subspace WJ. Thus,

PJf = Z 0.1,): <PJ,k (x), Cu: = (f .9014) (2-19)
k=—oo

QJf = Z dJ,ki/1J,k(x)a dJ,k =(f,'¢/1J,k) (2.20)
k=—oo

where (0,0) denotes the Euclidean inner product. Since WJ_1 is the orthogonal
complement of VJ_1 in VJ ,

PJ—lf =PJf -QJ—1f (221)

From equation (2.19), it can be seen that

CJ—r,k =(PJ—1f,¢J—Lk) (2-22)

Substituting (2.21) in (2.22),

15

 

Also. it ca

Similarly.

in the equz

leads to th

These eqL

aJSOIithm i
Consider a
b} a "6th

equaIlOn (:

W~‘C‘Ior of Cc

l (X)
cJ—l,k=‘\7§' Z: cJ,jaj—2k (2.23)

j=—oo
Also, it can be shown that
1 0° '
dJ—l,k =— E cJ,j(_1)j aN—l—j+2k (2.24)
J2 -__.
j— 00
Similarly, substituting the relation
PJf = PJ—lf + QJ—rf (2.25)
in the equation
Cm = (PJf,sDJ,k) (2.26)

leads to the relation
1 0° 1 0° 1:
CJ,k=— Z cJ—l,jak—2j+— Z: dJ—l,j(_1) aN—l—k+2j 0-”)
J2 .__ 45 -__
j— 00 j— 00
These equations (2.23), (2.24) and (2.27) form the basis of the Mallat transform
algorithm implemented using the signal-processing concept of filter-banks as follows.
n—l

= Z 0J,k 901,1:
k=0

Consider a periodic, finite dimensional firnction f E V J, represented

 

 

by a vector of coefficients c J , where n = 2‘] is the period of the function. Then,

equation (2.23) in the discrete periodic form represents a circular convolution of the
vector of coefficients with a discrete filter of the form

.. l T
h=——- 0 0 0 a _ a 2.28
Jilao N r 2 01] ( )

l6

 

followed

vector)-

discrete 5

followed
discrete v
single-stag

of such fili

followed by a dyadic down sampling (keeping only every other entry of the resulting
vector). Similarly, equation (2.24) is represented by the circular convolution of another

discrete filter

1 T
—_-—a _ o o o — a ._ —a _ 2.29
g ,—2[N1 00 N3 N2] ( )

followed by dyadic down sampling. Thus, the analysis phase (decomposition) of the
discrete wavelet transform can be illustrated graphically as shown in figure 2.3 for a
single-stage transform. Further decompositions are implemented by applying a cascade

of such filter-banks recursively to the coarse-scale approximation coefficients, c j- .

 

 

 

 

 

‘r——’ *3 ———q 12 +—> arJ_l

 

 

 

CJ——-’

 

 

 

|——D 42h ——u 12 ._, cJ_1

 

 

 

 

 

 

Figure 2.3: Single-stage discrete wavelet decomposition

In the figure, following the stande notation, 4: represents the convolution operator

illustrated in figure 2.4 and i represents the down sampling operator defined as shown in

figure 2.5; all the indices and arguments are evaluated modulo n.

 

x(i) ——> *h .__., y(i)=Zh(k)x(i—k)

 

 

 

Figure 2.4: The convolution operation

17

 

Similarly,
Wavelet t

as illustra

and

Here T d.

denOte

x(n) a} ytn)=x(2n)

Figure 2.5: The dyadic down sampling operation

 

Similarly, the discrete periodic form of equation (2.27) leads to the inverse discrete
wavelet transform and is implemented by a sequence of convolutions and up samplings

as illustrated in figure 2.6. The filters in this case are

l T
hzfilao 01 02 div-2 aN—l 0 01 (230)
and
1 T
g=:/—§[aN—I —aN_2 aN_3 01 —a0 0 "° 0] (2.31)

 

 

dJ—1——h 12 ——> *g

 

 

 

 

 

 

CJ

 

 

CJ—1——N T2 *—> *h

 

 

 

 

 

 

Figure 2.6: Single-stage discrete wavelet reconstruction

Here T denotes the up sampling operator and can be defined as shown in figure 2.7 and

619 denotes the addition operator.

l8

 

1n the p:
hence all
strictly

convolut
particula:
input vet
inverse t
However
influence
“1115me
Using Var

31C.

2.6 Co

VNRIVnal

The explic
onh'the E
accurate

“(fights c

 

x(n) ——> 12 L——->y(n): y(2n)=x(n)

 

 

 

Figure 2.7: The dyadic up sampling operation

In the preceding discussion it is assumed that the associated functions are periodic and
hence all the indices and arguments are evaluated modulo n. It is noted that this is not
strictly necessary and one could define the same transformations using linear
convolutions rather than cyclic convolutions. In the periodic case, the formulation is
particularly clean because the total number of terms in the transformed vectors and the
input vector are always the same, the transformation matrix is square with a simple
inverse that has an interesting structure and can be efficiently calculated using an FFT.
However, there is the additional problem due to aliasing in the periodic case, this is the
influence of the terms in one end of the function affecting the other end due to cyclic
transformations. These problems are sometimes overcome by extrapolating the functions
using various techniques such as zero padding, symmetric padding, reflective padding,

etc.

2.6 Constructing the Basic Scaling Function and

Wavelet

The explicit use of the scaling firnction or wavelet is rare in most applications (one uses
only the scaling and wavelet filter coefficients); however, they might be required for
accurate function evaluations and visualization. In general, scaling functions and

wavelets do not have a closed-form solution. Instead, they have to be computed

19

 

recursiw

written t

Evaluati:

;(ij)=<

In matrix

01'

In Order tC

recursively from the dilation equation (equation 2.1). The dilation equation can be

written explicitly as

<p(x) = aocp(2x)+ a190(2x— l) + - - - + aN_1cp(2x — N +1) (2.32)

Evaluating the above expression for all integer values x = j E Z , it can be shown that

<p(j)=0 for j<0 and j>N—1. Thus, the only remaining equations are

 

<p(1)=a090(2)+an0(1)+0290(0)
cp(2)faosa(4)+ar<p(3)+ach(2)+03w(1)+a4s0(0)

 

 

 

 

 

(2.33)
<p(N—2)= aN_3cp(N——1)+aN_2<p(N—2)+aN_1<p(N—3)
cp(N—I)= aN_1<p(N—1)
In matrix fornr, this can be written as
rno o o o o o " 90(0) ' ' 90(0) ‘
a2 a1 a0 0 0 0 (0(1) <p(1)
a4 a3 (12 0 0 0 30(2) <p(2)
2 I Z 2 I I I = (2.34)
0 0 0 “iv-3 “iv-4 “Iv-5 $(N—3) iPiN-3)
0 0 0 div-1 div-2 div-3 MAI-2) 90(N-2)
_o o o o o aN_1J_cp(N—1)J “(0(N—1)‘
Ol’
M024) (2.35)

The vector of integer values of the scaling firnction (ID) is then given by the eigenvector

of the matrix M corresponding to the eigenvalue 1, i.e. the solution to the system

(M—I)<p= o (2.36)

In order to uniquely determine a solution to the above system a normalizing condition

arising from equation (2.13) yields the condition

20

Thus. the V
(2.36) no

half-integer

This proces

points L
2!!

wavelets.

Z <p(j)=1 VjEZ (2.37)
k=—oo

Thus, the values of the scaling function at integers are given by the solution to equation

(2.36) normalized by equation (2.37). Once these are computed, the values of <p(x) at

half-integer points can be determined using the dilation equation as
x 00
“OH = Z ak(,0(x—k) (2.3s)
k=—oo
This process is repeated as many times as necessary to find the values <p(x) at all dyadic

J

points (5;; j,n E Z]. Figure 2.8 shows some commonly used scaling functions and

wavelets.

21

 

 

44441141

11144414441

 

 

 

 

 

 

 

 

' g fimction and wavelet

Daubechies-4

 

 

 

 

 

 

vvvv‘vvvw'

riervav

Ii 'f}’h>

       

DD DD Di iiiib

AAAAAAAAAAAAA

 

 

1.85

0.st

 

'1

g firnction and wavelet

'1

Daubechies-6

 

' f v v v v r

v v ‘ v w v v v r v '

 

141444411111 4444111414144414411

1
A A A A A A A A A A

A A A .-

 

 

 

 

 

.
A

A

. A
. A
. A
. A
v A
. A
. A
. A

r. 5 0 8 1 5 2
a g - u g
0 0 1
. _

. t .................. L
T A
a A
v A
j A
n

A

A

A

 

 

'J.

Coiflet-Z scaling function and wavelet

 

 

 

4411141441441444414144

A A

A

A

A A

A A _ A

A

itblfihitiiifii’hiii’

 

 

 

 

A A

l A A

A

 

 

Coiflet-8 scaling function and wavelet
Figure 2.8: Some commonly used orthogonal scaling functions and wavelets

'9

22

 

2.7. F 1
Projec

In certair
poly-nomiz
multi—resc
approximz
illustrated

analysis ar

Let 1"”
respeCIIVQI'
v . FOr

Element fur

2.7. Function Approximation Using Orthogonal

Projections

In certain applications, functions expressed in some commonly used basis (such as
polynomials) need to be approximated in a suitable scaling function basis so that the
multi-resolution analysis described earlier could be applied to them. Here, the
approximation of functions in a scaling firnction basis using orthogonal projections is
illustrated using a basis of bi-linear (hat) functions used commonly in finite element

analysis and the Daubechies D6 scaling functions.

 

Figure 2.9: Orthogonal projections from wavelet to finite-element spaces
Let V w and V h be finite-dimensional spaces of periodic functions in 1.2 (1R) , spanned
respectively by D6 scaling functions cpw and bi-linear (finite-element) shape functions
h w w Ph w h - - -
(,0 . For any f EV , let f EV be the prOJectron onto the space of finite
element firnctions such that the L2 -norm of the error e, (e = f w - Ph (f w )) between the

original firnction and the projected function is minimized, i.e.,

n—1 n—1
W” = Z fkw sot". P” (1"): Z fi'soi' (2.39)
k=0 k=0

23

 

The Static

In the vet

Here C
wavelet a

integrals a
the COmpL
Problem; f
‘3‘ a! [27]
transforms

as

Consider n'

Find f)? 6 IR that

 

 

 

n—1 n—1 h h 2 (2.40)
min f 23 fkwwi" (y)—§: fk 90k (y) dy
k=0 k=0
The stationary point of equation (2.40) is the solution to
n—l h n-l h h h
f 23/13"in 90j03’=f ka c.01.: ‘pjdy (2-41)
k=0 1:20
In the vector form this can be written as
Ct" = Nih (2.42)
4:) rh _—. N“Cr“’ (2.43)

Here C and N are block-circulant matrices comprising of the inner products of the
. w h h h ‘
wavelet and hat functions, f 80k cpjdy and f 90k cpl-dy , respectrvely. Values of these

integrals are commonly referred to as connection coeflicients. Using the scaling relation,
the computation of these quantities is usually reduced to the solution of an eigenvector
problem; for details regarding the computation of these connection coefficients, see Latto
et a1 [27], Dahmen and Michelli [9], Kunoth [26]. Thus a projection matrix that
transforms vectors in the wavelet space to that in the finite element space can be written
as

P” = N‘ 1C (2.44)

Consider now, the projection, Pw ( f h) of a given finite-element firnction f h E V b onto

the wavelet space. As before, the statement of the projection problem can be written as

n—1 n—1
Vf” = 2 that e V”. P” (f”)= Z) gin" (2.45)
k=0 k=0

24

 

The statio

(Note: (4‘)

(Compare
that transfr

written as

 

Findgi" 6 R that

 

n—l h h n—l 2 (2.46)
min] 2 ft 90k (”-2 giver? (y) dy
The stationary point of equation (2.46) is the solution to
n—l h h n—l
f 21m spray= f Eire}! wydy=ff (2.47)
k=0 k=0

 

 

(Note: (up, k e [0, n — 1]} forms an orthonormal basis for V‘” ). In vector form this is,

CTrh = r‘” (2.48)

(Compare the left hand sides of equations (2.41) and (2.47)). Thus a projection matrix

that transforms vectors in the finite-element space to that in the wavelet space can be
written as

P” = CT (2.49)

where, the matrix C is the same as the one defined earlier.

2.8 Wavelets In Multiple Dimensions

In the preceding discussions it is assumed that the associated firnctions are one-
dirnensional. In higher dimensions, the wavelet and scaling functions are defined as
tensor products of the respective functions in 1-D. The space of two-dimensional

functions at a scale J is given as
VJ =VJ®VJ (2.50)
where V J is the corresponding space of 1-D functions. For example, a two-dimensional

scaling firnction is defined as follows:

25

 

All the
functions

fixed (fine

It follows

J—l can

Thus. a mu

notation as,

“here the c'

and the deta

 

 

‘PJ,k1 (X, y) = <PJ,k (10%,! (y) (251)
All the concepts discussed can then be extended to this tensor product of scaling

functions. Thus a two-dimensional periodic firnction, u (x, y) , can then be expressed at a

fixed (fine) scale J as

N-l N—l
uJ(x,y)= PJ“(xy)= Z ”k190J,kl(xJ’)=Z 141%ka )‘PJ,I()’) (2-52)
k,l= o k,l=0

It follows from using (2.4) in (2.50) that the coarse-space and the detail space at scale
J —1 can be expressed as

VJ—r = VJ—1®VJ—1

(2.53)
WJ—r = (WJ—r ®WJ—1)€9(VJ—1 ®WJ-1)€B(WJ—1 ®VJ—1)

Thus, a multi-resolution representation of this function can be expressed using a compact

notation as,
—1 J—l 2'"-l3
“J(x}’ =2: fir/90M: xy)+ +2: 2: Zum,ki¢m,(ki (’0’) (2.54)
k,l=0 m=jk,l=0r=l

where the coarse-scale basis function,
‘10],kl (x,y) = Soj,k (1090),! (y) (255)

and the detail basis functions are denoted as

T/Jrln,k1(x,y)= "r/Jm,k (x)‘pm,l (y)
who (x.y)=<pm,i (44..., (y) (2.56)
thine (350’) = ¢m,k 0‘)va (Y)

26

 

It should
reduces ti

decompo:

 

Origin;

Figu

 

 

It should be noted that each reduction of a discrete, periodic two-dimensional firnction
reduces the size of the coarse-scale vector of coefficients by a factor of 4. A two-stage

decomposition of a two-dimensional firnction can be illustrated as shown in Figure 2.10.

 

 

 

 

 

<PJ—2 ® iPJ—z
<PJ—2 8’ SOJ—Z
i”:- ‘ 110-2 69 <PJ—2
¢J— ® ill _
<PJ ®<PJ 2 J 2
Original function at scale J Two—scale representation of the

firnction
Figure 2.10: A multi-scale representation of two-dimensional functions

27

 

Mor

This chat
the analy:
resOiutiOn
and behax
method Ca
blocks to ‘
«Sign for
efficiency

are strong].
bodies thjS

analySis of Si

Chapter 3

Model Reduction in Elastostatics

This chapter deals with methods to consistently reduce the size of some models used in
the analysis of large structural systems. The reduction is accomplished using a multi-
resolution analysis applied to differential operators to average the fine-scale properties
and behavior while retaining the coarse-scale information. In the discrete form this
method can be used to construct small stiffness matrices that can be used as building
blocks to construct larger systems. An application of this scheme is proposed to be in
design for crash-worthiness, where intensive computations demand computational
efficiency. It is noted that the models of structural systems for crashworthiness analysis
are strongly nonlinear and involve other phenomena such as contact between rigid
bodies; this work is only a starting step towards such a comprehensive scheme for the

analysis of such models by considering the analysis of structures in linear elastostatics.

28

This C
introdu'
problerr
multi-re:
section L
the displ
proposed
in section
suitable fc
parameters
Finally, nu

Presented ii

3.1. Sca

LCI S repre

be mOdeled

Let SJ be

 

Here the in d)

discretlZed 8

Upon dlSCre‘

 

This chapter is arranged as follows: The notion of scales in structural systems is
introduced in section 3.1. This is followed by a discussion of the discretized elasticity
problem at a fixed (fine) scale in section 3.2. A model reduction technique based on a
multi-resolution analysis of the material distribution is presented in section 3 .3. In
section 3.4, a model reduction strategy based on a periodic multi-resolution analysis of
the displacements is presented. This is followed by comparisons between the two
proposed reduction schemes as well as that with the classical homogenization technique
in section 3.5. Section 3.6 deals with computing a reduced stifliress matrix that is
suitable for assembly with other stiffness matrices. A technique to compute fine-scale
parameters using the reduced solution, called augmentation, is presented in section 3.7.
Finally, numerical examples that illustrate and compare the methods discussed are

presented in section 3.8.

3.1. Scales And Material Properties In Elasticity

Let S represent a structure of interest. The relevant behavior of S in linear elasticity can

be modeled by a linear system of the form

L(u) = f + boundary conditions (3.1)

Let S J be a model of the structure that incorporates all of its details up to a scale J .
Here the index J is used to denote the level of refinement at which the structure has been
discretized and it is chosen such that all the relevant features of the structure are resolved.
Upon discretization, the elasticity equation (3.1) becomes

LJuJ = g (3.2)

29

 

where l
displace:
open for

of the st
More exr

to the eqt

accounts 2
in finite el
of certain

is usually c

decide whj

 

where L J is a stimress matrix, f J is the force vector and “J is a vector of unknown
displacements. The specific choice of the basis functions used in the discretization is left
open for the time being. We are interested in whether we can reduce the size (dimension)

of the stiifiress matrix in (3.2) and still retain all the relevant information present in "J.
More explicitly, the question is: Is there an operator H J _ k such that the solution '6 J_ k
to the equation

31—): fiJ—k = fJ—k (3-3)

accounts for all the relevant information present in u‘] ? This question is a familiar one
in finite element analysis, where this is usually accomplished using a static condensation
of certain unwanted degrees of freedom, e. g., see “filson [38]. However, this procedure
is usually done as an ad hoc scheme, i.e., there is no general set of rules that help one to
decide which degrees of freedom to condense out and which to retain. Also, in some
cases, there is no physical meaning or basis associated with the condensed degrees of
freedom. Here we consider an approach that is more closely related to the methods of
periodic homogenization used in the computation of effective properties of rapidly
varying materials. We associate the different scales of detail in the structure S with the
scales present in the geometric layout of the material in the structure. This is facilitated

by assuming that the material distribution in the structure can be expressed as
Eijkl (x)=p(s)E3k1 (3.4)
where Egkl is a known, reference material tensor and p: IR" —+ (O, l] (n is the spatial

dimension of the model) characterizes the spatial distribution of the material. In this

context, p is commonly known as the efi‘ective density function. This material model is

30

 

frequently
Thus, we 2
density fiJn
We
an assembl

3.1.

 

frequently used in solving structural topology optimization problems, e.g.. see [33].

Thus, we associate the scales present in the model to the scales present in the relative
density function.

We are interested in structures whose material distributions can be expressed as

an assembly of several sub-domains (called here substructures) as illustrated in figure

3.1.

Y2

 
 

X
2 (X1129

x2
i—> x1

x1
Figure 3.1: A structure built from substructures of different scales.

Here 9: U06, where no is a substructure. Any point x in QC can be expressed as
x = Tic + y , where Ye is the global coordinate of a reference point in the substructure and

y is a coordinate local to the substructure. Then we can express the material distribution

as

p(x)=p(fc,y) (3.5)
Also, it may be possible (and desirable) to express the stiffness matrix associated with the
structure as an assembly of matrices corresponding to each substructure. Thus the
problem of finding a reduced model of the structure becomes a problem of finding
reduced stifliress matrices corresponding to each individual substructure. The concept of
a structure being divided into substructures is similar to the notion of super-elements in

stande finite-element analysis.

31

The
homogeniz
a mixturec
with a fret
infinitesima
(weak limit
scales in w
the averag
displacemer

described ir

3.2. ml

The Plane-s

Such that

Where 5 ( ll

 

E is the e!

elastic tense

 

 

The above formulation of the problem has obvious similarities to that in periodic
homogenization, where the goal is to find the effective material tensor that corresponds to
a mixture of materials. The mixture is characterized by a cell that is repeated periodically
with a frequency 1/ e, (e —i O) , i.e., material variation is assumed to take place at
infinitesimal scales. Such problems are solved using asymptotic expansion techniques
(weak limits) and the result is a homogenized material tensor. In the present problem the
scales in which the material is distributed in Q are of finite dimensions and the result of
the averaging process is an operator (a stifl’ness matrix) that relates loads and
displacements in the reduced scale. The construction and reduction of the operators are

described in the following sections.

3.2. The F ine-Scale Elasticity Problem

The plane-stress elasticity problem on a prescribed domain 9 = Uflc seeks u E V (9)

such that

fEe(u)e(v)dQ= ftvdI‘ Vv€V°(Q) (3.6)
Q 1‘t

where 5(u) is the strain tensor associated with the displacement u; t is an applied
traction on the boundary 1"; V is a space of kinematically admissible displacements and

E is the elastic tensor defined on 9. Using the material model described earlier the

elastic tensor within each substructure is expressed in the form

E(y)=chO (3-7)

32

 

where EU

a substruct

problem (1

where the
we assume
with jump
called the .s

discretizati
QC lS I’CSOl
is the V'aIUe

that resolv
from subs:
“ye
discretiZed l
discretizEltio

e'g‘w 118ng f

where E0 is a reference material tensor, pc characterizes the material distribution within
a substructure and y E QC is a coordinate system local to the substructure. The elasticity

problem defined on Q is now:

ZfE(y)e(u)e(v)dy=ft.vdr Vv€V0(Q) (3.3)
c Q Pt

where the sum is interpreted in the sense of assembly. In order to facilitate computations

we assume that pc is resolved with sufficient accuracy by a piecewise constant firnction

with jump discontinuities at Cartesian grid lines spaced Sc = 2"J LC units apart (Sc is
called the scale of the discretization), for some positive integer J called the level of the

discretization. LC is the length of a side of the substructure. In practice, the geometry in
QC is resolved by a digital image composed of 2J x 2‘] pixels of size Sc xSc , where pic],
is the value of pc at the center of the pixel (LI). Sc and J are a measure of the scale

that resolves pc. We denote this by writing pc (y) E p3 (y). Both Sc and J may vary

from substructure to substructure.

We refer to (3.8) as the fine-scale problem when all the substructures are
discretized at their finest scale. In the typical problems of interest such a fine-scale
discretization results in a system with too many degrees of freedom for emcient analysis,
e.g., using finite elements. The large size of the fine-scale problem requires that the fine-
scale stifliress matrix associated with each substructure be replaced by an equivalent
stifiress matrix of a smaller dimension obtained through some consistent process of

reduction. Here we propose two such reduction strategies: one based on a multi-

33

resolutic
multi-ref
inexpenS
particula
displacer
resolutior
rigorous

relatively

3.3 t!
Distril:
All appro.‘
USlng a “1

15 resolved

Where the 1
Pixel (k, I)
and detail
implemEme'

Splits timer"
1

resolution analysis of the material distribution functions pc and the other based on a

multi-resolution analysis of the displacements u. The first approach is computationally
inexpensive but crude. It is used where there is no necessity to go from the solution at a
particular discretization level to another, i.e., there is no consistent procedure linking the
displacements at different levels of discretization. The second strategy, based on a multi-
resolution analysis of displacements, is a more consistent procedure and provides a
rigorous link between the displacements at various levels of discretization, but it is

relatively a computationally expensive approach.

3.3 A Multi-Resolution Analysis of the Material

Distribution

An approximation scheme based on a representation of material distribution function pc

using a wavelet expansion is presented. If the material distribution in a substructure QC

is resolved by a 2‘] x 2‘, digital image, we can write

N—l

pc(y)=p5 (y)= E kaz 9011:! (y) (3.9)
k,l=0

where the functions 90.1,“ are 2D Haar scaling fimctions, i.e., piecewise constant over the
pixel (k,I) and N = 2J. A wavelet expansion of this function splits pc into coarse-scale
and detail functions, ,5( y) and i)”( y) , respectively. The transformation is easily
implemented using a 2-D discrete wavelet transform. More specifically, the transform

splits functions pc 6 VJ , a space of dimension 22] = N 2 , into functions 5 6 VJ_1 and

34

 

fiefij_p

dimension
coarse Spa

decomposr

The coarse

tensor prod
Of pc OVe,

filnction [7

lenSOr

 

 

 

 

[)6 WJ_1, where WJ__1 is the orthogonal complement of VJ_1 in VJ. Space V J_1 is of

dimension 220—1) 2 N 2 /4 , i.e., each application of the wavelet transform reduces the
coarse space by a factor of four. After several applications of the wavelet transform the

decomposition is of the form

N —l
y)= Z: pawn (y)
k,1=0
(3.10)
n—l J —l 2'" —l
= Z EWJ-My yl+ )m2 2 Z (’W'ltmfluy
k,l=0 -1 k, 1-0 r-l
The coarse-scale function is of the form
:20?“ also, “(y (3.11)

and n=2f . Here the 2D scaling functions $131610) and the wavelets rpm (y) are
tensor products of their corresponding 1D firnctions. Coefiicients 'p'k, are average values
of pc over pixels of size A = 2‘] LC in an equally spaced grid of size nxn. The

fimction 70'( y) is nothing but an arithmetic average of p] and in view of (3.7), the elastic

tensor
EA(y)=/7(y)E° (3.12)
is an arithmetic average of E (y) over the substructure.

A substructure made of material E A would over-estimate the stiflhess of the original

substructure, as EA (y) is an upper bound for E (y). To compensate, we shall look for a

harmonic average of E (y). For this purpose we define the firnction

35

whe

coar

by l

CONE

isal

used

a11d I]

This E

a fine

N—l

950’): Z gilt <PJ,k1(y) (3-13)
k,l=0

where 0,{, i l/kal . A multi-resolution analysis of this fine-scale fimction produces a

coarse representation

n—l
9 (y) = Z: 91d $11101) (3.14)
k,l=0
from which we can define
n
20) = 2 gm can (y) (3.15)
k,1=0

by letting Bk! i 1/ 5k]. The firnction p is a harmonic average of p] and the

corresponding elastic tensor is then

EH (y)=p(y)EO (3.16)
is a harmonic average of E (y) over the substructure. A substructure made of material

E H would under estimate the stiffness of the original substructure. E A and E H are

used to define an efi‘ective material for substructure QC as the weighted average
Ef- =aEA +(l—a)EH, a€[0,1] (3.17)
and the corresponding material distribution function p; (y) as
Pi(y)=afi(y)+(1-a)e(y), 06M (3-18)

This average approximates the material properties of the substructure when reduced from

a fine level J to a coarse level j < J . Using these properties and a finite element

36

discretization of equally spaced elements of size A = 2’1 LC an effective stiffness matrix

K5- of dimension 2(n+1)2 x2(n+l)2 can be constructed by assembling n2 standard

finite element matrices per substructure, i.e.,
K3. : Z (p;)“k° (3.19)

Here, [to is the stiffness matrix of a four-noded square finite element made of material
E0. The reduced substructure can be interpreted as a super-element whose stifiiress
matrix K3- can be used as a building block in an assembly of other substructures

following the usual rules for assembly.

3.4 A Mum-Resolution Analysis of the Displacement

The matrix K j constructed in the earlier section (3.19) is based on a multi-resolution

analysis of the material distribution alone and thus it is only an approximation to the
effective properties of a substructure. This approximation is surprisingly accurate and is
sufficient for problems with simple nricrostructures. However, a more consistent analysis
may be required when detailed representations of the small-scale effects are of interest.
One such analysis is based on numerical homogenization using a multi-resolution

analysis of the displacements.

3.4.1 A Periodic Mum-Resolution Reduction Scheme

The variational form (or weak form) of the plane elasticity problem defined on a single

periodic substructure discretized at its finest scale can be written as

37

fE(y)e(u)e(v)dy= ffvdy, Vu,vEVp(Qc) (3.20)
O

9c c

where
Vp(nc)= {uEH1(QC):uich—periodic} (3.21)
and the body force f has mean-value of zero in (26. This problem is solved using a

wavelet-Galerkin technique, outlined briefly here (for further details about the technique

refer Diaz [14], DeRose and Diaz [1 l] and the references therein). The tensor of elastic

properties E (y) and the trial and weight fimctions u, v are approximated as follows:

2-’ —1
Ed): 2: Ed vino)

k,I=0
2’ —1

ad): 2 ad some) (322)

k,l=0
2’ —1

”(y) = 2 Va WHO“)

k,l =0
where the coefiicients Elf, , u]{, and vi, are wavelet coeflicients associated with pixel

(k,I) in a J -level discretization of substructure QC. The material tensor is approximated

using the piece-wise constant Haar scaling functions cpflaa' (y) as the fimction E (y) has
to satisfy only minimal continuity requirements whereas the displacement functions nwd

to satisfy suflicient smoothness conditions. Daubechies (D6) fimctions 90( y) are used to
approximate displacements. Using the approximations (3.22) into (3.20) yields a

2N2 x 2N2 linear system of equations

LJIIJ ZfJ (3.23)

38

where the mesh size that resolves pc is N x N , N 2 2J and L is of the form

L? L?
= yx W (3.24)
LJ LJ
In equation (3.24)
xx _ J,ll L0,l,0 J,33 0,1,0,l
LJ,klmn — ZEN Cqulmn +qu Cqulmn
P q
)9! _ J,33 0,1,1,0 ,12 l,0,0,l
LJ,klmn -2qu Cqulmn +EPI‘I Cqulnm
1” q (3.25)

yx __ J,33 l,0,0,l ,12 ,l,l,0
LJ,kImn _ZEP¢I Cqulmn+Eliq Coqulmn
p,q

11);)“ :2 J,22 ,l,0,l +Egzl33CL0,l,0
, mn

 

 

P9 qulmn qulmn
P")
where
' .11 ,12 ‘
EA, Big 0
J _ ,12 ,22
Em — 2;, m 0
,33
0 a. 4
The terms

Catti= x

 

 

f ‘p-IZP (xlDo‘cpiJc 0005905,», (1061*
Q

 

 

f¢£q(y)03905,1(y)09.05,n (y)dy
Q

are called connection coeflicients and depend only on the choice of basis functions. Di
denotes a derivative operator (several algorithms are available to compute these

coefficients, see Kunoth [26]). L J is the stifl‘ness matrix associated with a material

distribution pc discretized at the finest scale Sc. It is the most detailed representation of

39

the substructure but it is too large to be useful in computations involving many

substructures. We reduce L J using a multi—resolution analysis that splits “J into

coarse-scale and fine detail components 6 and ii respectively.

Consider the system represented in equation (3 .23) and let the associated space of

functions be VJ, i.e., LJ:VJ—>VJ and “J,fJEVJ. Define W to be the

transformation

w : VJ _, VJ”l EBWJ_1 (3.26)
that maps a function in the space VJ into functions in the space VJ—1 and its orthogonal
complement WJ_1. The orthogonal transformation

iJ—1
iiJ—r

PJ—luJ

WIIJ =
QJ-1“J

(3.27)

 

 

 

 

represents the splitting of a vector “J into its coarse scale component i J_1 and the

details frJ_1, where the operators PJ_1 :VJ —+VJ_1 and QJ_.1:VJ -> WJ"l are the
coarse and detail projection operators respectively. The equation (3.23) can now be
expanded as follows:

LJ-r CJ—l E4

I'J—r

5.1—1
fiJ—r

WLJwT(WuJ) = (3.28)

 

 

 

 

 

 

T
CJ—r AJ—r
where,

AJ—r == QJ—ILJQJ—l I WJ_1 —* WJ_1
CJ—l = PJ—rLJQJ—r 3WJ_1 —’ VJ—l (3-29)

LJ—r = PJ—rLJPJ—l ZVJ—l —+ VJ_1

40

Here 6-1 and fJ_1 are the coarse-scale and detail components of the external force.

The coarse scale component of the displacements iiJ_1 is found by block Gauss

elimination, which yields the equation
—1 T ._ - -—1 ~
(LJ—r —CJ-1AJ-1CJ—1)“J-1 = fJ—1_CJ—1AJ-lfJ—l (3.30)

For the class of operators of interest here it can be shown that the operator AJ_1 is
indeed invertible. This can be proved for our prototype problem in 2D elasticity as
follows (The result can be extended to other cases easily). A 2D elastic stiffness matrix
L J is positive semi-definite with two zero-eigenvalues that correspond to the two rigid
body modes (translation in the x-direction and the translation in the y—direction). These

two modes are constant functions and can be represented exactly at any scale, i.e., the

detail components associated with these modes are always zero. This means that the

operator Q J_1 is orthogonal to these rigid body modes. Let 8 denote the subspace
spanned by the rigid body modes. From the positive-semi-definiteness of L J we know,

that for all non-zero vectors x in the space V J that are orthogonal to the rigid body

modes
xTLJx >0, x:(x¢0,xEVJ\6} (3.31)
Using the orthogonality of the wavelet transform operator, (3.31) can be written as
(Wx)T (WLWT)(WX)>O, x:{x¢0, xEV‘I \9} (3.32)

Using equation (3.28) this can be expressed as

T LJ—r CJ—l

Wx T
CJ-r AJ—r

Wx

 

 

 

 

>0, x:{x=:0,xEVJ\9} (3.33)

 

 

41

Also, any vector in 6 has no detail components, i.e.,

VxEG, Wx=

 

x (3 34
0 ‘ )

0
Consider yEVJ such that Wy=[ l, where vEWJ-l. It can be seen fi'om equation
v

(3.34) (and using the linearity of the concerned operators) that y is orthogonal to any

vector in the subspace 9. Now, for all such vectors y ,

yTLJy > o (3.35)

Using the same approach as in equation (3.33), we have

T LJ—l CJ—l 0

T
CJ—r AJ—r

0

V

> o, Vv 6 WJ‘1 (3.36)
V

 

 

 

 

 

 

Alter carrying out the multiplication we have
vTAJ_1v > o, Vv 6 WJ—1 (3.37)

Thus, fi'om equation (3.37) we can see that the matrix A J_1 is positive-definite and thus
invertible. More generally, Engquist and Runborg [17] show that the operator AJ_1 is
bijective for a class of elliptic operators L J obtained from bilinear forms.
The matrix

HJ—l i LJ—r — CJ—iAilrci—r (3 38)
is the effective stifliress matrix at level J —1 associated with a periodic patch of identical
substructures QC discretized at level J . Further reductions are possible by applying the

reduction procedure recursively to obtain a sequence of effective stifl‘ness matrices

HJ_2 , HJ_3, etc. It should be noted that the each reduction of the stiifiress matrix

42

results in a matrix that is relatively much denser than the original matrix. This is

illustrated in figure 3.2 where the dark regions denote non-zero entries.

   

LJ

Figure 3.2: Structure of LJ and HJ_k
If the structure considered involves forces that are slow-varying in nature, i.e., f = 0,
then the above set of matrices (L J,{H J_k }) form a complete description of the

substructure at various scales.

The obtained efi‘ective stifliress matrices operate on vectors that result fiom a
wavelet discretization of the associated functions. This creates some difficulties not only
in the application of boundary conditions but also in the assembly of stifi‘ness matrices
from several substructures. Moreover, the strain energy of an assembly of substructures
is not the sum of strain energies of individual substructures since the associated basis
functions in each substructure extend beyond the boundaries of the substructure. This is
due to the fact that the support (=5) of the considered scaling fimctions (D6) is greater
than 1. It might be advantageous to convert these wavelet efl‘ective stifi‘ness matrices
into equivalent nodal stiffness matrices that operate on the nodal values of the associated

forces and displacements defined on the substructure (similar to a bi-linear finite-element

43

stifiiress matrix). Also, these nodal-matrices would be more portable and can be
incorporated into general-purpose widely available finite element codes. This conversion

is described next.
3.4.2 Transformation to Finite-Element Form
This section deals with finding an approximation to 111- that acts on finite-element

instead of wavelet spaces. Introduce V” and V” as finite-dimensional spaces of

periodic functions on DC spanned respectively by D6 scaling functions (cp) and bi-linear

finite element shape functions (cph) shown in figure 3.3.

   

'2

4
A D6 Scaling Function (2D) A 2D Bi-linear (hat) Function

Figure 3.3: Basis functions in the wavelet and finite-element spaces respectively

Representative functions in these spaces are of the form

n-1 n-1
h h
uw = Z ulcisokz and uh = Z um (3.39)
k,l=0 k,l=0

where n=2j . Here it” ={ui’1} are wavelet coefficients and uh ={ullr'1} are

displacements at nodes on an nxn uniformly spaced grid with spacing A: 2—j LC.

Three approaches to convert a periodic stiffness matrix, H J- , into an equivalent nodal

stiffness matrix in the (bi-linear) finite-element space are presented.

3.4.2.1 Method I

For a given periodic stifliress matrix H (corresponding to a periodically tiled domain)
and a force f w in V w with mean value zero (i.e., orthogonal to the two rigid body

modes) there exists a unique displacement u” E V w (and associated coefficients n”)

such that

Hu‘” = r” (3.40)
Let f h be the orthogonal projection of f " onto Vb computed as

r” = Phr’” (3.41)

where Ph is an orthogonal projection matrix that depends only on the choice of the two

basis fimctions (as described in chapter 2, section 2.7). We now look for a transformation

matrix Q that transforms the vector of wavelet displacement coefiicients uw , into a

vector of nodal displacements uh , in such a way that the work of the loads 1'" and fh is
the same, i.e.,
u” = Qu‘” (3.42)
and
rWTu‘” = rhTu” (3.43)
This is illustrated in figure 3.4.

45

 

Figure 3 .4: Conversion from wavelet to finite element spaces using method I

The work due to the force and displacement in V h can be written as

thuh = fWTPhTQuw (3.44)

Combining equations (3.43) and (3.44) we see that
—1
Q = [PM] (3.45)

Now, we look for a matrix E that relates the forces and displacements in Vh as

Ku” = r” (3.46)

Using equations (3.42) and (3.45) in (3.40), we have
-—1
Eu" = r" => HPhTuh = [Ph] r” => PhHPhTuh = r” (3.47)
Thus the equivalent finite element matrix is

17: = PhHPhT (3.48)

3.4.2.2 Method II
Here we start with a finite-element firnction (force) f h E V h and denote f w to be its

orthogonal projection onto V” who’s coefiicients are computed as

46

r” = owh (3.49)
where P” is an orthogonal projection matrix. The objective here is to find a

transformation matrix R : V h —+ V w that transforms a finite-element displacement into an

equivalent wavelet displacement and in turn the associated nodal stiffness matrix such
that the work of the loads 1'” and f h is the same, i.e.,

u‘” = Ru” (3.50)
and

thuh = rWTu‘” (3.51)

This process is illustrated in figure 3.5.

 

Figure 3.5: Conversion from wavelet to finite element spaces using method 11

The work associated with u”, f w E V W can be expressed as

{”11” =rhTPWTRu” (3.52)

Combining equations (3.51) and (3.52) we see that the works are the same only if

R =[PWTJ—1 (3.53)

47

As before, we look for a finite element stiffness matrix that relates uh and I” as

in” = r” (3.54)

Using equations (3.50) and (3.53) in (3.40), we have

-1
Eu” = r” => 11hr”) uh = owh

_1 -1 (3.55)
4(1)”) HlPWT] uh =r"
The equivalent nodal stifliress matrix is then
— w —1 WT —l
K =[P ] H[P ] (3.56)
3.4.2.3 Method III
Here we introduce a matrix B , defined such that
{a ( .) 1")
— y) e u e u dy (3.57)

0c
The right hand side of equation (3.57) is the elasticity bilinear form with finite-element
(hat) trial functions and wavelet weight functions (compare to 3.22). The elastic tensor

E is an average of the material in the substructure and is assumed to be of the form
E=p3(y)E° (3.58)

where E0 is a reference elasticity tensor and p3 is a distribution of relative densities that

is obtained by averaging the fine-scale distribution of relative densities to the coarse-scale
using (3.18) (It should be noted that this is the only place where an average of the
material is used in this method that is based on a multi-resolution analysis of

displacements). The matrix B is one that transforms a finite-element displacement

uh E Vh into a wavelet body force f W E V w, i.e., f w = Bub is the wavelet body force

48

tl

{1'3

WC

and

that corresponds to a finite-element pre-strain 5(uh). If the basis firnctions used as the

trial and weight fimctions are interchanged, matrix BT is the one that transforms a given

(wavelet) displacement uw E V” into a finite-element body force f h E V h .
Having defined this operator B , we look for an operator Q : V w —’ V k that

transforms a given wavelet displacement u” E V" into a nodal displacement uh E V h

such that the work due to the transformed forces and displacements is the same as the

work of the corresponding firnctions in V” , i.e.,

uh = Quw (3.59)

and

thuh = rWTu‘” (3.60)

This process is illustrated in figure 3.6.

 

Figure 3.6: Conversion from wavelet to finite-element spaces using method H1

The work corresponding to the displacements and forces in V h can be written as

f ”uh = uWTBQuw = f WTH-IBQuw (3.61)

49

From eql

Using eQI

Thus, the

Tl
equivaleni

substrucn

3.4.3 Cc

In this se
matrix of

1- A:

 

sul

(di

 

From equations (3.60) and (3.61) we see that

Q = B-lflj (3.62)

Using equation (3.62) in the equation I—(uh = f h we have,
fin" = r” 4:) Ron" = 13%” e» I‘m—Inn” = BTu‘” (3.63)
Thus, the equivalent nodal stiffness matrix is then given as

R = 3711“]; (3.64)

The matrix K obtained using either of these approaches is the “finite element
equivalent” to Hj. It relates nodal degrees of fieedom to nodal forces in a periodic

substructure reduced from level J to level j .

3.4.3 Computational Aspects

In this section, the computational procedure involved in obtaining a reduced stifliress
matrix of a substructure using the multi-resolution analysis of displacements is outlined.

l. Assemble a fine-scale wavelet stimiess matrix of a substructure, L J , as

NJ
L J = Z p‘gl0 (sum interpreted in the sense of assembly)
e=l

where N J 2 2J x 2J is the number of pixels in the fine-scale discretization of the

substructure, l0 is a pre-computed wavelet “element” stiffness matrix of a pixel

(dimension 50 x 50 for D6 scaling function) with a reference material tensor and

p! is the value of the relative density firnction at a pixel e.

50

For]:-

.N
0

C1

DJ

at

Tl

ha

end [00p

de

an.

 

Usi

Var

 

For j = J to J —k +1 (k is the number of reduction levels), do

2. Compute the wavelet decomposition of the stiffness matrix at level j

 

 

 

N2-

1

2 W21) 2
H - 2N - 2
'CI-l 1“ T

 

 

 

 

 

 

 

3. Compute the Schur’s complement to obtain the efi‘ective reduced stiffiress matrix

at level j -l

-1 T
Hj—l = Lj—l _Cj—IAj—le-l

J

. . . 3N3- 3N2- . . N}- .
This mvolves the solution of a TX T system of equatrons With —2— right

2 2
N,- 3N]-

2
3N-
_x_ J

2
N.

2

hand sides and the multiplication of a matrix with a
matrix.

end loop

n
4. Assemble B= prgbo , where n: ZJ-k x2‘]—k , p3 is an averaged relative
e=l

density distribution at the reduced scale (obtained as shown in equation (3.18))

and b0 is a pre-computed “element” matrix of size 50 x8 that is constructed
using D6 and Hat firnctions as the weight and trial fimctions respectively in the

variational form of the elasticity equation.

5. Compute I—( = BTHjlkB

51

End.

The mo
This my
as matrix
stiffness 1
computat
finest-sea]
reference

Steps 4 ar

are much

fine~SCale

relatively ,

 

The matrix H J_k is positive semi-definite (it has two rigid body modes). It is
inverted by adding a matrix 50 to remove the rigid body modes. Here, 5 is a

1 1 o 0
suitable scalar penalty and Q=pr, where p: 0 0 1 1 is a

matrix of rigid body modes.

End.

The most computationally intensive step is computing the Schur’s complement (step 3).
This involves the solution of a system of equations for multiple right hand sides as well
as matrix multiplication. At the end of each reduction stage (steps 2 and 3) the size of the
stiffness matrix (and hence the number of equations) reduces by a factor of 4 and thus the
computations progressively reduce after each reduction step. The first step, involving the
finest-scale stiffness matrix, can be implemented efficiently with a pre-computed
reference “element” stiffness matrix and using sparse assembly. The computations in

steps 4 and 5 are performed on matrices in the reduced-scale; the sizes of these matrices

22k

are much smaller (factor of , where k is the number of reductions) compared to the

fine-scale stiffness matrix. It should be noted that the reduced stifliress matrices are

relatively much denser than the fine-scale matrices, as illustrated in Figure 3 .2.

52

3.5 C

In this :

respect

3.5.1 l
displa
C onsidc
to an h

the spa

nxn‘ i
matrix

KE 2 t
diagong

basesf

eigenve

3.5 Comparison of the Model Reduction Schemes

In this section, comparisons of the proposed model reduction schemes are presented with

respect to each other and with the classical homogenization scheme.

3.5.1 Comparison between MRA of material distribution and MRA of

displacements

Consider K : V —+ V and K : V —> V to be the reduced stiffiress matrices that correspond
to an MRA of material distribution and an MRA of displacements, respectively, and V is
the space of kinematically admissible functions. Let the dimensions of the matrices be
nxn, where n is the number of degrees of freedom in the models. Define E to be a
matrix of eigenvectors and Q to be a diagonal matrix of eigenvalues of K , i.e.,

KE = 9E. Similarly, define I}. to be the matrix of eigenvectors and fl to be the

diagonal matrix of eigenvalues of K , i.e., Ki) = {If}. Since E and i) are orthonormal

bases for V , any eigenvector e,- of K can be expressed as a linear combination of

A

eigenvectors (3;) of K , i.e.,
f

The two stiffness matrices have two zero-energy (rigid body) modes that correspond to a

translation in the x and y directions respectively. Let i) be the matrix of eigenvectors of

K other than the rigid body modes, dim (E) = n xn —— 2 . This can be expressed as

A

F: = EB (3.66)

f1

5?

as

DU]

where the entries in matrix B are of the form

Bi- =(a;,éj) (3.67)

and E,- and E,- are colurrm vectors of I73 and B respectively.
Any (displacement) vector u that is orthogonal to the rigid body modes can be written as,
u = Ea (3.68)
for some ai E R . Using the orthogonality of Ii) , we can write the strain energy
associated with this displacement it using the material MRA model as
uTKu = «TETKEa = (ITS-ill (3.69)
where Q is a diagonal matrix of the non-zero eigenvalues of K. The displacement u
can be expressed in terms of the eigenvectors of K as
u = 13:13.: (3.70)
where the entries of B are as defined in (3.67). Now, the strain energy associated with u
computed using the MRA of displacements can be expressed as
uTKu = «TBTETKEBa = 6713712136 (3.71)
The relative error in strain energy obtained using the two models can then be expressed

as

IuTKu — uTKu laT [fl — BTSAIB] a

 

 

 

 

— _ (3 .72)
luTKu laTflal

 

It can be shown using Rayleigh’s principle that the eigenvalues of the matrix in the

numerator of the right hand side of equation (3.72) form an upper bound on the relative

54

error of

rigid bo«

Here, th

five difi?

Figure 3

obtained

eight)

(All =

 

 

error of the strain energies, i.e., for a normalized displacement u that is orthogonal to

rigid body modes,

 

luTKu—uTKul __ l T .
max T , *Smax eig(I—ST B QB) (3.73)
MS; '11 Kul
Ila:

Here, the reduced stiffness matrices of a periodic tiling of a substructure for a choice of

five difi‘erent material distributions (Figure 3 .7) are computed and compared.

(a) (C) (6)
Figure 3 .7: Choices of material distributions for the comparison

 

Figure 3.8 shows the plots of the eigenvalues of the matrix (A = I —§-1BT§B)

obtained for the various choices of the material distributions chosen.

 

 

 

 

 

-A- (a)
as .
O V
0.25 -v- (d) ' 4
-o- M
a 0.2L 711
a ' ‘- ‘ ...‘ ~ .-
'5 ‘.
11
5: 0.15 -

0.1

0.05

 

 

 

l

 

Figure 3 .8: Bounds on the maximum relative error between the strain
energy computed using material MRA and displacement MRA

55

The r 6
spacing
are 121
usually
spatial

error i1
combin

onthel

3.5.2 1
Home

Here,

displace

mfinites.
15 assun

base cell

Where

represent:
Wording,I

 

The reduced stiffness matrices are at a discretization that corresponds to a uniform
spacing of degrees of freedom in a 8 x 8 grid. The dimensions of these stiffiress matrices
are 128x128. The figure 3.8 shows the first 64 eigenvalues of the matrix A, as we are
usually not concerned with the higher energy modes that correspond to rapidly varying
spatial displacements. The value of the plots at each ‘i’ corresponds to a bound on the
error in the strain energy associated with a displacement that can be expressed as a linear
combination of the first ‘i’ modes. It can be seen that for the first few modes, the bound

on the error is small but it increases as more number of modes are included.

3.5.2 Comparison between MRA techniques and Classical
Homogenization
Here, the averaging schemes based on the MRA of material distributions and
displacements are compared to the classical homogenization method.

In classical homogenization, the material distribution is assumed to vary at
infinitesimal scales and an asymptotic analysis is used to compute effective properties. It

is assumed that the material distribution is characterized by the periodic repetition of a

base cell (1’ ). The displacement field it is expanded in an asymptotic series

u=uo(x,y)+€ul(x,y)+62u2(x,y)+~- (3.74)

where
y = 16/8 (3.75)
represents the local (microscopic) coordinates and x is the global (macroscopic)

coordinate. It can be shown (see Bensoussan et al [5]) that the effective material tensor is

given as

56

where \

periodic 1

The straill

Where E

50. Thu

$1

a periodic

Equation

 

 

k1

 

1 3x
EW— _ Y —f 5W —E,-qu-—’i- dY (3.76)
Y ay‘l

 

where xfg’ is the solution to the so-called cell-problem defined on the infinitesimal

periodic characteristic cell, given as follows

8
faqu— -—X—” div—tor: fEiJ-Hade VvEV'
6y 6y 6y,- (377)
V: {vzv rs Y-periodic}Y
The strain energy form of (3.76) can be written as
_H 50 80 _1 o * 0 *

Y

where E is the material tensor and e. is the resultant strain due to the applied pre-strain,

50. Thus, we can define the strain energy associated with a (finite size) substructure QC ,

corresponding to a periodic repetition of infinitesimal cells Y , as

(PH g -}17£(Eiqu (€3- —e;-)(egq —e;,q))dY theas(Qc) (3.79)

Similarly, the strain energy associated with applying a constant pre-strain (so) on

a periodic tiling of a (finite-size) substructure QC can be defined as
t t
(I) Q f(E,-qu (83- ‘51)“qu -epq))d§2 (3.80)
9c

Equation (3.80) can then be written as

‘1’: IE iqusije pq +fE queijem 2f}; 0109595 N (3°81)

57

The ela:

substruc‘

Using (3

The first

strain an
prescribe
second 1
displacen

The strair

Where K

bang that

 

The Strair

Where R

SCaJe SOlu

 

The elasticity problem associated with the application of a constant pre-strain on a

substructure no, is defined as: Find if, such that
1 mp.) .11: is. o-eu*r.1) so

Using (3.82) in (3.81), the strain energy becomes

<I>=fEiqu 51]qu —fE,~J-pqe;je W (3.83)

The first term on the right hand side of equation (3.83) does not involve the resultant

strain and can be computed exactly without solving any elasticity problem (as 80 is a

prescribed constant strain). The model reduction schemes are used to compute the
second term using an approximation of either the material distribution or the
displacements at a coarse scale, i.e., MRA of material or MRA of displacements.

The strain energy using a multi-resolution analysis of material is then

M *T "
<I> =fE ”megs pq —u Kcu (3.84)

where KC is a reduced stiffness matrix obtained as shown in (3.19), the only difference
being that the substructure involved is assumed to be periodic.

The strain energy using a multi-resolution analysis of displacements is

D *T- *
ch =ch megs m —u Ku (3.85)

where K is a reduced stiffness matrix obtained as shown in (3 .64) and u‘ is the coarse-

scale solution to a constant pre-strain elasticity problem (It should be noted that the body

58

force (T.1

 

material
T
homogerj
various 1
these me]
to a pres
a particu
would be

material .

 

material

directioni

Strains 5‘

 

force (right hand side of 3.82) in this case is computed using an equivalent, reduced
material tensor).

The comparison between the multi-resolution reduction schemes and classical
homogenization is carried out by applying constant pre-strains to substructures with
various material distributions and comparing the resulting strain energies when either of
these models is used to compute the strain energy. The stifl‘er material, when subjected
to a prescribed constant pre-strain will result in higher strain energy. In order to say that
a particular material is stifl‘er, this result has to hold for all possrble strains. However, it
would be meaningful only to consider strains whose principal directions coincide with the
material axes, as only then would the strain energy be maximmr, i.e., the best use of the
material is only when it is oriented in such a way that it coincides with the principal

directions of the applied strain (see Pedersen [31]). Thus, among all such (normalized)

strains 80 =

 

 

, the strain energy produced by the strain with ,8 = O is the highest.
77

Here, substructures of unit sizes are considered with five choices of material distributions
and the strain energies as a result of applying constant pre-strains of the form

0 l

5‘:

 

0
1, where r) E [—1,1], are computed.
77

Figures 3.9 (a) to (e) show the plots of strain energy (<I>) versus the parameter 7) ,

computed using the various methods discussed here.

59

2.00 Ti,
1.80 «
1.60 «
1.40 .
1.20 -
1.00

 

 

 

 
 

 

 

 

 

 

 

0'80 - + Fine-Scale
05° " —I— Homogen-tion
0.40 + + Disp MRA j=3
0.20 ~— -x-— Mat MRA j=3
0.00 I r n

 

-1.00 -0.50 0.00 0.50 1.00
Figure 3 .9(a): Strain energy comparisons for material distribution (a)

 

 

 

 

   
  

2.00 - (I)

1 8° - + Fine-Scale

1 60 -I— Homogenization
' + Disp MRAj=3

“'0 ‘ —x— Mat MRA j=3

1.20 d

1.00 .

0.80 -

 

 

 

0.60
0.40
0.20

now r I T 7 l n
-1.00 -0.50 0.00 0.50 1.00

 

 

Figure 3.9(b): Strain energy comparisons for material distribution (b)

60

 

 

 

  
 

 

 

 

 

 

0.80 + Fine-Scale
0.60 - —I— Homogen-tion
0.40 q + Disp MRA i=3
0.20 q -—>(—— Mat MRA i=3
0.00 . 1 . . ’7

 

 

-1.00 -0.50 0.00 0.50 1.00

Figure 3 .9(c): Strain energy comparisons for material distribution (c)

 

2.00 “ Q

1.30 i + Fine-Scale
1.60 — —I— Homogenization
1_4o . + Disp MRA F3

1 20 J —)(— Mat MRA i=3

 

 

 

 

 

 

    

0.00 r 1
-1.00 -0.50 0.00 0.50 1.00

 

Figure 3 .9(d): Strain energy comparisons for material distribution (d)

61

hrthe at
discretiz
anmnnc
energy'i
xmemr
uxd F
propertic
dumm
underesfi
complica
loss of ;
essentiall
matrices

the effect

3.6 (

Matri:

The redl
DOt 3’61 t
build a r

be exp!

\

In the above plots, the reduction process involved starting fi'om a 64 x 64( J=6 ) element
discretization and reducing it to a 8 x8 ( j=3 ) discretization. For simple material
structures such as the perforated substructures in Figure 3.9 (a), (b) and (c), the strain
energy in the reduced model is essentially indistinguishable from the energy in the filli-
scale model, regardless of whether the material based or displacement based MRA is
used. For comparison, the graphs also show the energy in a substructure whose elastic
properties are those of a periodic mixture at infinitesimal scales, i.e., the result from
classical homogenization methods. In all cases, using such efl‘ective properties will
underestimate the strain energy (and hence the stiffness) of the substructures. For more
complicated structures, such as in Figure 3.9 (d), the loss of resolution results in some
loss of accuracy. However, even in this case the two reduction processes result in
essentially identical approximations of the strain energy and the reduced stifi‘ness
matrices still provide a more accurate of the strain energy in the firll-scale structure than

the effective properties obtained from classical homogenization.

3.6 Computing the Non-Periodic Reduced Stiffness

Matrix of a Substructure

The reduced matrix I? , obtained using a multi-resolution analysis of the displacements is
not yet ready to be used as a super-element and assembled with other such matrices to
build a model of a complex structure. This is because, while the area of the structure can

be expressed a union of smaller areas that correspond to smaller substructures,

Q = UCQC , the total stiffness matrix of the structure cannot be expressed as an assembly

62

of pen}

interpret
that cha:
For this
different
Ieclmiqu
a single
surround
substruct
layout is

3.10.

 

Fic

 

of periodic stiffness matrices of individual substructures, K(Q)¢Zl_(c(flc) (sum
C

interpreted in the sense of assembly). However, I_( can be used to construct a matrix K
that characterizes a substructure as a single entity in a non-periodic arrangement.

For this we turn to the well-known computational scheme used in the solution of partial
differential equations on arbitrarily shaped domains called the fictitious domain
technique, (e.g., see Bakhalov and Knyazev [1], Glowinski et al [21]). The properties of
a single substructure that is not part of a periodic arrangement can be obtained by
surrounding the substructure by weak material of a sumciently low strength so that the
substructure is essentially unaffected by its surroundings. A periodic arrangement of this

layout is characterized by a fictitious substructure QFc =Qc U000 as shown in figure

3.10.

n U
n n
U u

 

 

Fictitious Substructure QFc Periodic arrangement of fictitious substructures

Figure 3.10. Substructure surrounded by a weak fictitious domain

63

For exar

fine—soak

the fictiti

Thus. th
substruct

to level

substruct

Using thi

Where B

Recall rh
the COTTe
OV'er the
SubStrUct

3 See Sec.

 

For example, if the material distribution within 0, is resolved by 2’ x 2’ pixelss, the

fine-scale problem in “Pa involves 2‘]+1 x 21+l pixels and the material distribution in

the fictitious substructure is defined as follows

U €<<1 iajEQOC

Thus, the fine-scale problem on a substructure is now extended to one on the fictitious
substructure and the reduction procedure discussed eariier is performed from level J +1

to level j +1 to compute a reduced stiffness matrix corresponding to this fictitious

substructure, denoted as H i +1.

Using the approach in section 3.5.2.3, we compute the nodal matrix as

—1
—F F F F
K 1+1 :3 1+1 [11 1+1] 131-+1 (3.87)
where Bf“ is defined as
uWTBfHuh = fEfH (y) €(uw)s(uh)dy (3.88)
0

Recall that when this method was illustrated with periodic stiffness matrices and vectors,
the corresponding material tensor used in the definition of the matrix B was the average
over the substructure (as given in (3.58)). In the case of a fictitious substructure (i.e.,
substructure surrounded by a weak fictitious material) the average material tensor at scale

j+l is defined as

 

§ See section 3.2

This is c
original
resulting
The pres
analysis
the actuz

both the?
0f the fit
is definet

state of c

 

 

Where (

dOma'm E

058 (y)

The Den

in fight".

)= p§3(y)E0 if yEQc

E . (
1'1 y .

(3.89)

This is done so that the terms corresponding to the non-periodic stiffness matrix of the
original substructure (without the fictitious domain) are the only non-zero terms in the
resulting matrix obtained using (3.87) and thus can be directly obtained.

The presence of the fictitious domain in the reduction process using the multi-resolution
analysis of displacements causes the terms that are associated with the boundary between
the actual substructure and the weak fictitious domain to have properties that are due to

both these materials. These edge eflects need to be compensated using a judicious choice
of the filnction p§3 , which is the only adjustable parameter in this process. Here p53 (y)
is defined such that a patch of 3 x 3 isotropic, homogeneous substructures reproduces a

state of constant strain exactly. We assume that Rig (y) is of the form
PiB(y)=a§B(y)p§ (y) (390)

where 053 (y) is a correction factor to compensate for the presence of the fictitious

domain and p; (y) E] for a homogeneous substructure. The piecewise constant function

033 (y) is of the form

n
0‘58 (y) = Z 01:190ka (y) (3.91)
k,l=0

The patch is subjected to three prescribed constant strain fields 4,551,531 as illustrated

infigure3.ll.

65

3....1-C-I1-I-I
I I

---l

 

 

 

_____ ____ ___.. ... ’2: ::flP--r----’

 

 

 

 

 

 

--1--1--1
I
I
I
I
I
nun-uI-II-n'

I

I

I

I

I

In...
~

‘-
‘~
s
f.

 

 

 

 

 

 

 

 

 

 

 

 

 

--- -————— ———— ———— II-J

 

L... ———— ———— —————

Figure 3.11. Constant pre-strains applied to a patch of substructures

The coefiicients a“ are chosen to minimize the function

       

2
s” rc+| ”1 —e({” I (3.92)

I 11

where e , e 51”

and are the strains resulting from the imposed pre-strains.

K5511 is a matrix of size 2(4n2 x4112) for n = 2} , fi'om this we extract a sub-matrix K3-
of dimensions 2((n+l)2 x (n+l)2), in the case of approach 1]], these are the only non-

zero terms in K 1+1 This matrix characterizes the behavior of a substructure QC as a

single entity in a non-periodic setting when reduced fi'om level J to level j < J . This

can now be used as a building block in an assembly of other substructures to construct a

complete reduced model of a structure.

3.7 Computing Fine-Scale Stresses By Augmentation

The solution to the reduced system of equations obtained using either of the earlier
described methods represents the coarse behavior of the original system. This is useful

and may even be suficient for certain design problems that involve only large-scale

66

behavior. However, it may be necessary to compute some small-scale effects at certain
locations (substructures) of interest, such as local stresses at locations of possible stress-
concentrations such as sharp corners or abrupt changes in geometry or material
properties. In this section, a method to compute the small-scale stresses (where
necessary) given the coarse solution is presented.

This problem can be compared to that of computing the small-scale stresses in classical
periodic homogenization. In the case of classical homogenization, the large-scale
(coarse) solution is used to compute a coarse-scale strain field. Since this method
assumes that the small-scale is microscopic, the coarse-scale strain field at any point in
the domain may be thought of as a constant strain over the (infinitesimal) region occupied
by the microstructure. The small-scale displacements in the case of classical
homogenization are computed as the solution to the cell problem (3.77) for a constant
pre-strain, where the value of constant pre-strain is the value of the large-scale strain at a
location of interest. The important difi‘erence in the present scheme is that the small-
scale problem is solved on a domain of finite size rather than on an infinitesimal cell.
Thus the small-scale displacements cannot be computed simply by applying a constant
pre-strain over the substructure, as in classical homogenization. The large-scale strains in
the substructure are in general not constant. In the case of the model reduction using a
multi-resolution analysis of the displacements, it is possible to obtain the fine-scale
displacements at any substructure using a consistent augmentation procedure. This is not
possible in the case of the multi-resolution analysis of the material.

Let c be a particular substructure of interest whose small-scale stresses need to be

computed and let uj E V; be the large-scale nodal displacement in this substructure

67

obtained fiom the solution of the reduced system. The multi-resolution analysis

discussed earlier assumes periodicity of the associated fimctions in their respective

domains and also deals with functions in a wavelet basis (V w ). Thus, in order to use the
multi-resolution procedure to compute the fine-scale displacements, we first need to find

an equivalent, periodic large-scale nodal displacement function and convert it into an

. . - W
equlvalent filnctlon m V j .

3.7.1 Computing the periodic large-scale nodal dlsplacements

It is assumed that the fine-scale strain in a substructure, refine , can be expressed as the
sum of a periodic fine-scale strain in an equivalent periodic substructure and a large-scale
non-periodic strain, i.e.,

Efine = eljjne + Ex/olgrrse (3.93)

where the subscripts P and NP denote periodic and non-periodic, respectively. This
assumption in fact means that the non-periodic component of the strain is large-scale in
nature. In addition, it is assumed that the large-scale displacement in the substructure,
u j , can also be expressed as a sum of a periodic component and a non-periodic
component, i.e.,

P NP
This assumption is illustrated in figure 3.12, where a typical non-periodic displacement

filnction is shown to be composed of a non-periodic component and a periodic

component.

68

 

 

 

 

 

A periodic function (uf)

21> 619

 

 

 

 

 

A non-periodic filnction (uj)

 

 

 

A non-periodic fimction (14f? )

Figure 3.12. Decomposition of a general non-periodic fianction into a coarse,

The coarse-scale strain energy in the substructure can be expressed as
UC=%S{Eje(uj)s(uJ-)dy (3.95)

Using (3.94) this can be expressed as
U6 =51)ij e(uf)e(uf)dy+%({ F3]- e(aj.VP)e(a§/P)ay

+1 516("i)€("5ypldy

QC

(3.96)

It should be noted that the augmentation process does not add any strain energy into the

system, it merely computes the fine-scale displacements associated with the same coarse-

69

scale strain energy. It can be shown that the strain energy associated with a periodic
stifl’ness matrix at a level J is the same as the strain energy computed from a consistently
reduced stifi‘ness matrix at level J-k, i.e., the multi-resolution reduction scheme conserves
strain energy. Thus, it would be useful to express the strain energy purely as the sum of
one due to periodic displacements and another due to non-periodic displacements. This
would require the last term in the right hand side of (3.96) to be equal to zero, i.e., the

periodic component of the displacement needs to be orthogonal to the non-periodic

component with respect to the energy inner product, (0,0) E , defined as
ME = f Ej5(°)5(')dy (3.97)
QC

such that the strain energy due to the coupled (periodic and non-periodic) terms is zero.

Thus, we look to decompose the displacement 111- according to equation (3.94) such that
P P
(uj ,u}’ )E = o (3.98)
Using equation (3.94), we can write this as

<uf,(uj —uf)>E =0 (3.99)
1731- e(uf)e(uf)dy =17ij €(uj)e(uf)dy (3.100)
[u]? T D -u - = I? 5(uj)€(uf)dy (3.101)

In the discretized form, equation (3. 100) becomes

70

[uflT Kin? = [uflT Djuj (3.102)

where K]- is the periodic nodal stifiress matrix of dimensions 2112 x 2n2. The matrix
D] (of dimensions 2112 x2(n+ l)2) could be thought of as one that transforms a vector

of non-periodic nodal displacements into a periodic nodal body force. A candidate a? is
such that it satisfies the system of equations
I—( -u’-’ = D -u- (3 103)
J J .l J '

This ensures that equation (3.102) is satisfied and that the total strain energy can then be
expressed as the sum of the strain energy due to periodic displacements and that due to

the coarse, non-periodic displacements.
The periodic displacement vector u? obtained by solving equation (3.103) is one that

corresponds to the nodal values of the displacement filnction in a periodic substructure.
In order to apply the multi-resolution analysis to this filnction it needs to be expressed in

an equivalent wavelet basis, this conversion is described next.

3.7.2 Conversion to a Wavelet Basis

The conversion of a periodic filnction in the finite-element (nodal) basis Vh into an
equivalent filnction in a wavelet basis V‘” is proposed using the reverse of the process

described in section 3.4.2.3. According to this scheme a periodic nodal displacement a;
can be converted into the equivalent wavelet displacement filnction 171- E V jw as

— —l P

71

where 11,- is a reduced stiffness matrix of the substructure in a wavelet basis and B j is a

matrix that transforms a given periodic nodal displacement into an equivalent body force
in a wavelet basis. Once this coarse-scale displacement filnction in the wavelet basis has
been computed, it can then be refined using the reverse of the multi-resolution analysis
discussed in section 3.4.1 to obtain the fine-scale displacement filnction. This process is

described next.

3.7.3 Periodic Mum-Resolution Refinement

Recall the multi-resolution reduction scheme discussed in section 3 .4. 1. As expressed in
(3 .28) a single stage wavelet transform applied to a system represented by the equations

L j+1ll j+l = j+l yields the partitioned system

T

 

 

 

 

 

 

Li Cf 6,- _ ’1
r ~. T".
Cj Aj ll] f]

Thus, given the coarse-scale displacement ii]- , the detail components III of the
displacement at scale j can be obtained from
{if =Aflij—A;‘C§fij (3.105)
where A j and C J- are as defined in equation (3.29). Under the assumption that the force
is slow-varying (i.e., f = 0), the detail component is given as
(U = —A;1C§fij (3.106)
The displacement at scale j +1 can then be obtained using the inverse wavelet transform,

W—l, as follows

72

_ —1 if
uj+1=W (3.107)
fir
This process (equations (3.106) and (3.107)) is recursively repeated several times to

obtain the displacement at finest scale J .

Once the fine-scale displacement has been computed, the total fine-scale strain in the

substructure is then obtained by adding the fine-scale periodic strain, €(uJ), to the non-

periodic coarse strain, €(llj-VP ), i.e., using

a, =e(uJ)+e(u}”’) (3.108)

The fine-scale stresses can be obtained using this strain and the material distribution in

the given substructure
”J 00:15.! (y)€J (y) (3.109)

where E J represents the elastic tensor.

73

3.8 Numerical Examples

This section presents numerical examples that illustrate and compare the schemes

presented in this chapter. (Images in this dissertation are presented in color.)

3.8.1 Example 1

The first example considered is a square plate with circular perforations of various sizes
subject to a uniform compressive load at the tip as shown in figure 3.13. Here the solid
material is assumed to have a Young’s modulus of 0.91 and the weak (void) material with

a Young’s modulus of 0.045 with a Poisson’s ratio of 0.3 in both.

/

E = 0.91

E = 0.045

     

 

 

A
V
I:

 

Figure 3.13. Geometry and boundary conditions for example 1
This structure is modeled at a fine-scale that resolves all the features of the material

distribution using a commercial finite-element software. The total number of degrees of
freedom in the model is 74,498. The compliance of the structure is 0.671. The maximum

displacement is at the center of the right edge and is of magnitude 1.44. Figure 3.14

74

shows the distribution of Von Mises stress from this fine-scale model. The maximum

Von Mises stress in the structure is 2.69 near the center perforation.

it

 

 

 

 

Figure 3.14. Von Mises stress distribution using the fine-scale model of
the structure in example 1. (Number of degrees of fieedom = 74,498)

Figure 3.15 shows the assembly of 33 substructures used in constructing the reduced
model of the structure. The substructure in the center is modeled at a fine-scale
corresponding to a 64x64 pixel discretization and reduced to a scale corresponding to a

16x16 discretization. All the other substructures are also modeled at a fine-scale that

75

corresponds to a 64x64 discretization but arereduced to onethatcorresponds toan 8x8
discretization. It is noted that if all the substructures were to be modeled in the given fine
scales then the displacement across the boundary between the central substructure and
those surrounding it would not be continuous, i.e., additional constraints would need to
be applied in order to enforce the continuity. However, the reduction is done is such a
way that the displacement across substructure boundaries in the reduced models are
continuous. The total number of degrees of freedom in the reduced model is 4802, i.e., it

is approximately %6 of the fine-scale model.

Reduction:
64x64—+16xl6

Reduction:
64x64 —» 8x8

 

Figure 3.15. Assembly of substructures for example 1

76

 

 

 

 

 

 

 

 

 

 

Figure 3.17. Von Mises stress distribution in reduced model using displacement MRA

77

 

' 0.8
0.6
0.4

m

Figure 3.18. Von Mises stress in substructure (A) using the reduced model

.25

 

 

 

 

 

 

Figure 3.19. Von Mises stress in substructure (A) after augmentation

l' I.

 

    
 

Mn = 1.053-01

Figure 3.20: Detail of the Von Mises stress near substructure (A) obtained using
the fine-scale model

Error Max. om Error

2.69

 

Table 3.1: Main results for example 1

The reduced model using an MRA of material yields a compliance of 0.667 and a
maximum displacement of 1.43. The reduced model using the MRA of displacements
yields a compliance of 0.675 and a maximum displacement of 1.44. These values are

very close (less than 1% difference) to that obtained using the fine-scale model. Figures

79

3.16 and 3.17 show the distribution of the Von Mises stress using reduced models
constructed using the material MRA and the displacement MRA respectively. The
maximum Von Mises obtained from the reduced model using material MRA is 2.21 and
that from the reduced model using displacement MRA is 2.25. These values are almost
17% less than that obtained using the fine-scale model. Figure 3.18 shows a more
detailed look at the stresses in the central substructure (A). The augmentation procedure
is performed on the solution in this substructure to compute the fine-scale stresses. The
augmented stress distribution is shown in figure 3.19. The maximum stress obtained
afier the augmentation is 2.57, only 4% ofi‘ from the fine-scale. It is noted that the
coarse-scale results obtained from the two reduction procedures are in general not too
difi‘erent. However, in the case of displacement MRA, it was possible to carry out an
augmentation procedure and compute the stresses with a greater accuracy at the desired

location.

80

3.8.2 Example 2

The next example considered is the structure shown in figure 3.21. As before, the solid
material has Young’s modulus of 0.91 and the weak material (void) has Young’s

modulus of 0.045. Poisson’s ratio is 0.3 in both cases.

f=2 E=0m

E = 0.045

 

 

 

 

V

A
V7

Figure 3.21. Geometry and boundary conditions for example 2

The fine-scale model of the structure consists of 163,840 degrees of freedom and is
computed using a commercial finite element sofiware. The compliance of the structure is
6.96 and the maximum displacement is 7.66 and is observed at the center of the top edge.
Figure 3.22 shows the distribution of Von Mises stress in the structure. The maximum
stress (as expected) is observed near the sharp comers and the magnitude of the

maximum stress obtained is 4.72.

81

    
  
 
 

> 4.04a+00
< 411494-00
< 3.37e+00
< 2.7De+00
< 2,029+!!!)
< 1,353+00
< 83740-01
< 4280-05

Max = 4.729+00
Min = 4280-05

      

Figure 3.22: Von Mises stress distribution fi'om fine-scale model of the structure
with 163,840 degrees of freedom

Figure 3.23 shows the arrangement of 50 substructures used in constructing the reduced
model of the given structure. The substructures in the outer periphery of the structure
with large perforations are modeled at a fine-scale that corresponds to a 64x64 uniform
discretization and reduced to an equivalent of a 16x16 discretization. The other smaller
substructures are also modeled using the same fine-scale discretization but are reduced to
an equivalent of an 8x8 discretization. The total number of degrees of fieedom in the

reduced model is 10,240.

82

Reduction:

 

Figure 3.23. Arrangement of substructures for example 2

The reduced model using material MRA predicts a compliance of 6.7 and a maximum
displacement of 7.35. The compliance and maximum displacement obtained from the
reduced model using displacement MRA are 6.7 and 7.4 respectively. These values are
very close to those obtained fiom the fine-scale model. Figures 3.24 and 3.25 show the
Von Mises distribution obtained from reduced models build using the material MRA and
displacement MRA respectively. The maximum Von Mises stress from the material
MRA reduced model is 3.19 and that from the displacement MRA reduced model is 3.22.
These are approximately 30% less than the maximum stress obtained from the fine-scale
model. Figure 3.26 shows the coarse-scale stress in substructure (A). The fine-scale
stresses in (A) are computed by the augmentation procedure and are shown in figure 3.27.
The maximum stress after augmentation is 4.2. This is still ofi~ fiom the fine-scale result

by about 11%.

83

 

 

 

Figure 3.24. Von Mises stress distribution from a reduced model
using material MRA

 

 

   

Figure 3.25. Von Mises stress distribution from a reduced model
using displacement MRA

 

 

 

 

Figure 3.26: Von Mises stress in substructure (A) from the reduced
model using displacement MRA
F45
4

3.5

 

 

 

   

Figure 3.27: Von Mises stress in substructure (A) alter augmentation

      
 

< 4280—05

Max 2 ave-+00
Min = 4280-05

Figure 3.28: Detail of the Von Mises stress in a region around
substructure (A) from the finescale model

Max

7.66

7.35

7.4

 

Table 3.2: Main results for example 2

In this case, the maximum stress even afier augmentation has greater than 10% error as
compared to the fine-scale. This is due to the fact that sharp corners are locations of high

stress gradients and the discretization level at the coarse-scale is not enough to capture

86

these gradients. Locations of high stress gradients usually need to be modeled at a much
finer discretization than the other regions
3.8.3 Example 3

In the last example the material distribution is as shown in figure 3.29. The bottom edge

of the frame is clamped and the top edge is subject to a uniform unit load.

 

 

 

Figure 3.29. Geometry and boundary conditions for example 3

The compliance of the structure obtained from a fine-scale model using a commercial
finite element software is 0.765. The maximum displacement is 1.86. Figure 3.30 shows
the distribution of the Von Mises stresses obtained fiom a fine-scale model of the
structure. The maximum Von Mises stress is 0.83. The assembly of substructures used
in the construction of a reduced model of this structure is illustrated in figure 3.31. The
solid rim is not reduced and is modeled in accordance with the continuity requirement
across the substructure boundaries. The interior of the structure is modeled using 12
substructures modeled at a fine scale corresponding to a 64x64 mesh and reduced to a

8 x 8 mesh.

87

      
 

< 3.110-03

Max = 8309—111
Mr) == 3.110-03

 

Figure 3.30. Von Mises stress distribution from a fine-scale model

Rim Not Reduced

 

Figure 3.31. Arrangement of substructures for example 3

88

 

 

 

   

Figure 3.32: Von Mises stress distribution fi'om a reduced model
using material MRA

 

 

 

 

Figure 3.33: Von Mises stress distribution fiom a reduced model
using displacement MRA

F025

“0.2

'0.15

 

 

   

Figure 3.34: Von Mises stress in substructure (A) from the reduced

model using displacement MRA
H0. 8
' 0.7

‘0.6
‘0.5

: 30.4

0.3

 

 

 

Figure 3.35. Von Mises stress in substructure (A) alter augmentation

90

 

 

 

Figure 3.36: Detail of the Von Mises stress in substructure (A)
obtained from the fine-scale model

The reduced model using material MRA yields a compliance of 0.724 and a maximum
displacement of 1.71. The corresponding results fi'om the reduced model using
displacement MRA are: compliance of 0.756 and a maximum displacement of 1.90.
Figures 3.32 and 3.33 show the distribution of Von Mises stress obtained fiom the
material MRA model and the displacement MRA model respectively. The maximum
Von Mises stress from the material MRA reduced model is 0.61 and that from the
displacement MRA model is 0.69. These values of maximum stresses are observed at the
bottom comers. The maximum stresses from the reduced model are approximately 20-
30% less than that obtained from the fine-scale model. The augmentation procedure is

then carried out on substructure (A) and the resulting mardmum stress computed. Figures

91

 

3 .34 and 3.34 show the coarse-scale stress in substructure (A) and the augmented stress
respectively. The maximum Von Mises stress afier augmentation is 0.817, only 2% off

from the fine-scale result.

 

 

 

 

 

 

Max Disp % Error Compliance % Error Max. om % Error
Fine-Scale 1.86 - 0.77 - 0.83 -
Material 1.71 8% 0.72 < 7% 0.61 27%
MRA
Displacement 1.9 2% 0.75 < 3% 0.69 17%
MRA '
Augmented - - - - 0.82 < 2%

 

 

 

 

 

 

 

Table 3.3: Main results for example 3

It is noted that the maximum stress fi'om the reduced models are observed at the bottom
comers of the outer rim. However, in the fine-scale, the bottom corner stress though
considerably large, is not the maximum. The maximum stress in the fine—scale model is
observed in the interior substructures (A). In the reduced models, the stresses in the
substructure (A) are much less than that in the fine-scale model (compare figures (3.34)
and (3.36)). After augmentation the maximum stresses are only 2% less than the fine-
scale result. But, the augmentation over estimates the locations of the maximum stresses,
i.e., the locations of maximum stresses afier augmentation are spread out over a larger
area than as computed from the fine-scale result. Nevertheless, as seen by comparing
figures (3.34) and (3.35), there is a tremendous increase in the accuracy of the magnitude
of maximum stress computed in substructure (A) using the reduced model after

augmentation.

92

 

Chapter 4

Multi-Scale Layout Optimization Of

Structures

In this chapter, some strategies for the optimal layout design of structural systems using
materials with finite-scale heterogeneities are presented. The standard methods that are
commonly used to solve these problems are either homogenization-based techniques that
involve materials with infinitesimal heterogeneities (microstructures) (refer Bendsoe and
Kikuchi [2], Diaz and Bendsoe [12]) or methods that use fictitious material models (see
Bendsoe [3], Rozvany et al [33,34]). Structural systems involving many finite-scale
heterogeneities yield very large systems of equations that are not suitable for the kind of
iterative solution schemes that are required in the case of optimization problems. The
main goal in this chapter is to propose strategies that incorporate the model reduction
techniques discussed earlier into a problem of layout optimization of structural systems

where finite-scale features can be accounted.

93

This chapter is arranged as follows: section 4.1 gives a brief introduction to the
problem of layout optimization of structures and the standard techniques used to solve
these problems. The next section, section 4.2, introduces the problem of optimizing the
layout with finite-scale materials. There, a simplified version of this problem using
perforations as a prototype for heterogeneities is presented and the associated sensitivity
analysis is outlined. Section 4.3 discusses the formulation of compliance minimization
problems using perforated substructures. Section 4.4 illustrates the solution of
compliance minimization problems using fixed layouts of substructures. The next
section, section 4.5, deals with the dependence of the optimal layout on the arrangement
of various sizes of substructures. Here, the formulation and solution of compliance
minimization problems with varying layouts of substructures are discussed. Finally,

numerical results are presented in section 4.6.

4. 1. Background: Topology Optimization of Structures

Topology optimization problems in general seek to find the optimal layout of a structure
in a prescribed design space using a given amount of material, subject to constraints on
the response of the structure under prescribed loading conditions. Typical objective
firnctions involve: mean compliance, total mass, eigenvalues. Typical constraints include
volume and stress (e.g. see Suzuki and Kikuchi [36], Diaz and Kikuchi [13], Duysinx and
Bendsoe [l6], Haber et al [23]).

A typical topology optimization problem can be expressed as a problem that seeks

an optimal distribution of the elasticity tensor Eijkl (x) over the design domain, 0 , by

writing

94

where Em is a reference tensor and x (x) is an indicator filnction for the part 9'" of O
that is occupied by material, i.e.,

1 if xefl’"
x(x)= 14.2)
0 if x€Q\Q’"

An approach to solving such optimization problems using finite elements results
in each point, x , in the domain having the discrete choice of having material or no
material, i.e., this distributed parameter optimization problem is formulated using a
discrete valued parameter firnction. The solution of this type of problems requires the use
of discrete optimization algorithms. However, such an approach would be unstable with
respect to the choice of elements and the discretization mesh, as the distributed problem,
in general, does not have a solution unless composite materials are introduced (see Kohn
and Strang [25], Murat and Tartar [29]). An alternative solution to this problem was
introduced by Bendsoe and Kikuchi [2]. According to this scheme, rather than
determining the mixture of two materials at the macroscopic level, the mixture is allowed
to occur at an infinitesimal scale. This leads to a problem formulation using material
with a microstructure, i.e., a material with microscopic perforations of different sizes
controlled by the introduction of a parameter called effective density (volume occupied by
material in a characteristic unit-cell), which may vary continuously fi'om 0 to l, the two
limiting cases being the void and a solid material and the intermediate densities

correspond to a composite material. The relation between the effective density and the

material tensor, Erjkl (x), is determined through the use of a homogenization method,

95

where the material distribution at the microscopic level is used to detemrine the effective
properties of the material at the macroscopic level. Figure 4.1 shows a structure that is

composed of a periodic composite microstructure.

 

\

 

microstructure

characteristic

\cill

 

 

 

 

Figure 4.1: A structure with a composite microstructure

Another approach to solve the problem is by introducing an artificial density

function 11(x), x69, 0<u(x)gl, p>l and defining the elasticitytensoras

Ey'kl (x) = [#(x)lp 51ij (4-3)
This model is known as the Solid Isotropic Material with Penalization (SIMP) model (see
Rozvany et a1 [33, 34]) and yields results with fictitious materials of low stifliless for
intermediate densities for sufiiciently large value of the penalty factor, i.e., it forces the
material at each point to have an efi‘ective density that is either close to being solid or

void. Figure 4.2 illustrates a typical topology optimization problem and its solution.

96

 

 

 

e—m—s
:3

 

 

 

Figure 4.2: A typical topology optimization problem and solution

4. 2. Optimal Design Using Macro-Scale Heterogeneities

The work presented here is concerned with obtaining optimal mixtures of two materials
(solid and void) at finite scales (macro-scale), i.e., scales that may even be comparable to
the dimensions of the structure. Furthermore, it is assumed that the boundary of the
desired structure is known a-priori. This is ditferent from the case of standard topology
optimization problems where finding the boundaries is part of the problem. The problem
is now reduced to one of finding the optimal arrangement of finite-size heterogeneities on
a prescribed domain. The heterogeneities may be of any type; however, to illustrate the
idea, this discussion deals with perforations as the macro-scale heterogeneities.

A typical structure with finite-scale perforations is illustrated in Figure 4.3. The
Optimization problem seeks to find the locations and the sizes of various perforations in
the given domain. The implementation of this kind of an optimization problem using
standard techniques such as finite element methods requires the domain to be suitably
discretized and this usually involves adapting the geometry of the structure to a
conforming mesh. Since the layout of the material keeps changing during the course of

the optimization process, repeated re-meshing of the structure would be required at each

97

stage of the optimization problem, making this process not only computation-intensive

but also difficult to automate.

 

Figure 4.3: A typical structure with macro-scale heterogeneities (perforations)

One possible approach to solving such problems is proposed here, using the sub-
structuring idea described in the previous chapter. The proposed method seeks to build
structures in a given design domain 0, that can be expressed as the union of regular
substructures, DC, in such a way that all the perforations lie in the interior of the

substructures, i.e.,

n = U12, and 611, map = e (4.4)

C

where 0” represents the perforated portions of the domain. The substructures can be of
different sizes and in general each substructure can include more than one perforation.
One could construct libraries (databases) of such substructures (stifi‘ness matrices) of
various sizes and types of perforations. The optimization problem would then be

approximated into one that seeks to build a structure as an assembly of an optimal

98

selection of substructures chosen from a library of perforated substructures of various
sizes. As a result of the assumption on the design domain (4.4), the optimal structure
obtained as an assembly of selected substructures cannot have overlapping perforations
and the geometry of the perforations are limited to those in the pre-computed library of

perforated substructures. An example of such a structure is illustrated in Figure 4.4.

O

 

Figure 4.4: A structure built from an assembly of substructures with
circular perforations

4.2.1 Building a model using reduced substructures

Here, some relevant features of the construction of a model (stifliless matrix of a
structure) using reduced stiffiless matrices of substructures are presented. Consider a
structure assembled from square substructures of various sizes and with various sizes of
circular perforations as shown in Figure 4.5. The structure is represented using five

difi‘erent substructures of three difl‘erent sizes, L , 2L and 4L. Let Jc be the level of the

discretization required to resolve the material distribution in a substructure c, such that

99

 

. . . L
the discretization corresponds to a umforrn mesh wrth spacrng equal to Sc = Tc , where
2 C

LC is the dimension of the substructure.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

10L 1 4L 1 ‘2L.

0 O O 0 if?

J=7 ' "
( O O 00 1:5 c=s

Q

Q J=5

8L 0 Q 0 021 F) F,
O O V
A J 5 0:4

00 :5 4H—

00 ’ #35
c=5

 

 

 

 

 

c = 2
Figure 4.5: Assembly of substructures

Substructures of the same size (e.g., l, 2 and 3, 4) may require a different resolution
scale, depending upon the size of the perforations in them. In the figure, it is assumed

that substructures 2, 4 and 5 are discretized at the same level Jc = 5. However, the

spacing of the degrees of freedom in each of these substructures is: S2 = -‘%= g,
2
2L L L . . .
S4=——=— and S5=—=— While substructures 1 and 3 are discretized at
25 16 25 3

different levels, J]: 7 and J3 = 6, the spacing in these substructures is the same,

S1=%=% and S3 -£=%. Notice that in the fine scale problem when Sc is
2

26
difi‘erent in two adjacent substructures, a one pixel — one element finite element

discretization is in general not confirming and would require additional constraints in

order to enforce continuity across substructure boundaries. An example of such a non-

conforming mesh is shown in Figure 4.6.

Indicates non-
conforming

boundary

 

Figure 4.6: A non-conforming assembly of discretized substructures

The reduction of each substructure is performed such that the displacement across
boundaries of adjacent reduced substructures is continuous, thus avoiding additional
constraints to enforce inter-substructure compatibility. This requirement leads to the
following condition on the reduced discretization level of a substructure,

A, =Lc/zjc =A (4.5)
for c=l,2,---,Nc. For the substructures shown in Figure 4.5, this implies that
jl = j; = j5 +2 and j3 = j4 = j5 +1. It is emphasized at this point that these relations

do not determine the finest scale at which a substructure is to be modeled.

4.2.2 Constructing libraries of perforated substructures

Here, the construction of libraries (databases) of stiffness matrices and material
distributions associated with perforated substructures is discussed. A particular type of

(pararneterizable) perforation is chosen and a set of material distribution filnctions

101

(k) . . . .
{p5} associated with a discrete set of desrgn vanables in a substnrcture of srze L are

constructed at a suitably chosen fine scale J . If a reduction based on the MRA of

material distributions is to be used, then the library consists of reduced material
. . . c (k) . . . . . .
drstnbutlons {pf} at various levels 1, obtalned usrng the process described in section
3.2. If an MRA of displacement is used, the library consists of reduced stifi‘ness matrices
c (k) , , . .
{Kj} (sections 3.3, 3.4), at levels j=J,J—l,---,3 (1:3 13 the smallest level

possible when using D6 wavelets). This library is comparable to the library of effective
material tensors used in some topology optimization problems (see Suzuki and Kikuchi
[36]) and needs to be constructed only once and can be re-used in other problems with
little difi'rculty. The resolution and the number of different perforation diameters in the

library depend on the choice of the finest discretization level, J. An interpolation
scheme is then used to approximate K3- as a continuous function of the design variables

in the specified range of allowable values the design variables can take. To illustrate,
Figure 4.7 shows the largest, intermediate and smallest perforation that can be modeled at
resolutions corresponding to scales J = 4 , J = 5 and J = 6. The spacing between
nodal degrees of fieedom would be the same after reduction, regardless of the starting
level. Thus the choice of J only affects the offline computations and not the

computations within the optimization process.

102

 

a = 0.96875 a = 0.5 a = 0.03125

Figure 4.7: Substructures with circular perforations for various finest levels

4.2.3 Sensitivity Analysls

The sensitivity analysis using substructures with centered circular perforations is
illustrated here. The analysis can be extended to substructures with other kinds of
perforations easily. In the case of substructures with circular perforations, the design
variable is the diameter of the perforation in a substructure given by the following
relation

dc = aLc, 0 S a 3 am (4.6)
where LC is the dimension of the substructure and ozmax is a prescribed bound on the

size of the perforations. When the reduction of the substructures is based on an MRA of

material distribution, entries in the library are pixel values of pi, obtained using a

103

(J — j )-level reduction of the fine-scale material distribution firnction p J (a(k)). The

pixel values (1);)” are obtained fiom an interpolation of the entries in the library.

The effective stiffness matrix [(3- (a) is then computed using these interpolated

reduced material distribution functions. The gradients of the stifl‘ness matrix with respect

to the design variables is computed as

BKC- 1 6K0- 1 n
J _ J _ Z 0
6 LC 0 c k,l=1 ’

 

a C.
p] is computed using the interpolation filnction and k0 is the element

 

where p' =

stifl‘hess matrix of a solid element (reference).
When the reduction is based on an MRA of displacements, the effective stifi'ness

1k)
matrix K;- (o) is computed by interpolation in a of the entries in the library {K5} .

This interpolated firnction is similarly used in the computation of the gradients of the

stifi‘ness matrix with respect to the design variables,

 

Bdc LC Ba '
6K6-

1

where the derivative with respect to a , — , is computed using the interpolation
a

firnction.

104

4.3. Compliance Minimization Problems

While the proposed procedure can be applied using any geometry of (parameterizable)
perforations, here the formulation of the compliance mininrization problem is illustrated
for two particular cases: substructures with centered rectangular perforations and
substructures with circular perforations. It is a common practice in standard topology
optimization problems to consider characteristic cells with rectangular voids (see Suzuki
and Kikuchi [36]). This is usually done so that firlly void regions can be modeled by
letting the perforation extend to the whole cell (void cells), see Figure 4.2. In the present
case, the exterior boundary of the design domain is known a-priori and we are only
interested in the distribution of the perforations inside this boundary. The choice of the
shape of the perforation in this case may be dictated by other considerations such as case
of implementation, etc. Figure 4.8 shows two typical perforated substructures and the

associated parameters.

 

 

   

 

 

 

 

 

 

Figure 4.8: Substructures with rectangular and circular perforations

A compliance-minimization problem for a library of substructures of centered
rectangular perforations can be written as follows: For each designable substructure c,

find comma, that

105

subject to 21L: — acbc) gymeas (Q)

C
aIIIinLc S ac COS 66: b6 COS 6c 5 LC (49)
aminLc S ac sin 90, be sin BC 3 LC
0 3 BC < 1r/2
Ku = f

The design variables in this case are the dimensions of the perforations and the
orientation of the perforations with respect to the substructures.

Similarly, the optimization problem for a library of substructures of centered
circular perforations can be expressed as follows: for each designable substructure c ,

find dc that
minimize CéfTu

subject to E

C
OSchamec
Ku=f

2 d
Lc—7r—46— _<_7meas(fl) (4.10)

 

 

The design variable in this problem is the diameter dc of the perforation in each
substructure. In these problems, f is a prescribed load vector, 7 is a prescribed volume
fraction, 0 < 7 <1, and the bounds am“ and amax on the size of the perforations are
given data.

Compliance minimization problems (or any other problems) solved using the
proposed scheme can be divided into two broad categories: fixed and variable
discretization of the domain. In the problems of the first kind, the discretization of the

design domain (2 into substructures is prescribed a-priori and remains unchanged

throughout the optimization, i.e., only the size of the perforation in each substructure

106

needs to be determined. Figure 4.9 illustrates a possible starting and final layout of

perforations in an optimization problem using a fixed layout of substructures.

 

Figure 4.9: Optimization using fixed layout of substructures

In problems of the second kind, the discretization of the design domain into
subsnuctures varies at each stage of the optimization, i.e., at each step the size of each
substructure and size of the perforation are determined. Figure 4.10 illustrates the
optimization using variable layouts of substructures. However, in both types of
problems, since (consistently) reduced substructures are used the number of degrees of
freedom in the reduced model of the structure is always the same regardless of the

number of substructures used or their sizes.

 

    

 

.. 1' {-11

Initial layout

   

Figure 4.10: Optimization using variable layouts of substructures

107

 

4.4 Optimization Using Fixed Layouts 0! Reduced Substructures

In this type of problems the size LC of each substructure is prescribed a-priori, i.e., the

optimization problem seeks to find only the sizes of the perforations within a prescribed

assembly of substructures. Furthermore, from (4. 5), in order to guarantee compatibility

across substructures, LC is of the form,
Lc=2m0L (4.11)
where me is a non-negative integer and L is a prescribed dimension in the problem, see

Figure 4.6. The exponent mc may vary fi'om substructure to substructure. For
simplicity, it is assumed that the geometry in all the perforated substructures is resolved
by a discretization of the same level, J . Thus, the fine scale at any substructure is

Sc = LC /2J (4.12)

i.e., the material in each substructure is modeled using 2'] x2] (pixels) elements. In the
reduced system, the degrees of fi'eedom are spaced such that displacements are

continuous across boundaries. This determines the reduction level jc for each
substructure (see equation 4.5). If the nodal spacing after reduction is
A = L / 2”'0 (4.13)
for some positive integer m0 , then a perforated substructure is reduced from level J to
level
jc =m0+mc (4.14

Clearly, m0 must be such that for all substructures jc 3 J .

108

A possible (fixed) layout of substructures is shown in Figure 4.11. Three sizes of
substructures are used in this layout, L , 2L and 4L. In order to satisfy the continuity
requirements across substructure boundaries the reduction levels in the substructures are
relating according to the following rule: if the substructures of size 4L are reduced to a
level j, then the substructures of sizes 2L and L are reduced to levels j —l and j—2
respectively. The lowest allowable level for any substructure is j = 3 when D6 scaling
functions are used. Thus, the substructures of size 4L can be reduced at the most to a

level j = 5 in order to satisfy the previous criterion.

4 16L —__§

 

 

l

16L

 

 

 

 

 

 

 

 

 

—>l|<——>I|<—
2L L

k—H
4L

Figure 4.11: A design domain and a possible layout of reduced substructures

109

4.5 Optimization Using Variable Layouts Of Reduced

Substructures

4.5.1 Perforated Substructures and Layout Dependency

In problems with variable layouts of substructures, one needs suitable criteria to choose
between assemblies of smaller substructures and large substructures. This necessitates a
comparison between the stiffness properties of a single perforated substructure with an
assembly of smaller substructures keeping the total amount of material in both equal.
Here, an assembly of four substructures, each of size LXL, is compared with a single
large substructure of size 2L x 2L . Two tests are performed: a prescribed traction on the
boundary test and a prescribed constant pre-strain test.

The first test involves computing the resulting strain energies for various arbitrary
tractions applied on three edges and constraining all the degrees of freedom on one edge,

as shown in Figure 4.12.

11“ HF
Q J; U9: g

1/2‘

 

 

 

 

 

 

 

 

 

 

 

\ x
a 311‘ Pgan Up“ I‘gBQ

Figure 4.12: Boundary traction test
In a structure made up of several substructures subject to some prescribed loading, this
test seeks to answer the question: whether removing a large substructure with a single
perforation and replacing it with four small substructures (of half the size and with the

amount material being the same in both) would make the structure stiffer or alternatively,

110

removing a patch of four small substructures and replacing them with a single large
substructure makes the structure stifi‘er. For the same tractions, the stiffer substructure
would result in a lesser defamation and hence the resulting complementary strain energy

of a stifl‘er structure would be lesser. The random tractions are chosen arbitrarily fi'om a

uniform random distribution in lRlo’l]. The resulting strain energies in each case are
shown in Figure 4.13. The results show that a configuration of four substructures, of
size L each, is stifi‘er (lesser strain energy) that a substructure of size 2L with a single

perforation and having the same amount of material.

 

 

 

 

 

33
5
G:
l-I-l
E
l:
(I)
8
E
Q
E
O
0
who
Random Boundary Tractions

Figure 4.13: Comparison of strain energies in substructures with different

number of perforations for various boundary tractions

The second test involves the application of constant pre-strains of the form,
0 0 . . . .
e = ,nE[—l,l], (see sectron 3.x). For the same prescribed strain a strfi‘er

 

 

0

111

structure would result in higher strain energy. As before, the resulting strain energy due
to these pre-strains as obtained from a single substructure of size 2L and a patch of four
substructures of size L are shown in Figure 4.14. It can be seen fiom the figure that the
strain energy in the patch of four smaller substructures is always greater than that in the
large substructure and thus conforming that the patch of four substrucmres of size L each
is indeed stifi‘er than a single substructure of size 2L for the same amount of material in

both.

 

 

 

 

 

Strain Energy

 

 

 

 

Figure 4.14: Comparison of strain energies in substructures with difl‘erent
number of perforations for various constant pre-strains

It should be noted that these results are local to a particular substructure, i.e., they
assume that the change in the layout of a substructure does not alter the overall stress
distribution in the entire structure. However, this may not be true in general.

Nevertheless, this result may be used as a criterion in an updating scheme to determine

112

the optimal arrangement of substructures that result in the stifi‘est structure. This needs to

be done iteratively and is discussed in the following sections.

4.5.2 A Dlvldlng Approach to Optlmlzatlon with Variable Layouts

In this approach the structure is initially built entirely using as many large substructures
as possible and each update of the layout corresponds to a subdivision of a large

substructure into four smaller ones. This process is illustrated in Figure 4.15.

 

 

 

 

 

 

 

 

 

 

 

Figure 4.15: Illustration of the possible evolution of the layouts using the dividing
approach

Here, additional constraints such as a bound on the total perimeter of the perforations

(PM) or a limit on the number of substructures of a particular size are introduced in

order to have a suitable stopping criterion for the layout updating.
This process may be summarized as follows. (Here, the method is illustrated using a

bound on the perimeter of the perforations as the stopping criterion.)
(0). Start with an initial layout of substructures {LC } , c = lto N0 , with perforations of

suitable sizes such that the volume constraint is satisfied, e.g.,

dc =2Lc 1:1

7r

113

where dc is the diameter of a perforation, 7 is a prescribed volume fraction and LC

is the size of a substructure.

(l). Solve a fixed-layout optimization problem (4.10), i.e., for c =1toN0 , find optimal
dc
(2). Sort the substructures according to decreasing order of magnitude of strain energy in

each, such that
ungiui Zunguj, if i <j
where, i and j are two substructures in the sorted list whose corresponding stiffiress

matrices are Ki and Kj with displacements u,- and uj

Set N=1v0 and P=erdc

C

(3).For each substructure c in the sorted list, c =1 to No

If the perforation is of a significant size and the perimeter constraint is not violated
upon subdivision, divide the substructure into four smaller substructures with
perforations of half the size of the original perforation. Each such division increases
the total perimeter of the perforations by 7rd and the number of substructures by 3.
The requirement that the perforation be of significant size is because the division of

solid (or almost solid) substructures does not improve performance.

114

Let D=dc and Let L=Lc
( ifD>aminL
(if (P+7rD)SPmax

dc = dN+1 = dN+2 = dN+3 =

Nlhw|b

Lc = LN+1 = LN+2 = LN+3 =

N=N+3
P=P+7rD
Kendif
kendif

 

 

where amin is a suitably chosen minimum (relative) size of a perforation (e.g., 0.05 ).
Else, go to step 4 (break loop).
(4). Solve a fixed-layout optimization problem (4.10) to determine the optimal sizes of
perforations, i.e., d , 0 =1 to N (where, N > No)
The resulting layout may be used again as a starting layout and the entire process is
repeated recursively till the stopping criterion is met, i.e., any further division results in a
violation of the perimeter constraint.

A method that directly follows from the reverse idea of this approach is discussed next.

4.5.3 A Merging Approach to Optimization with Variable Layouts

In this approach, the starting layout consists of purely small substructures and each
update of the layout corresponds to merging four small substructures to create a bigger

substructure. This is illustrated in Figure 4.16.

115

 

Figure 4.16 Illustration of the possible evolution of the layouts using the
merging approach

This approach is more flexible in terms of feasible starting layouts. In the previous
approach, there is only one way to divide a substructure into four equal smaller
substructures. However, in this approach, each substructure may have up to four ways of
merging with neighboring substructures to make up a large substructure. Thus, the
number of possible layouts of substructures in this approach is a lot more than that using

the earlier approach.
This approach may be summarized as follows. (As before, the method is illustrated using

a bound on the perimeter as the stopping criterion.)
(0). Start with a layout of small substructures {LC}, c =1toNo , with perforations such

that the volume constraint is satisfied, e.g.,

dc =2Lc 1.17.

71'
where dc is the diameter of a perforation, 7 is a prescribed volume fraction and LC

is the size of a substructure.

( l). Solve a fixed-layout optimization problem (4.10), to find the optimal dc for

C=1tON0

116

(2). Group neighboring substructures (of the same size)§ as possible candidates to be
merged. A substructure can be a part of up to four groups depending on whether it is

located in an edge or in the interior of the domain. Create a table G , of dimensions

NG x 4 , that consists of all possible groups of four neighboring substructures, where
NO is the number of groups.

(3). Sort the groups according to increasing magnitude of total strain energy of the

substructures in the groups

Set P=Z7rdc and N=No

C

(4).For each group in the sorted table, g =1 to N6 , if the total perimeter exceeds the

allowable limit and if the maximum deviation of the perforation size fiom the mean is
within a prescribed bound (and the perforations are of significant size), merge the

four substructures in the group, i.e.,

 

[if P> P11m
{d}g ={d,-,dJ-,dk,d,}
(if mean({d}g)>aminLg and ::(({;};))—1 <6

 

D=fl£+fi+£+fi)

P=P—tr((d,- +dj +dk +d,)—D)
N=N—3

Kendif

\endif

 

 

 

’ relevant for recursive restarts

117

where, {i, j, k,l} is any group in G , am“ is a prescribed bound on the minimum
relative size of perforations, L8 is the size of the substructures in the group, 6 is a

prescribed bound on the maximum deviation of the perforation sizes allowable in a
group to be merged, e. g., 30% and D is the new size of the merged perforation. Each
such merging reduces the number of substructures by 3 and the reduction in the
perimeter depends on the size of the perforations merged.

(5). Solve a fixed-layout optimization problem (4.10), to find the optimal dc , for
c=1toN. Here, N<N0.

These steps are repeated recursively till the stopping criterion is met, i.e., the perimeter
satisfies the prescribed bound.

As before, the sorting of the groups and the formation of the groups themselves
must take into account the significance of the perforation sizes in the substructures to be
grouped in order to avoid unnecessarily merging solid substructures. In the previous
approach the sizes of the updated perforations are exactly half the size of the original
perforation in a substructure. In this approach, the occurrence of a group of four
neighboring perforations of the same size cannot be expected in general. It should be
noted that in this approach the starting layout has a large initial value of total perimeter of
perforations and this keeps reducing as the optimization progresses. This is in contrast to
the previous approach where the perimeter is initially less than the allowable value and
keeps increasing as the optimization progresses. Also, an important difference in the two
approaches is that the division is performed on substructures in the decreasing order of
strain energy and the merging is done on groups of substructures in the increasing order

of strain energy. In these methods, the main idea behind the changing of the layouts is

118

the knowledge that a perforated substructure of size 2L with perforation diameter 2d is
weaker than a patch of four substructures of size L each with a perforation diameter d.
Thus, at locations of high strain energy, it is better to have smaller substructures than
large ones. This also explains the use of the bound on the maximum deviation of the size
of the perforations in a group to be merged, since the criterion results from comparing a
patch of substructures with the same size perforations with a large substructure with a
perforation of double the size of the small perforations. The same criterion cannot be
used when comparing a patch of substructures with very different sizes of perforations in
each with an equivalent larger substructure with the same amount of material and so
merging of substructures with very different perforation sizes are avoided by prescribing

a bound on the maximum deviation on the perforation size in a group to be merged.

4. 6. Examples

Some numerical examples that illustrate the proposed optimization schemes for the
compliance minimization problem are presented in this section. Here the schemes are
illustrated using substructures with circular perforations. In all the examples shown, the
solid material has a Young’s modulus of 0.91 and the weak material (void) has a Young’s

modulus of 0.045 and the Poisson’s ratio is 0.3 in both.

4.6.1. Example 1

The first example considered is a simple rectangular design domain that is clamped at the
sides and two loads are applied at the centers of the top and bottom edges as shown in

Figure 4.17.

119

 

 

Figure 4.17: Problem description for example 1

The fixed layout of substructures is also shown in the above figure. The domain is
modeled using 200 substructures of equal size, as shown. The allowable volume
fraction of solid material is 0.6. The optimal layout obtained is shown in Figure 4.18.
The compliance of the structure with uniform sized perforations is 35.2. The compliance

of the optimal layout obtained is 25.2.

 

Figure 4.18: Optimal layout for example 1

120

 

4.6.2. Example 2
This example illustrates the effect of multiple load cases. The design domain is a
rectangle pinned at the bottom corners with three loads at the bottom edge applied one at
a time as shown in Figure 4.19.

24

 

V

A
‘

 

1 I 1

Figure 4.19: Problem description for example 2
The domain is modeled using 192 substructures of uniform size as shown. The volume
fraction of the solid material allowed is 0.6. The optimal layout obtained using multiple
load cases with the objective function being the mean compliance of the three load cases

is shown in Figure 4.20. The compliance of the optimal layout is 53.8.

 

 

I I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4.20: Optimal layout using multiple load cases for example 2

121

rv

Figure 4.21 shows the optimal layout obtained when all the three loads are applied

simultaneously. It can be seen that this layout is quite difi'erent from the previous layout.

 

Figure 4.21: Optimal layout using single load case for example 2

4.6.3. Example 3

The third example is an L-shaped domain with the bottom edge clamped and a unit tip
load applied on the center of the right edge. This domain is modeled using a uniform
layout of 192 substructures of equal size as shown in Figure 4.22 and using an arbitrary

layout with two sizes of substructures defined as shown in Figure 4.23.
16

 

Figure 4.22: Problem description with uniform layout for example 3

122

 

Figure 4.23: Problem description with multi-size layout for example 3

The optimal layout obtained using the uniform layout is shown in Figure 4.24 and that

obtained from the layout with multi-size substructures is shown in Figure 4.25.

 

Figure 4.24: Optimal layout using single-size substructures for example 3

123

 

Figure 4.25: Layout using multi-size substructures for example 3

The compliance of the layout with uniform size substructures is 82.8 and the compliance
of the layout with substructures of two difi’erent sizes is 86.3. Clearly, this arrangement
of substructures of two sizes is not the optimal arrangement. However, such layouts
involving larger substructures are necessary when additional constraints such as those on
the total perimeter of perforations in the structure are prescribed. The arrangement of
layouts involving multiple size substructures needs to be done iteratively and is the
subject of discussion in the following examples.

The examples presented so far require the layout of substructures to be prescribed
a-priori and that it remains unchanged throughout the optimization. The next two
examples illustrate the proposed schemes for optimization using variable layouts of

substructures.

124

4.6.4. Example 4

This example is the well-known short cantilever beam problem (see Bendsee et al [4]).
The standard topology optimization solution to this problem is shown is Figure 4.2. The
design domain is a rectangle that is clamped at the left edge and a tip load at the right
edge as shown in Figure 4.26. The volume fraction prescribed is 0.7. The layout of

substructures is not prescribed but rather it is part of the problem.

 

e———16 a

 

10

 

 

 

 

E

Figure 4.26: Problem description for examme 4

This problem is solved using the two proposed approaches for problems with variable
layouts of substructures, i.e., dividing approach and merging approach. Substructures of

two sizes, L and 2L are used in the layout optimization. A perimeter constraint
PM = 207.3 is prescribed, which is the average of the perimeters of the optimal layouts

of perforations when uniform discretization of the domain using of substructures of size L

and 2L are used. .

125

WARM

The starting layout of substructures for the dividing approach (a uniform layout of
substructures of size 2L) and the optimal layout of perforations for this arrangement of
substructures are shown in Figure 4.27. The compliance associated with this layout is

37.4. The total perimeter of the perforations is 138.2.

 

Figure 4.27: Starting arrangement of substructures for optimization using the
dividing approach for example 4 and the corresponding optimal layout of
perforations
Figure 4.28 shows the sequence of steps in the dividing approach. According to the
algorithm described in section 3.5.2, the substructures shown in figure 4.27 are sorted
according to the strain energy and are divided in sequence till the total perimeter of the

perforations reaches the prescribed limit.

126

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4.28: Sequence of steps in the dividing approach for example 4

127

The layout of substructures after successive divisions (step 18) and the optimal structure
with this as a fixed layout are shown in Figure 4.29. The compliance of this structure is

34.2 and the perimeter is 201.8.

 

Figure 4.29: Optimal arrangement of substructures using the dividing
approach for example 4 and the corresponding optimal layout of perforations

Merging Approach

The starting layout of substructures (a uniform layout of substructures of size L) and the
optimal structure using this initial layout are shown in Figure 4.30. The compliance of
this structure is 34.5 and the perimeter of the perforations is 276.5. Notice that the
compliance of this structure is already very close to that of the optimal structure using the
dividing approach. However, the perimeter of the perforations is much larger than the

prescribed limit.

128

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

II

 

Figure 4.30: Starting arrangement of substructures using the merging
approach for example 4 and the corresponding optimal layout of
perforations

Figure 4.31 shows the sequence of steps in the merging approach starting fiom the
optimal structure obtained using the initial layout. As described in section 3.5.3, the
substructures in figure 4.30 are first arranged into a list of possible candidates for
merging in an increasing order of total strain energy. Then, the groups in the list are
sequentially merged till the total perimeter of the perforations falls below the allowable

limit.

129

 

J<

at. ..

 

 

 

4
if
in.

 

.F

taunt -... - .

 

(10)

   

 

 

 

 

 

(12)

 

 

 

 

Figure 4.31: Sequence of steps in the merging approach for example 4

130

‘ I 1
. vs
Iii.
a. .. .3.
I ‘rl.
Vt .- .
I a
.
an; o
a I."
6‘. .4.
e

‘I 1. I n

 

 

 

 

 

 

 

 

The final layout of substructures obtained alter successive merging and the optimal
structure using this layout are shown in Figure 4.32. The compliance of this structure is

33.6 and the total perimeter of the perforations is 207.1.

. . ...
0.....o...

0.... - '0'.
oeoo- .0;

Figure 4.32: Optimal arrangement of substructures using the merging
approach and the corresponding optimal layout of perforations for example 4

 

The optimal structure obtained using the merging approach is stifl‘er than that obtained
using the dividing approach. However, the perimeter in the optimal layout using the
merging approach is slightly higher than the dividing approach, although both satisfy the

constraint on the perimeter.

131

 

4.7.5. Example 5

The last example is a slight variation of the classical MBB beam problem (see Duysinx
and Bendsee [16]). Here the design domain is a slender simply supported rectangle with

a uniform unit load on the top edge as shown in Figure 4.33.

illlllllllllilll ll

 

 

 

 

 

 

 

 

 

 

‘ 30 >

Figure 4.33: Problem description for example 5

As before, this problem seeks to build the optimal structure in the given design space
using perforated substructures of two sizes L and 2L. Here, to illustrate other possible
stopping criteria for the layout updating, a constraint on the number of substructures of
size 2L is prescribed as 11. This is used as the stopping criterion in the merging and the

dividing approaches to determine the optimal arrangement of substructures.

Merging Aggmach

The initial layout of substructures for the merging approach and the optimal structure
using this layout are shown in Figure 4.34. The compliance of the structure is 25.6. The

total perimeter of the perforations is 325.

132

 

 

 

 

 

III III

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4.34: Initial layout for the merging approach and optimal structure
using this layout for example 5

A sequence of layouts obtained using the merging approach is shown in Figure 4.35. The
arrangement of substructures afier successive merging and the corresponding optimal
structure using that layout is shown in Figure 4.36. The compliance of that structure is
25.4. The total perimeter of the perforations is 250. Although this perimeter is not used
at this point to determine the termination of the merging process, it is used later to
compare the optimal structures obtained using the dividing approach with two possible

termination criteria, i.e., number of substructures and total perimeter.

133

 

Figure 4.35: Sequence of layouts obtained using the merging approach for

example 5

134

 

Figure 4.36: Optimal arrangement of substructures and the corresponding
optimal structure obtained using the merging approach for example 5

Dividing Approach
The starting layout for the dividing approach and the optimal structure for this layout are
shown in Figure 4.37. The compliance of this structure is 29.9 and the total perimeter of

the perforations is 177.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4.37: Starting layout of substructures for the dividing approach and the
corresponding optimal structure for example 5

135

000- “““ 0000
......

O. r- O.
OOOoo- one...

9. '- O.

......ooo eeeeeee 00.9.09...
0... OOOO
OOOOOOOOOOO....
--tosooooo.oo
. .QOOOOO.....O

00......
0:00.00. 000000;;

 

Figure 4.38: Sequence of layouts obtained using the dividing approach for
example 5

136

A sequence of steps using the dividing approach is shown in Figure 4.38. The dividing
process continues till the number of substructures of size 2L becomes equal to the
prescribed value of 11. The final layout in Figure 4.39 is then used as the starting guess
for a fixed-layout optimization problem. The layout of substructures and the
corresponding optimal structure are shown in Figure 4.40. The compliance of this

structure is 25.3.

0- - 000
. COO.
.0 -

 

Figure 4.39: Optimal arrangement of substructures and the corresponding
optimal structure obtained using the dividing approach with prescribed number
of substructures of size 2L for example 5

It can be seen that the optimal layout of substructures and hence the optimal structure
obtained using the dividing approach when the number of substructures of size 2L are
prescribed is the same as that obtained using the merging approach.

However, if instead of the constraint on the number of substructures, a constraint on the
perimeter of the perforations (Pmax = 250) is used as the termination criterion in the

dividing approach, the dividing process terminates at step 24 shown in Figure 4.39. (The

perimeter constraint prescribed is total perimeter of the optimal structure obtained using

137

the merging approach.) The number of substructures of size 2L in this layout is 17. This
layout is then used in a fixed-layout optimization problem to determine the optimal
structure for this arrangement of substructures. The arrangement of substructures and the
corresponding optimal structure is shown in Figure 4.40. The compliance of this
structure is 25.7. This is higher than the previously obtained optimal layout using the

constraint on the number of substructures of size 2L .

0 ° ' “.OOOOOOOO . O
a e .

 

Figure 4.40: Layout of substructures and corresponding optimal structure obtained
using the dividing approach with a perimeter constraint for example 5

138

.J

 

Chapter 5

Concluding Remarks

5.1 Summary

An approach for the analysis of structural systems in linear elasticity using an assembly
of consistently reduced substructures (smaller components) was presented. Two
strategies for constructing reduced models were introduced: one based on a multi—
resolution analysis of material distributions and the other based on a multi-resolution
analysis of displacements. For relatively simple geometries, it was seen that the two
methods produce results that are not too different in the coarse-scale. However, using an
augmentation procedure it was shown that some parameters such as stresses could be
computed more accurately at certain desired locations in the case of the approach using a

multi-resolution analysis of displacements.

139

The concept of structural analysis using reduced models was then applied to the
problem of layout optimization of structures in which finite size heterogeneities of
multiple scales could be accounted. This involved constructing libraries of reduced
substructures with suitable heterogeneities (illustrated using perforations here) of various
sizes. The dependency of the optimal layout on the arrangement of substructures was
discussed. Two approaches for optimization of layouts with finite size heterogeneities
were indicated: one using fixed discretization of the design domain into substructures and

the second in which the arrangement of the substructures was part of the problem.

5.2 Conclusions

0 The concept of analysis of large structures using assemblies of reduced
substructures proved very convenient and successful in problems such as layout
optimization of structures.

0 A particular advantage is derived from the fact that the reduction of a given
substructure needs to be performed only once but the results are reusable. This
facilitates the construction of libraries of pre-computed reduced stifliress matrices
to be used as needed.

0 The use of reduced substructures of difl'erent sizes and reduction levels
accomplishes the goal of accounting for heterogeneities of different sizes and
scales.

0 As the large-scale parameters, such as compliance, are almost the same when

computed using either of the two reduction schemes the computationally

140

expensive MRA of displacements need only be used in cases where certain small-
scale parameters such as local stresses are required.

Finally, It has to be acknowledged that many serious issues need to be addressed
before this strategy can be usefirl in problems such as those arising in
crashworthiness design, such as: how to deal with dynamic behavior? — splitting
of ii into coarse-scale and details — how to accommodate nonlinear material

behavior and non-linearity associated with large deformations?

5.3 Areas of Future Work

An immediate extension of the model reduction scheme presented is to
incorporate harmonic loading. This involves the computation of a reduced mass
matrix as well as a reduced stiffness matrix. This is of interest in the design of
wave-guides and in acoustics.

Similarly, a direct extension of the presented layout optimization scheme is the
problem of optimal arrangement of beadings or corrugations on plates for
improved dynamics properties such as natural frequencies.

An approach to suitably combine the dividing and merging approaches in the
optimization using variable layouts of substructures.

Extension of the proposed schemes for three-dimensional problems in elasticity.
Incorporation of transient dynamics, non-linearity arising from material properties

as well as fiom large displacements.

141

BIBLIOGRAPHY

142

Bibliography

[1] NS. Bakhalov and AV. Knyazev, Fictitious domain methods and computations of
homogenized properties of composites with a periodic structure of essentially
different components, in G. I. Marchuk (editor), Numerical Methods and
Applications, CRC Press, Boca Raton, FL, (1994), 221-266

[2] MP. Bendsoe and N. Kikuchi, Generating optimal topologies in structural design
using a homogenization method, Computer Methods in Applied Mechanics and
Engineering, 71 (1988) 197-224.

[3] MP. Bendsoe, Optimal shape design as a material distribution problem, Structural
Optimization, I (1989), 193-202.

[4] MP. Bendsoe, A.R. Diaz and N. Kikuchi, Topology and generalized layout
optimization of elastic structures, in Topology design of structures, edited by MP.
Bendsoe and C .A Mota Scares, Kluwer Academic Publishers, (1993), 159-205.

[5] A. Bensoussan, J .L. Lions and G. Papanicolaou, Asymptotic analysis for periodic
structures, North-Holland, (1978).

[6] G. Beylkin and N. Coult, A multi-resolution strategy for reduction of elliptic PDEs
and eigenvalue problems, Applied and Computational Harmonic Analysis, 5, (1998),
129-155.

[7] M. Brewster and G. Beylkin, A multi-resolution strategy for numerical
homogenization, Applied and Computational Harmonic Analysis, 2 (1995), 327-349

[8] RR. Craig and M.C.C. Bampton, Coupling of substructures for dyncunic analyses,
AIAA Journal, 6 (7), (1968), 1313-1319.

[9] W. Dahmen and CA. Micchelli, Using the refinement equation for evaluating
integrals of wavelets, SIAM Journal of Numerical Analysis, 30, (1993), 507-537

[1011. Daubechies, Ten Lectures 0n Wavelets, CBMS-NSF Series in Applied
Mathematics, 61, SIAM, Philadelphia, 1992.

[I 1] G. C - A. DeRose, Jr. and AR Diaz, Single scale wavelet approximations in layout
optimization, Structural Optimization, 18 (1999), 1-1 I.

143

[12] AR. Diaz and MP. Bendsoe, Shme optimization of structures for multiple loading
conditions using a homogenization method, Structural Optimization, 4 (1992), 17-22

[13] AR Diaz and N. Kikuchi, Solutions to shape and topology eigenvalue optimization
problems using a homogenization method, International Journal of Numerical
Methods in Engineering, 35, (1992), 1487-1502

[14] AR. Diaz, A wavelet-Galerkin scheme for analysis of large scale problems on
simple domains, International Journal of Numerical Methods in Engineering, 44,
(1999), 1599-1616

[1 5] M. Dorabantu and B. Engquist, Wavelet-based numerical homogenization, SIAM
Journal of Numerical Analysis, 35 (1998), 540-559

[16] P. Duysinx and MP Bendsoe, Topology optimization of continuum structures with
local stress constraints, International Journal of Numerical Methods in Engineering,
43, (1998), 1453-1478

[17] B. Engquist and O. Runborg, Wavelet-based numerical homogenization with
applications, In T.J. Barth, T.F. Chan and R Haimes, editors, Multi-scale and Multi-
resolution Methods: Theory and Applications, 97-148, Vol. 20, Lecture Notes in
Computational Science and Engineering, Heidelberg 2001 , Springer-Verlag.

[18] M. Frazier, An introduction to wavelets through linear algebra, Springer-Veriag,
New York, (1999).

[19]M.I. Friswell, S.D. Garvey and J.E.T. Penny, Model reduction using dynamic and
iterated IRS techniques, Journal of Sound and Vibration, 186(2), (1995), 311-323

[20] AC. Gilbert, A comparison of multi-resolution and classical one-dimensional
homogenization schemes, Applied and Computational Harmonic Analysis. 5 (1998),
1-3 5

[21]R. Glowinski, T.W. Pan, R.O. Wells and X. Zhou, Wavelet and finite element
solutions for the Neumann problem using fictitious domains, J oumal of
Computational Physics, 126, (1996), 40-51

[22] RJ. Guyan, Reduction of stiffitess and mass matrices, AIAA Journal, 3 (2), (1965),
380.

[23]R.B. Haber, C.S. Jog and MP. Bendsae, A new approach to variable topology shape
design using a constraint on perimeter, Structural Optimization, 11, (1996), 1-12

[24] W.C. Hurty, Dynamic analysis of structural systems using component modes, AIAA
Journal, 3 (4), (1965), 678-685.

144

[25] RV. Kohn and G. Strang, Optimal design and relaxation of variational problems,
Communications in Pure and Applied Mathematics, 39, (1986), 113-137 (part I),
139-182 (part H), 353-377 (part 11])

[26] A Kunoth, Computing refinable integrals, Technical report ISC-95-02-MATH
(1995), Institute for Scientific Computation, Texas A&M University, College
Station, Texas, USA.

[27] A Latto, H.L. Resnikofl‘ and E. Tenebaum, The evaluation of connection coefiicients
of compactly supported wavelets, in Y. Maday (editor), Proceedings of F tench-USA
Workshop on Wavelets and Turbulence, Princeton University, Springer, New York,
(1992).

[28] 8G. Mallat, A theory for multi-resolution signal decomposition: fire wavelet
representation, Communications in Pure and Applied Mathematics, 41 (7), 674-693
(1988)

[29] F. Murat and L. Tartar, Optimality conditions and homogenization, in A Marino et al
(editors), Nonlinear Variational Problems, Pittman Advanced Publishing Program,
Boston, (1995), 1-8

[30] S. Pecullan, L.V. Gibiansky and S. Torquato, Scale Eflects on the elastic behavior of
periodic cmd hierarchical two-dimensional composites, Journal of the Mechanics and
Physics of Solids, 47 (1999), 1509-1542.

[31]P. Pedersen, 0n optimal orientation of orthotropic materials, Structural
Optimization, 1, (1989), 101-106

[3 2] H.L. Resnikoff and RD. Wells, Wavelet Analysis: The scalable structure of
information, Springer-Veflag, New York, (1999).

[33] G.I.N. Rozvany, M. Zhou and T. Birker, Generalized shape optimization without
homogenization, Structural Optimization, 4, (1992), 250-252

[34] G.I.N. Rozvany, M. Zhou, T. Birker and O. Sigmund, Topology optimization using
iterative continuum-We optimality criteria (C 0C) methods for discretized systems,
in Topology design of structures, edited by MP Bendsoe and CA Mota Soares,
Kluwer Academic Publishers, (1993), 273-286

[3 5] P. Seshu, Substructuring and component mode synthesis, Shock and Vibration, 4 (3),
(1997), 199-210.

[3 6] K. Suzuki and N. Kikuchi, A homogenization method for shwe and topology
optimization, Computer Methods in Applied Mechanics and Engineering, 93 (1991)
291-318

[3 7] R0. Wells and X. Zhou, Wavelet solutions for the Dirichlet problem, Numerische
Mathematik, 70, (1995), 379-396

145

[3 8] EL. Wilson, The static condensation algorithm, International Journal of Numerical
Methods in Engineering, 8 (1974), 199-203

[39]E.L. Wilson, M-W. Yuan and J .M. Dickens, Dynamic analysis by direct
superposition of Ritz vectors, Earthquake Engineering and Structural Dynamics, 10,
(1982), 813-821

[40] R-J. Yang, Metarnodeling development for vehicle fiontal impact simulation,
(DETC2001/DAC-21012), Proceedings of the 27th Annual ASME Design
Automation Conference. Pittsburgh, PA, Sept 9-12, (2001).

146

   

un]rurrrrnmrgrrrjru