'I

 

‘ .16.
(VT.-
.

.4“?

.1 . .J
9.1%..
.3:
Ian? . _ .
u. .. .4... £9. , ‘
U C'. o. t. «a i , y
a“! ”Manhunt. 4 . a
.f , ‘ ‘ ‘ . .
bluuv»
1.3.. I).
Aukfit‘lavfh 3,}:
‘ A» .. A n
_ .

‘Iutv. 1.. .-
2.!

.13....

3:.
hi

tit-1.!
.1. x. .n

w.-
lust-1t

k + .

“3%” H

gal
‘ .1 (.
,fim

any...

a‘ h 301 I
,.ll' ‘5--
. 112.4
I- 1".»
!

 

‘4]

“3.513 LIBRARY
Michigan State
University

This is to certify that the
dissertation entitled

A METHODOLOGY FOR MATERIAL DESIGN APPLIED TO
POROUS MEDIA WITH FLOW

presented by

DEEP BANDYOPADHYAY

has been accepted towards fulfillment
of the requirements for the

PhD degree in Mechanical Enfieeing

 

 

41/ We“ &W

F Major Profassor’s Signature
X /2 I /o?

Date

MSU is an afiirmative-action, equal-opportunity employer

 

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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5/08 K./Prolecc&Pres/CIRCIDateDue indd

A METHODOLOGY FOR MATERIAL DESIGN APPLIED
TO POROUS MEDIA WITH FLOW

By

Deep Bandyopadhyay

A DISSERTATION

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

DOCTOR OF PHILOSOPHY
Mechanical Engineering

2008

ABSTRACT
A METHODOLOGY FOR MATERIAL DESIGN APPLIED TO
POROUS MEDIA WITH FLOW
By

Deep Bandyopadhyay

Two methodologies to design the microstructure of porous materials
are presented in this work. The methodologies are based on shape and
topology Optimization and allow to identify layouts of the microstructure
of a periodic cells using criteria such as maximizing the effective
properties, minimizing the energy losses, or maximizing the mixing of a
dispersed solute. The porous materials studied are made of a mixture of
a heterogeneous solid matrix with its void filled with fluids. There exist
two relevant length scales in the model materials:. a microscale, which is
associated with the pores, and a large scale associated with the overall
part.

The governing equations for the microscale and the large scale are
related through the effective properties. These effective properties are
derived using the theory of homogenization. Expressions derived using
homogenization for effective properties such as permeability, dispersivity
are computed using finite element analysis and validated with
experimental results. Shape and topology Optimization are used to find

the optimal shape and layout of the microstructure.

In shape optimization, an algorithm is developed to find an optimal
shape of the pores in the microstructure for a given criterion for single and
multiphase flow. Three different shapes for the solid region of the
microstructure were analyzed: circular, elliptic and rectangular
geometries. The macroscopic fluid flow and solute transport equations
were solved based of the effective properties computed for the Optimized
microstructure. The velocity and solute fields were compared with those
computed from a microstructure form by square array of identical porosity
as of the Optimized microstructure. The result showed that the optimized
microstructure has a significant improvement in reducing energy loss
during fluid flow and increasing mixing of the dispersed solute.

Topology optimization is then used to design porous media for two
different objective functions: minimizing dissipation power and
maximizing dispersive power. The governing equations were solved using
finite element analysis and the sensitivity is computed using an adjoint
problem based of the approach. The results were evaluated by comparing
the macroscopic fluid flow for the optimal microstructure with the flow
obtained for a microstructures formed by square array with same porosity

as the optimized microstructure.

ACKNOWLEDGMENTS

I wish to express my sincere gratitude to my advisor Dr. André Bénard
and my committee members Dr. Charles Petty, Dr. Alejandro Diaz, and
Dr. Craig Somerton for their unfailing support in my PhD study. Their
advice on technical matters, encouragement, motivation and constructive
criticism were greatly acknowledged.

This work is supported by I/UCRC—MTP. This support is gratefully
acknowledged. The implementation of the method of moving asymptotes
used here was provided by Prof. Krister Svanberg from the Department of
Mathematics at KTH in Stockholm. I thank Prof. Svanberg for allowing
us to use his program.

I would also like to thank Department of Mechanical Engineering,
Michigan State University for their financial support in time Of need. I am
also grateful to College of Engineering, Michigan State University which
has provided the best research resources and facilities. I wish to thank the
entire staff of Mechanical engineering, 0188, DECS and Graduate School
for their valuable timely assistance.

Finally, I wish to express my heartfelt gratitude to my parents Biplab
and Krishna Bandyopadhyay and my in laws Madhukar and Meena
Meshram for the constant support and encouragement. Words cannot
express how thankful I am to my dearest wife Manisha for her love and

care. I would also like to thank my sister Joyeta Chatterjee, brother in

iv

laws Bipul Chatterjee and Mayur Meshram and last but not the least to my
dearest nephew Kabir. Without their support this dissertation would not

have been possible. It is to my family I dedicate this dissertation.

TABLE OF CONTENTS

 

 

 

 

LIST OF TABLES ..................................................................................... ix
LIST OF FIGURES ........................................................................ x
LIST OF FIGURES .................................................................................... x
NOMENCLATURE ................................................................................ xvi
CHAPTER 1 1
INTRODUCTION ...................................................................................... 1
1.1 General .............................................................................................. 1
1.2 Solute Transport and Fluid Flow in Porous Media ........................... 4
1.3 Research Approach ........................................................................... S
1.4 Material design techniques ............................................................... 7
1.5 Dissertation layout ............................................................................ 8
CHAPTER 2 10
THEORY FOR SINGLE PHASE FLUID FLOW IN POROUS MEDIA 10
2.1 Introduction ..................................................................................... 10
2.2 Derivation of Effective Permeability Tensor .................................. 12
2.3 Computation of Effective Permeability Tensor .............................. 19
CHAPTER 3 22
THEORY FOR SOLUTE TRANSFER IN POROUS MEDIA ................ 22
3.1 Introduction ..................................................................................... 22
3.2 Derivation of Effective Dispersion Tensor ..................................... 22
3.3 Computation of Effective Dispersion Tensor ................................. 29
CHAPTER 4 36
DESIGN OF POROUS MEDIA USING SHAPE OPTIMIZATION TO
OPT IMIZE EFFECTIVE PROPERTIES ................................................. 36
4.1 Introduction ..................................................................................... 36
4.2 Identifying the Periodic Cell ........................................................... 36
4.3 Formulation of the Optimization problem ....................................... 41
4.3.1 Maximize effective permeability ......................................... 41
4.3.1.1 Example Problems .................................................................... 43
4.3.2 Maximizing the effective dispersivity ................................. 54
4.3.2.1 Example Problems .................................................................... 55
4.4 Discussion ....................................................................................... 68
CHAPTER 5 70

 

vi

DESIGN OF POROUS MEDIA USING SHAPE OPTIMIZATION TO
MAXIMIZE “DISPERSIVE POWER” AND MIN IMIZE DISSIPATION

 

 

 

POWER ..................................................................................................... 70
5.1 Introduction ..................................................................................... 70
5.2 Optimization Of “dispersive power” and dissipation power ........... 70
5.3 Example Problems ...................................................................... 74

5.3.1 Optimization of dispersive power and dissipation power for a
given base dimension and pressure gradient ................................. 77
5.3.2 Optimization of the “dispersive power” and the dissipation
power for a given base dimension ................................................ 83
5.4 Discussion ....................................................................................... 97

CHAPTER 6 103

TOPOLOGY OPTIMIZATION OF FLUID FLOW IN POROUS MEDIA

TO MINIMIZE DISSIPATION POWER ............................................... 103
6.1 Introduction ................................................................................... 103
6.2 Formulation of the optimization problem ..................................... 106

6.2.1 Minimizing Dissipation power .............................................. 108
6.2.2 Implementation Issues ........................................................... 112
6.3 Example Problems ........................................................................ 114
6.3.1 Minimizing Dissipation Power .............................................. 116
6.4 Discussion ..................................................................................... 127

CHAPTER 7 131

TOPOLOGY OPTIMIZATION OF SOLUTE TRANSPORT IN

POROUS MEDIA TO MAXIMIZE “DISPERSIV E POWER” ............. 131
7.1 Introduction ................................................................................... 131

7.2.1 Maximizing “Dispersive power” ........................................... 136
7.2.2 Implementation Issues ........................................................... 139
7.3 Example Problems ........................................................................ 141
7.3.1 Maximizing “Dispersive power” ........................................... 143
7.4 Discussion ..................................................................................... 146

CHAPTER 8 147

THEORY FOR MULTIPHASE FLUID FLOW IN POROUS MEDIA 147
8.1 Introduction ................................................................................... 147
8.2 Derivation of Effective Permeability Tensor using Mixture Model
............................................................................................................. 148
8.3 Identifying the Periodic Cell ......................................................... 154
8.4 Formulation of the Optimization problem ..................................... 158

8.4.1 Minimize dissipation power for multiphase fluid flow ...... 158
8.5 Example Problems ........................................................................ 159
8.5.1 Minimize dissipation power for multiphase fluid flow ...... 160
8.6 Discussion ..................................................................................... 171

vii

CHAPTER 9

 

SUMMARY AND CONCLUSION .......................................................

APPENDIX A

 

APPENDIX B -

 

APPENDIX C

 

REFERENCES

 

viii

174
174

180

190

195

200

LIST OF TABLES

Table 1.1 List of engineering disciplines where porous materials are used.
In many instances the microstructure of the porous material has a direct
impact on performance. .............................................................................. 2

Table 3.1 Parameters corresponding to the geometry shown in Figure
(3.2) and which are used in the solution of the dispersion tensor ............. 32

Table 4.1 Dimensions of the microstructures shown in Figure (4.8) ........ 51
Table 4.2 Dimensions of the microstructures as shown in Figure (4.13) .63
Table 5.1 Dimensions of the periodic cells shown in Figure (5 .2) ........... 80
Table 5.2 Dimensions of the microstructures as shown in Figure (5.8). ..90
Table 5 .3 Optimized value of normalized dissipation power and
normalized dispersion power for the microstructures as shown in Figure
(5.8) ........................................................................................................... 91
Table 8.1 Dimensions of the microstructures shown in Figure (8.5) ...... 165
Table 8.2 Optimized value of normalized dissipation power and

normalized dispersion power for the microstructures as shown in Figure
(8.8) ......................................................................................................... 170

ix

LIST OF FIGURES

Images in this thesis/dissertation are presented in color.
Figure 1.1. a-d. Examples of industrial applications of porous media. ...... 3
Figure 1.2 Flow chart providing an overview of the methodology . .......... 6

Figure 1.3. a-b. Examples of (a) shape Optimization and (b) topology
Optimization from (Jorgen, 2003). ............................................................. 7

Figure 2.1 Schematic showing the particle pore diameter for porous media
from (Kaviany 1991) ................................................................................. 11

Figure 2.2. Schematic of a porous material (QM) in (a) made of an

assembly of periodic microscopic cells (52m) (b). The interface between

the solid-fluid in each periodic cell is F31 and the interface between two
periodic cells is denoted by I‘cen. .............................................................. 13

Figure 2.3 Contour plot of the distribution of the longitudinal component
of the permeability tensor (pU) for round cylinder in a square periodic

domain. The results are presented in m ........ 21

Figure 3.1. Schematic of macrostructure (Q M ) (a) consisting Of periodic
microscopic cells (9.," ) (b). The interface between the solid-fluid in each

periodic cell is 1‘81 and the interface between two periodic cells is denoted
by reel]. ......... vo ............................................................................................ 2 4

Figure 3.2 Schematic Of a periodic microscopic cells. Dimensions of a, b,
l, h are tabulated in Table (3.2) for a typical cells. “a” is the only
parameter used in the optimization of a cylinder. ..................................... 31

Figure 3.3 Comparison of dimensionless longitudinal effective dispersion
tensor measured through experiments (data taken from Buyuktas and
Wallender 2004) with numerical results obtained in this work using
homogenization for a porosity of 0.8. ....................................................... 33

Figure 3.4 Comparison Of dimensionless longitudinal effective dispersion
tensor measured through experiments for a cubic array of a porosity 0.43
(Gunn and Pryce 1969) and numerical simulations obtained with square
array of particles of porosity 0.4 (Edwards et a1. 1991) with results
obtained in this work using homogenization. ........................................... 34

Figure 4.1 Schematic Of Of the arrangement of the cells in the macroscopic
domain studied consisting of periodic microscopic cells as shown in
dashed line. ............................................................................................... 38

Figure 4.2 Schematic of periodic microscopic cells with an elliptic solid
region. ....................................................................................................... 39

Figure 4.3 Schematic of periodic microscopic cells formed from a
rectangular solid region ............................................................................. 40

Figure 4.4 Schematic of the unstructured triangular mesh used for finite
element analysis. ....................................................................................... 44

Figure 4.5 Graph of the iteration history for circular solid region with
varied step size and initial conditions. ...................................................... 46

Figure 4.6 Graph of the iteration history for elliptic solid region with
varied step size and initial conditions. ...................................................... 47

Figure 4.7 Graph of the iteration history for rectangular solid region with
varied step size and initial conditions. ...................................................... 48

Figure 4.8. a-c. Schematic of the Optimal periodic cells Obtained after
optimization. Dimensions of the periodic cells are shown in Table (4.1) 50

Figure 4.9. a-c. Schematic of distribution of longitudinal permeability
(uUn) for microstructures shown in Figure (4.5) ..................................... 53

Figure 4.10 Graph of the iteration history for circular solid region with
varied step size and initial conditions. ...................................................... 58

Figure 4.11 Graph of the iteration history for elliptic solid region with
varied step size and initial conditions. ...................................................... 59

Figure 4.12 Graph Of the iteration history for rectangular solid region with
varied step size and initial condition ......................................................... 60

Figure 4.13. a-c. Schematic of the Optimal periodic cells Obtained after
Optimization with dimensions tabulated in Table (4.2) ............................ 62

Figure 4.14 a-c. Schematic of distribution of longitudinal dispersivity
(21 1) for microstructures shown in Figure (4.7) ........................................ 66

xi

Figure 4.15 a-c. Schematic of velocity for microstructures shown in
Figure (4.7) ............................................................................................... 67

Figure 5.1 Schematic of the unstructured triangular mesh used for finite
element analysis. ....................................................................................... 76

Figure 5.2 a-h Schematic of the Optimal periodic cells obtained after
optimization for the dimensions tabulated in Table (5.1) ......................... 79

Figure 5.3 Computed values of the objective functions for a parallel flow
configuration. ............................................................................................ 8 1

Figure 5.4 Computed values of the Objective functions for a cross flow
configuration. ............................................................................................ 82

Figure 5.5 Graph of the iteration history for the circular solid region with
a varied step size and initial conditions. ................................................... 85

Figure 5.6 Graph of the iteration history for the elliptic solid region with a
varied step size and initial conditions. ...................................................... 86

Figure 5.7 Graph of the iteration history for the rectangular solid region
with a varied step size and initial conditions. ........................................... 87

Figure 5.8. a-c. Schematic of the optimal periodic cells Obtained after
optimization with the dimensions tabulated in Table (5.2) ....................... 89

Figure 5.9 a-c. Contours of the longitudinal component of the
permeability (llUr 1) for the microstructures shown in Figure (5.8) .......... 94

Figure 5.10 a-c. Contours of the longitudinal component of the
dispersivity (Z1 1) for the rrricrostructures shown in Figure (5.8) ............. 95

Figure 5.11 a-c. Plots of the velocity vectors in microstructure
corresponding to Figure (5.8) ................................................................... 96

Figure 5.12 Schematic of the macroscopic layout with applied boundary
conditions .................................................................................................. 99

Figure 5.13 Comparison of the velocity fields for a domain with a
microstructure formed with the optimized elliptic cross section with the

velocity field for a square array of same porosity ................................... 100

Figure 5.14 a-c. Comparison of solute dispersion after time t = 0.02 sec
for the macroscopic domain formed with a microstructure of elliptic cross

xii

section cylinders (a) Obtained from shape optimization with the
microstructure obtained from a square array (b) of same porosity. The
magnitude of c is shown in right (c). ...................................................... 102

Figure 6.1 Flowchart showing the steps of the algorithm based on
topology Optimization for designing microstructures. ............................ 105

Figure 6.2. Schematic of a periodic microscopic cell (9," ). The interface

between the solid—fluid in each periodic cell is F51 and the interface
between two periodic cells is denoted by FCC" ........................................ 107

Figure 6.3 Plot showing the dependence Of a on q. ................................ 110

Figure 6.4 Schematic Of a periodic microscopic cells with a liquid zones
imposed on Zone 1 above to avoid a continuous solid region between the
cells. ........................................................................................................ 1 13

Figure 6.5 Schematic of the unstructured triangular mesh used for the
finite element analysis ............................................................................. 115

Figure 6.6 Schematic of microstructure layout for volume fraction 0.5 and
with macroscopic pressure gradient of lN/m2 in horizontal direction.... 117

Figure 6.7 Schematic of rrricrostructure layout for volume fraction 0.7
with macroscopic pressure gradient of lN/m2 in horizontal direction 1 18

Figure 6.8. Schematic of the microstructure layout for a volume fraction
0.2 with an equal macroscopic pressure gradient of lN/m2 in the
horizontal and vertical directions. ........................................................... 120

Figure 6.9. Schematic of the microstructure layout Obtained for a volume
fraction of 0.3 with an equal macroscopic pressure gradient of lNlm2 in
the horizontal and the vertical directions ................................................ 121

Figure 6.10. Schematic of a microstructure layout for volume fraction 0.4
with equal macroscopic pressure gradient of 1N/m2 in the horizontal and
the vertical directions. ............................................................................. 122

Figure 6.11 Schematic of the microstructure layout for a volume fraction

0.5 with equal macroscopic pressure gradients of lN/m2 in the horizontal
and vertical directions. ............................................................................ 123

xiii

Figure 6.12. Schematic Of the microstructure layout for a volume fraction

0.6 with equal macroscopic pressure gradients of 1N/m2 in the horizontal
and vertical directions. ............................................................................ 124

Figure 6.13 Schematic of the microstructure layout for volume fraction of

0.7 with equal macroscopic” pressure gradients of lN/m2 in the horizontal
and vertical directions. ............................................................................ 125

Figure 6.14. Schematic of the microstructure layout for a volume fraction

0.8 with equal macroscopic pressure gradient of 1N/m2 in the horizontal
and vertical directions. ............................................................................ 126

Figure 6.15 Schematic of the macroscopic layout with applied boundary
conditions ................................................................................................ 129

Figure 6.16 Comparison of the velocity profiles in the macroscopic
domain for a microstructure formed form a cross section obtained from
topology optimization with the velocity field obtained for a square array
of same porosity ...................................................................................... 130

Figure 7.1 Flowchart of the topology optimization methodology
employed ................................................................................................. 133

Figure 7.2 Schematic of a periodic microscopic cell. ............................. 135

Figure 7.3 Schematic of a periodic microscopic cells with the zones to
avoid continuous solid region ................................................................. 140

Figure 7.4 Schematic of the unstructured triangular mesh used in the finite
element analysis. ..................................................................................... 142

Figure 7.5 Schematic of microstructure layout for volume fraction 0.7

with equal macroscopic pressure gradient of lN/m2 in the horizontal and

vertical directions and a concentration gradient in the horizontal direction.
................................................................................................................. 144

Figure 7.6 Schematic of the microstructure layout for a volume fraction

. . . 2 . .
0.3 With equal macroscopic pressure gradient of IN/m In the horizontal
and vertical directions and a concentration gradient in the horizontal
direction .................................................................................................. 145

Figure 8.1 Schematic of the process of representing a two phase flow
with an equivalent a single phase flow model. ....................................... 149

xiv

Figure 8.2 Schematic of a homogeneous macroscopic domain studied
consisting of periodic microscopic cells. ................................................ 155

Figure 8. 3 Schematic of periodic microscopic cells with an elliptic solid
region ...................................................................................................... 156

Figure 8.4 Schematic Of periodic microscopic cells with a rectangular
solid region .............................................................................................. 157

Figure 8.5 Schematic of the iteration history for a circular solid region
with varied step size and initial conditions. ............................................ 161

Figure 8.6 Schematic of the iteration history for an elliptic solid region
with varied step size and initial conditions. ............................................ 162

Figure 8.7 Schematic of the iteration history for a rectangular solid region
with varied step size and initial conditions. ............................................ 163

Figure 8.8 a-c. Schematic of the optimal periodic cells Obtained after
Optimization ............................................................................................ 164

Figure 8.9 a-c. Schematic of velocity distribution Of fluid mixture for
microstructures shown in Figure (8.8) .................................................... 167

Figure 8.10 a-c. Schematic of pressure distribution of fluid mixture for
microstructures shown in Figure (8.8) .................................................... 168

Figure 8.11 Schematic of the macroscopic layout with applied boundary
conditions ................................................................................................ 172

Figure 8.12. Comparison of velocity field in a macroscopic domain for
microstructure form with the cross section obtained from shape
optimization to minimize energy loss for multiphase flow with that
formed by square arrays of the same porosity ........................................ 173

XV

m

TID- Q.

azr‘NN’PF‘D‘IOD

=3

(0) (1) (2) (3)
P ,P ’P 9P

Pe

NOMENCLATURE

Major axis

Characteristic pressure tensor

Minor axis

Volumetric concentration

Periodic pressure from perturbation expansion of c
Isotropic dispersivity

Equivalent diameter of the solid region
Dispersion tensor

Molecular diffusivity

Characteristic macro height
Characteristic micro height

Identity tensor

Isotropic permeability

Permeability Tensor

Characteristic micro length
Characteristic macro length

Periodic Vector

Normal vector

Porosity

Pressure

Periodic pressure from perturbation expansion of p
Characteristic pressure

Peclet number

Tuning parameter for penalty term
Time

Periodic velocity

xvi

u(0)

<

'-<><><1>‘g

u(1)

u(

1)

)

u(3)

Periodic velocity from perturbation expansion Of u

Characteristics velocity tensor
Any periodic velocity

Small scale

Large scale

Horizontal Coordinate

Vertical Coordinate

Penalty term
Dynamic Viscosity

Macroscopic domain
Microscopic domain
Solid region
Fluid region

Length scale ratio

Boundary

Objective function

Horizontal angle between two periodic cells

Design variable

Density Of fluid

xvii

CHAPTER 1
INTRODUCTION

1.1 General

A porous medium is a mixture of a heterogeneous solid matrix with its
void filled with fluids (Kaviany, 1991). A porous media has structural
properties such as elasticity, strength and have found numerous scientific
and engineering applications such as those listed in Table (1.1).

The term “porous materials” is usually reserved to materials such as
ceramic, fibers, concrete, or naturally occurring porous rocks. But a
broader use of the term porous media can describe a wide variety of
devices or components. Figure 1.1 shows few examples of porous media
such as heat exchanger (b) (Industrial Quick Search, 2007), material for
micro chip (a) (I-l-IP-O, 2007), filters (c) (Filtration and Separation
Buyers Guide, 2007), and static mixture ((1) (Cleveland Eastern Mixers,

2007).

Table 1.1 List Of engineering disciplines where porous materials are used.
In many instances the microstructure of the porous material has a direct
impact on performance.

 

 

 

Discipline Usage of Porous Material
Agricultural engineering Dealing with drainage and irrigation
Civil engineering Concrete is a porous medium

 

Environmental engineering Groundwater pollution by toxic liquids and

hazardous wastes

 

Chemical engineering Reactors, static mixers

 

Mechanical engineering Layout Of heat exchangers can be modeled
as a porous media and micro channel

cooling

 

 

Biomedical engineering Bones, lungs and kidneys

 

 

 

 

((1)

Figure 1.1. a-d. Examples of industrial applications of porous media.

1.2 Solute Transport and Fluid Flow in Porous Media

Various approaches have been presented in the literature to derive the
governing equations for fluid flow in porous media and compute the
effective properties. An effective property is a property of a material that
changes with factors such as modification in the micro structure or the
acting boundary conditions, whereas a physical property can be measured
or perceived without changing its identity of the material such as
viscosity, molecular diffusivity.

Some noted works to compute effective properties are self—consistent
methods by Kroner (Kroner 1978), statistical modeling techniques by
Kroner (Kroner 1986), averaging methods by Quintard and Whitaker
(Quintard and Whitaker 1988) and the theory of homogenization by
Bensoussan and Sanchez-Palencia (Bensoussan 1978, Sanchez-Palencia
1980). The last method is used in this work to derive an expression for the
effective properties. Effective properties of any material are computed
from solving and averaging the resulting quantity over the microstructure.

The theory Of homogenization provides a rigorous treatment of multi-
scale problems applicable to numerous differential equations with multi
scale features. It is assumed here that there exists two length scales (micro
and macro) which are very different in magnitude. When applied to derive

the governing equations for porous media, the microstructure is also

assumed to be periodic at the micro scale and homogeneous at the macro

scale.

1.3 Research Approach

The Objective of this work is to develop methodologies for the design
of the microstructure of a porous material. To achieve this, first a review
of previous works done on the mechanics Of single and multiphase flows
and solute transport in porous media is presented. The equations derived
from the theory of homogenization (Bensoussan 1978, Sanchez-Palencia
1980) are used as described by C. C. Mei (Mei 1992) and solved using
various commercial software. The effective properties of the given
material are computed and compared with various experimental and
numerical results. The method of moving asymptotes, as proposed by K.
Svanberg (Svanberg 1978, 2002), is used to find the optimal
microstructure of the porous using shape and topology optimization.
Figure 1.2 shows a flowchart of the research approach. Step by step
algorithms were developed using shape and topology Optimization to
identify a microstructure for a porous media that Optimizes the given set
of conditions such as maximizing dispersion of the solute or minimizing

energy loss for both single and multiphase flow.

 

Define the objective funcron

Define the
design
variables

 

 

 

 

 

 

 

 

Initial Guess

 

 

 

 

!
Compute the objective function
. using finite element or finite
difference software's

 

 

 

Y

Constraints
analysis

 

 

 

 

I

 

Sensitivity analysis

 

 

 

f l
Filtering the
sensitivities

 

 

 

 

No

I 1

Optimization
with MMA

 

 

 

 

l

1
Design
variables
update

 

 

 

 

     
 

Termination condition
achieved 7

Yes

1

 

Post-processing

 

 

 

Figure 1.2 Flow chart providing an overview of the methodology.

1.4 Material design techniques

Shape and topology optimization techniques are used to determine the
shape and layout of the solid region. Figure 1.3 shows examples of shape
and topology optimization. Here shape optimization (Figure 1.3a) is used
to find the optimal shape of the holes to avoid stress concentration under
given load (Jorgen 2003, Sigmund 2000). Topology optimization (Figure
1.3b), on the other hand, is a layout optimization, which gives an optimal

layout to avoid stress concentration under given load.

Shae Ortimization

 

 

(b)

Figure 1.3. a-b. Examples of (a) shape optimization and (b) topology
Optimization from (J orgen, 2003).

The algorithm for shape Optimization tries to find an optimal shape
from a given family of shapes, which will maximize or minimize the
given Objective function. The design variables are mainly the dimensions
of the periodic cell and the solid region such as the radius for circular
cross section, major and minor axis for elliptic cross section, and length
and height for the rectangular cross section.

Topology optimization, on the other hand, mainly focuses on the
layout Of the fluid and the non-fluid region. The domain is divided into
small domains called “pixels” and for each pixel the design variable p

controls the local permeability of the medium as follows:

( ) 0if xe SolidRegion (1-1)
x =
'0 lif xe Fluid Region

For a given domain and boundary condition the algorithm tried to find the
optimal layout of the solid and fluid region, which maximizes or

minimizes the given Objective function.

1.5 Dissertation layout

The layout of this dissertation is as follows:

In chapter 2 and 3, existing theories on fluid flow and solute transport
in porous media are presented along with expressions for the effective
properties of a porous material for single phase fluid flow and solute
transport in porous media. The governing equations are solved using the

finite element method and compared with experimental results.

Chapter 4 presents my findings for the optimal microstructures Of
porous media that maximizes the effective properties such as permeability
and dispersivity. Shape optimization is used.

In Chapter 5 I try to find the Optimal microstructure for porous media
using shape optimization that maximizes dispersion while minimizing
energy loss.

Topology optimization is used to find the nricrostructures for two
objective functions: minimizing dissipation power in Chapter 6 and
maximizing “dispersive power” in Chapter 7.

In chapter 8, I review the theory for multiphase flow in porous media
and solve the governing equations using finite differences to find the
optimal microstructure that will minimize dissipation power or energy
loss for multiphase fluid flow using shape optimization. It is to be noted
however that this problem is mathematically identical to the single phase
flow problem as presented.

Finally, in Chapter 9 I summarize and discuss the results Obtained

from various Optimization techniques shown in the previous chapters

CHAPTER 2

THEORY FOR SINGLE PHASE FLUID FLOW IN POROUS
MEDIA

2.1 Introduction

The fluid flow in this work is assumed to be incompressible and non
turbulent. The theory of homogenization, introduced by Bensoussan and
Sanchez-Palencia (Bensoussan 1978, Sanchez-Palencia 1980), provides a
general framework for deducing both the macro scale equations and the
effective properties for the dynamics of rigid porous media. The main
assumption in this theory is that there exist two well—separated length
scales: the micro scale =O(x) and macro scale LzO(X),

where l/ L = 8 <<1. In this study the porous macrostructure is assumed to

be homogeneous consisting of periodic micro cells having distinct solid
and fluid regions. Figure (2.1) shows the microscopic length scale for

different naturally found porous media.

10

 

 

 

 

 

Particle Size Diameter (m)
_x
o

 

 

 

" WA

Molecular

Figure 2.1 Schematic showing the particle pore diameter for porous media
from (Kaviany 1991).

11

2.2 Derivation of Effective Permeability Tensor

The macroscopic domain (9 ) as shown in Figure 2.2 (a) is assumed to

M

be homogeneous with periodic microscopic cells (9m) as shown in

Figure 2.2(b). The periodic cell has two distinct regions: a solid region

(9s) and a fluid region (52) such that Qm=quflf

f

and as m!) = 0. The interface between the solid-fluid in each periodic

f

cell is F51 and the interface between two periodic cells is denoted by Fee“.

12

 

 

j

DVD.
i

 

 

 

 

 

 

ibbfl
bbvvvvbv
vhvvvbyy

“\bbfib'bb

 

 

 

vvvhvbvw
thrivvvv

 

 

 

 

vvvbvvvv
vavyvvvv

m
v

 

 

 

"1

(I)
.—

       

1“cell

 

 

 

(b)
Figure 2.2. Schematic of a porous material (9 M ) in (a) made of an
assembly of periodic microscopic cells (Om) (b). The interface between

the solid-fluid in each periodic cell is F31 and the interface between two

periodic cells is denoted by I‘ceu.

For a rigid porous medium with incompressible Newtonian fluid of
constant density, the governing equation for the fluid flow is Navier -

Stokes equation (Mei 1991a, 1991b, 1992).

2 (2.1)

pV u—szpfu-Vu in “f

and V-u=0 in (If (2.2)

where is u is the velocity, p is the pressure, p is the dynamic viscosity of

the fluid and pf is the density of the fluid.

The boundary condition on the solid — fluid interface is the no slip

boundary condition given by

u = Corr [“3 (2.3)

I
For a low Reynolds number flow, the left hand terms in equation 2.1 are
equally important. Since there are two different scales, the pressure term
has two contributions: the global applied pressure which has the length
scale L and the local pressure variation due to the rrricrostructure which
has the length scale 1. For generality, let there be two comparable pressure
gradients. Then the global pressure must be much greater than the local
pressure by the factor 0(l/s). These two pressure variations are referred as
the driving pressure (global) and responding pressure (local) (Mei, 1992).
The multiple scale coordinates x (small scale) and X (x=aX where X is

the large scale) are introduced to relate the micro and macro length scale

14

respectively. In order to relate the vector and scalar quantities at different
scales, asymptotic expansions for u and p are used as follows

u=u(0) +5110) +£2u(2) +£3u(3) +... (2'4)

p: pm) +1510(1) +£2p(2) ”apes +... (2.5)

where u(0),u(l),u(2) and p(0), pa), p(2)... are 9m -periodic and e is

small parameter (the ratio of the small scale over the large scale). The

. 0 1 . . .
superscripts ( ), ( ) ...denote terms assocrated With a corresponding power

. O 1 . . . .
of e, 1.6. s , s m the asymptotic expansrons. The gradient operator

related the two length scales as a function of sis given as (Mei, 1995)

V = Vx + eV (2.6)

X
where the subscripts x and X represent the small and large scales
respectively.

Substituting the gradient operator and the perturbation expansions for u

and p in the Navier-Stokes equation and the incompressibility constraint,

then

2 2.7
”(VxT‘VX)(“(0)+£u(1)+£2u(2)+£3u(3)+...)— ( )

(Vx +eVX )[p(0) +£pa) +£2p(2) + €3p(3) +...) =
p(u(0) + sum + £2u(2) + £3u(3) + ...)(Vx + EV X )

(“(0) + sum + 8211(2) +€3u(3)) in O f

15

and (Vx +eVX )-(u(0) +aufl) +£2u(2) +£3u(3) +...) = Oin (2'8)

“f

From equation 2.7 and 2.8, keeping terms Of 0(80) then the following

problem is identified

2 (0) (1) _ (0) ° (29)
pVx u —pr ..pr m (If

and V «(0) :0 in Q (2.10)
x f

(0)

In Equation (2.9) p is the global pressure applied in the macroscopic

(1)-

scale, whereas p rs the local pressure and varies in rrricroscale only.

The needed no slip boundary conditions are derived by substituting the
asymptotic expansion of u in equation (2.3) as

“(0) == 11(1) = ...... = 0 on F (solid fluid interface), (2'1 1)

(0) (1) (2.12)

11 ,u andp

(0) )

,p(1 are Q -periodic.
m
“(0), p0) and p(0) can be written in dimensionless formfim), pa)

and [3(0) related by fi(0)=u(0) l<u(0)>, [3(1) =p(1)/ P and

[3(0) = pm) /P , where P is the characteristics pressure. The length scale x

and X are sealed with characteristic macroscopic length L. Equations (2.9)

and (2.10) can be written in dimensionless form as follows,

16

2.13
fl@w» ( )

PL szfi(0) -Vxfi(1)=vx fi(0) in Qf

and Vx 43(0) :0 in Q (2.14)

f

”@m»
where __ is a dimensionless number.

PL
Equation (2.9) relates the microscopic fluid flow with the macroscopic
pressure gradient. The two terms on the left hand side of equation (2.9)
are depended on the small scale where as the right hand side term

(

V pan is a function of large scale. Assuming that u 0) and pm depends

X

on the large scale pressure gradient, “(0) and pm can also be expressed in

terms of pm) from Darcy’s equation as described by (Benssoussan et a1,

1977, and Mei 1991a)

um) =-U-VXp(O) in of (2'15)
And

. 2.
p(1)=a-VXp(O) ran ( 16)

where U is the characteristic velocity tensor and a is the characteristics

pressure tensor obtained from the solution of

17

,quZU—an :1 in Of (2.17)

and

v .U=oin§2 - (2.13)
x f

The domain is the periodic cell and the boundary conditions applied are

no slip condition at the solid liquid interface and periodic boundary

conditions for U and a on Of. U and a in dimensionless form can be
written as 0 and a related by U :U/(Lzlfl) and fi=a/L, where L is the

characteristic macroscopic length. The length scale x and X are scaled

with characteristic macroscopic length L.

v 20—v 5:1 in n (2.19)
x x m
and
.‘ .. ' 2.2
Vx 0.0 m 52m ( 0)

The effective permeability of the porous media is computed from

K =y<U> in QM (2.21)
where the averaging operator is defined as
(r) = —1— j ran (2.22)
O Q m
m f
and in dimensionless form,
K = (u) m OM (2.23)

18

Equations (2.8 - 2.9) and Equations (2.19 — 2.20) are similar yet different
in many ways. In Equations (2.8 - 2.9) the solution for “(0) depends on the

boundary conditions applied on the macroscale (macroscopic pressure

gradient) and the microstructure whereas U is independent of the
macroscopic boundary conditions. Also for a computed value of U, “(0)

can be computed from equation (2.16). It is convenient and

computationally faster to solve equations (2.8 — 2.9) when flow
parameters such as “(0) and pm are of primary focus. On the other hand if

one needs to compute the effective permeability it is convenient to solve

Equations (2.19-2.20) as shown below.

2.3 Computation of Effective Permeabllity Tensor

In this section, the commercial software Comsol® is used to solve the
characteristic velocity equations (2.20) and (2.21), for a given
microstructure of dimension 0.002m with circular solid region Of radius

0.00062m. The domain was discretized using unstructured triangular

. . -5 .
mesh wrth an average element srze of 4c and 18 shown much further

below in Figure 4.2. The physical properties of the fluid are taken as

dynamic viscosity )2 to be 0.001 Pa.s, density of fluid pf tO be 1000

Kg/m3. The effective permeability is computed from Equation (2.22).

19

Figure 2.3 shows the longitudinal permeability for a square

2

microstructure in m. The effective permeability computed for the

l E
microstructures shown in Figure 2.3 is K=3.9(10_5)[€ 1"]m2,
c

where EC = 4.07(10"3 ). When compared with the intrinsic permeability

calculated from Kozeny-Carmen equation (Bird et al. 2001) , the results
are in reasonable argument. The intrinsic permeability is given by the
Kozeny — Carmen equation is given as

"2 (2.24)

k =—
kc 233

where n is the porosity and s is the specific surface which is expressed as
the ratio of the pore surface area to the total volume of the periodic cell.

From Equation (2.24) the intrinsic permeability computed for the
. . . -12 2 . . .
microstructure mentioned above rs kkc = 1.78(10 ) m . Permeability 1s
related to the intrinsic permeability not only on the viscosity but also on

the density of the fluid at the temperature of measurement by kk — _.
C

From the above equation the permeability if computed as kkc = 1.75(10'5)

m2 which is comparable with the results Obtained using Equations 2.20
and 2.21. Since the axial velocity is directly proportional to the
longitudinal permeability, Figure 2.3 also shows the contour for axial

velocity for an applied pressure.

20

 

Figure 2.3 Contour plot of the distribution of the longitudinal component
of the permeability tensor (qu) for round cylinder in a square periodic
domain. The results are presented in m .

21

CHAPTER 3

THEORY FOR SOLUTE TRANSFER IN POROUS MEDIA

3.1 lntroductlon

In many applications of porous media such as static mixers (solute
transport), heat exchanger (heat transport), material for microchip (heat
transfer), the dispersion of a solute or heat transfer is of critical
importance. Since the governing equations for solute transport are often
similar to the equations for heat transfer, the scope Of the work applies to
both solute and heat transport in porous media.

In this chapter, the governing equation for dispersion of solute in
porous media are presented using the theory of homogenization. As
discussed in the previous chapters, the theory of homogenization, as
introduced by Bensoussan et al. and Sanchez-Palencia (Bensoussan et al.
1978, Sanchez-Palencia 1980), provides a general scheme for deducing
both the macro scale equations and the effective properties for the

dynamics of rigid porous media.

3.2 Derivation of Effectlve Disperslon Tensor

In the analysis below, the macroscopic domain (.0 M ) as shown in Figure

3.1 (a) is assumed to be homogeneous with periodic microscopic cells

22

(52m) as shown in Figure 3.1(b). The periodic cell has two distinct

regions: a solid region (9.3) and a fluid region (Q ) such that

f

Qm =98 U52 andQs an =0. The interface between the solid-

f f

fluid in each periodic cell is F81 and the interface between two periodic

cells is denoted by FCC“.

23

 

 

 

 

 

 

 

vvvbvvvw
vwvvvvvfl

 

 

 

 

 

 

 

vuvhvvvv
QWDDUDVDD
vvvhvbvv
vwvvvvvv
“\vvvbvvvv
vbvvvvvv

 

 

 

V

 

E1

Cl)
.—

     

1.“cell

 

 

 

(b)
Figure 3.1. Schematic of macrostructure ( QM ) (a) consisting of periodic
microscopic cells (9," ) (b). The interface between the solid-fluid in each
periodic cell is 1‘51 and the interface between two periodic cells is denoted
by Feell-

24

The governing equations for the solute transport in porous media at the
macroscale can be derived using theory of homogenization and the
presentation of (Mei and Auriault 1989, 1991) is followed below.

For solute transfer, the governing equation for the volumetric

concentration c of the solute in the fluid region can be expressed as

g
a:

+(u)-Vc=V-[D-Vc] in Of (3.1)

where (u) the average velocity, t is time and D is the molecular

diffusivity of the fluid.
The asymptotic expansions for c and u are given as

u =u(0) +6110) +£2u(2) +£3u(3) +... (3'2)

c=c(0) +ac(1) +€2c(2) +£3c(3) +... (3'3)

where u(0),u(1),u(2) and C(O),c(1),c(2)... are (2", -periodic and e is

small parameter (the ratio of the small scale over the large scale). The

superscripts (0), (”...denote terms associated with a corresponding power

. 0 l . . . . .
of 8, Le. a , s 1n the asymptotic expansrons. The gradient operator In

terms of the two length scales is

V=Vx+eVX (3.4)

25

where the subscripts x and X represent the small and large scales
respectively.
Substituting the differential term and expansion for c and u in volumetric

concentration equation,

3c( (0)+€C(1)+£2c(2)+£3 6(3) .) (3.5)

+

 

a:
<(u(o)£u +£u (I )+ 8211(2) +£‘311(3)+...)>-(VJr + N X )

(C(O)+£c(l)+£2c(2)+£3c(3)+. ")=(Vx +£V X')

[D-(Vx+£VX)(c (O)+.¢:c(l )+€2c(2)+£3c(3)+. NJ] in Of

Equating the terms at orders from C(80) and 0(8') for steady state

condition, then a zeroth and first order problems can be identified

0(a)); “(0V do) = W 26(0) (3.6)
x x
003'): u(0)V C(1)+u(1)V C(0)+u(0)v J0) = (3.7)
x x X
DV (v c(l)+v C(O))
x x X

Following (Mei and Auriault 1989, 1991) cm can be expressed as a
function Of the macroscopic gradient of cm) with the proportionality

function N as

(l) =_N,ch(0) in am (3.8)

N is further described as any periodic vector satisfying the following

equations

26

Vx(D(I+Vx.N»—u(0)vx .N=u(0) —<u(0)>/n in am (39)

where n is the porosity, and

II+V -N)-n=0 on 1".
x (3.10)
Nis QM -periodic and (N>=0.

From Equation (3.9) it becomes clear that N depends on the fluctuating

component of “(0) (“(0) - <u(0)>/n ). Substituting the expression for c“)

from Equation (3.8) in Equation (3.7) and reorganizing the terms, then

(u)-V c‘°’=VX-lD-ch‘°’Jin OM (3.11)

X

where (u) =<u (O)>+<u (l)> and for weak inertia flow (where the

inertial force term is negligible) (um) == 0 (Mei and Auriault 1991). The

solution for cm) is Obtained by solving the macroscopic problem.
D is the effective dispersion tensor (Mei 1992) expressed as

D: D(z) in 12 (3-12)
m

where Z is “characteristic dispersivity tensor” given as,

T (3.13)
z:(l(v N+V NT)+I)-l[(u(O)N)+(u(O)N) Jinn
2 x x 2 m

The domain for Z is the periodic cell and the boundary conditions are

that of no slip on F5, and on the outer boundary 852m of an Qm -cell, the

27

boundary condition is

n-[VXcU—1)+VXC(1))=oon rm, (3-14)

where l is an integer used to represent the relevant scales, and C“) = 0

and C0) is Q", —periodic and the vector n denotes the normal vector to

0(0), 11 and D can be written in dimensionless form 5(0),1‘r and D

related by 6(0) = 6(0) / C , fi = u /<u (0)> and I) = D/ D. Equation (3.10) can

be written in dimensionless form as

<11 (0)>L (3.15)
. -0 e .0 '
_D___<u>,vxc( )=VX.[D.VXC()]1n QM

<.. ‘°’>L
where __ is a dimensionless number.

D
It can be noted that in equation (3.11) the volumetric concentration
equation is expressed over the macroscopic domain whereas the terms in
the dispersion is defined in the microscopic domain. Hence it can be said
that the solute transport in the macroscopic domain strongly depends on

the microstructure of the porous media.

28

3.3 Computation of Effective Disperslon Tensor

In this work, the commercial software Comsol® is used to solve using
finite elements the effective dispersion tensor given by Equation (3.11) for
given boundary conditions and geometry as shown in Table (3.1) and
Figure (3.2).

In the following examples, I (h = l) is the length Of the microscopic

cell and L (H=L) is the macro-scale domain length. The fluid properties

are the viscosity [1 and density pf . A macroscopic pressure gradient is

applied along the horizontal direction. Each periodic cell is identified by
the geometric parameters such as I, h, a, b, and 6 where l is the cell
length, h is the cell height, a is the major axis or length of the solid region,
b is the minor axis or height of the solid region, and 0 is the angle between
I and h in each cell.

Examples shown in this section limits 0 to 90° to represent a square
array. Given the physical properties and the values of a and b (major and
minor axis), equation (3.11) is solved to obtain the effective dispersion
tensor. The longitudinal (along the flow direction) effective dispersion
tensor (Dxx) for two-dimensional, spatially periodic, arrays of circular
cylinders is calculated over a range of Peclet number.

The Peclet number is a dimensionless number which relates the rate of
advection Of a flow to its rate of solute or mass diffusion or thermal

diffusion. For thermal diffusion it is equivalent to the product of the

29

Reynolds number with the Prandtl number, and the product of the
Reynolds number with the Schmidt number in the case of mass diffusion.

For mass diffusion Peclet number can be computed from

(u)d p n (3-6)

where dp is the “equivalent particle” diameter, (11) is the intrinsic volume

average velocity in the direction of the pressure field,n is the porosity.
For thermal diffusion the molecular diffusivity term D is replaced by
thermal diffusivity. The longitudinal effective dispersion tensor calculated
using homogenization method is compared with experimental (data taken
from Buyuktas 2003) and numerical (Edward 1991) results and shown in

Figures (3.3) and (3.4).

30

 

 

 

 

 

 

I

Figure 3.2 Schematic of a periodic microscopic cells. Dimensions of a, b,
l, h are tabulated in Table (3.2) for a typical cells. “a” is the only
parameter used in the optimization of a cylinder.

31

Table 3.1 Parameters corresponding to the geometry shown in Figure
(3.2) and which are used in the solution of the dispersion tensor

 

 

Symbol Value Unit
1/h l [-1
UL 0.002 [-I
a/b 1 [-l

n 0.8 and 0.4 H
AP/L 1 Nm-3

77 0.001 Pa-s

L 1 m

 

32

 

"5+5T e— A -)— .. .f W- i-

 

 

I Pfannkuch ,
' I
’1? "3+4 x Edwards 81 Richardson I, .
2 A Blackwell et al. /
E . . / '
0 R11 1 I. /

gee. a “a" , '/ _
E - - homogenization I I 1
TI l’x
.5 "5+2 ~ x ,
'5 I
3 /
z , .
g 0 / ‘
a 1E+1 ~ 0 / |
d o X I, '
§ 0 ’l
D _ X I

IE+0 )§-——1—)4—-)‘<——X-""’

1E'1 I T l “r I 1 i

1E-3 1E-2 lE-l 1E+0 1E+1 1E+2 1E+3 1E+4
Pe (Peclet Number)

Figure 3.3 Comparison of dimensionless longitudinal effective dispersion
tensor measured through experiments (data taken from Buyuktas and
Wallender 2004) with numerical results Obtained in this work using
homogenization for a porosity of 0.8.

33

 

1 .E+04

 

 

- - Using Homogenization (Square Array; / o

A porosity 0.4 ’X
5 1.E+03 _ o Gunn & Pryce (Cubic Array; porosity //
E 0.43) l“ O
3. )1: Edward, Shapiro, Brenner and / o
g Shapira (Square Array; porosity 0.4) ,’ o
_ I.E+02 4 ,'§I(
E /
'c
g [/x o
g 1.E+01 ~ A o
o >(
=£ //
a O
B 1.E+00 « X , % o
a . .. — -)lGl( elem-

1.E-01 . s .

1.5-01 1.E+OO 1.E+O1 1.E+02
Pe (Peclet Number)

 

1 .E+03

Figure 3.4 Comparison of dimensionless longitudinal effective dispersion
tensor measured through experiments for a cubic array of a porosity 0.43
(Gunn and Pryce 1969) and numerical simulations obtained with square
array of particles of porosity 0.4 (Edwards et al. 1991) with results

obtained in this work using homogenization.

34

Figure (3.3) shows the relative agreement of the results obtained using
homogenization method with various experimental results. Discrepancies
with the experimental results are attributed to the fact that the experiments
were done with disordered particles; whereas, this work is done using an
array of particles equally spaced with periodic boundary condition. Also
in homogenization, it is assumed that the Peclet number is of the order of
1 (low Reynolds number flow); and, hence, some significant variations
with experimental results for higher Peclet number are expected.

This can be further compared with the numerical analysis based on
Taylor dispersion theory as shown in Figure (3.4). In Figure (3.4), some
discrepancies are Observed with the experimental data (Gunn 1969). The
disagreement between the present result and those of the experiments is
due to the fact that their data from the experiments are for a three-
dimensional cubic array; whereas, the present work is for two-
dimensional, spatially periodic arrays of circular cylinders.

From the above comparison it appears that homogenization provides
very reasonable values for the dispersion tensor for a Peclet number Of

order one.

35

CHAPTER 4

DESIGN OF POROUS MEDIA USING SHAPE OPTIMIZATION
TO OPTIMIZE EFFECTIVE PROPERTIES

4.1 Introduction

Shape optimization is a well established method for design. Here the
problem is applied to material design and posed in such a manner that the
algorithm may find an optimal shape of the pores in the microstructure
that will optimize the effective properties. The optimization problem is
cast in the fashion of a standard shape optimization problem where
various shapes described parametrically are studied. The Method of
Moving Asymptotes as proposed by Krister Svanberg is used for
Optimization. The derivatives of the objective function with respect to the
design variables were computed using finite difference. The governing
equations needed to solve the effective properties were derived in

Chapters three and four and are solved using the finite element method.
4.2 Identifying the Periodic Cell

In the analysis, the macroscopic domain (Q M ) is assumed to be
homogeneous with periodic microscopic cells ( 52m ). The microscopic
cell has two distinct regions: a solid region ((23) and a fluid region (9

f)

36

such that Qm=quQ andQSnQ =0. P, l and L are the

f f

characteristic pressure, micro length and macro length scales. Two scales

are introduced: a small scale (xz 0(1)) and a large scale (X z0(L)).
Both scales are related by X = x/e.

The homogeneous macroscopic domain shown in Figure (4.1) is
formed of periodic microscopic cells as shown in Figure (4.2) and Figure
(4.3). Each cell geometry is identified by the following parameters: 1, h,
a, b, and 6. L is the cell length, h is the cell height, a is the major axis or
length Of the solid region, b is the minor axis or height of the solid region,
and 6 is the angle between I and h in each cell. Three different shapes for
the solid region such as circle (a = b), ellipse (Figure (4.2)), and rectangle

(Figure (4.3)) are considered below.

37

 

 

 

 

Y
Lx.

 

 

 

L

Figure 4.1 Schematic of of the arrangement of the cells in the macroscopic
domain studied consisting of periodic microscopic cells as shown in
dashed line.

38

 

Figure 4.2 Schematic of periodic microscopic cells with an elliptic solid
region.

39

 

I
l. ................... .1 1

Figure 4.3 Schematic of periodic microscopic cells formed from a
rectangular solid region

40

4.3 Formulation of the optimization problem

The aim is to find a microstructure of the periodic cell that maximizes
some function of the effective properties with the added constraint of
volume fraction using shape optimization. In the following section, the

Objectives functions are presented.

4.3.1 Maximize effective permeability

Filters constitute of the common applications of porous. It is common in
such an application to be concerned with the pressure drop across the filter
media. Hence the motivation behind the work in this section is to design a
porous media that maximize the permeability for a given porosity or
volume fraction of solid i.e. to seek a periodic cell that maximizes the
permeability with isotropic flow symmetry as proposed by Guest (Guest
2007) for a given porosity or volume fraction of solid. The rrricrostructure
is modified such that the fluid in a porous media only flows in the
direction of the applied pressure gradient, and hence follows a reduction
in the loss of energy in fluid flow.
The isotropic effective permeability is defined as

kI=K in (2 (4.1)
m

The complete formulation of the optimization problem is stated as
follows:

Find the Optimal value of6, l, a and b for a g fixed of h, that will

41

maximize: Isotropic Effective Permeability

where isotropic permeability is expressed as k - 8

Here k and 8p are computed from

i=1. i K- ‘4'”
21' =1 ”
and
. (1)1 2 (K -K )2 . g i (K )2] ‘4'”
p k i=1 ii (i+1)(i+1) i=1j=2 ij

where K1] is computed from equations (2.18).

subjected to: volumetric and geometric constraints

Constraint 1: 53—S- 2 V

9" fiac (4.4)
Constraint 2: a (é
hsin 0

Constraint 3: b <

 

where Os is the volume of the solid region, Om is the total volume of the

cell, and mec is value (0<mec<l) specifying the minimum solid volume.

The first constraint in equation (4.4) restricts the minimum volume of
solid region to mec. The next two constraints are imposed to insure the

continuity of each phase.

42

4.3.1.1 Example Problems

This section shows examples of shape optimization that maximizes
effective permeability as described in section 4.3.1. The solutions were
Obtained by casting the problem in the fashion of a standard shape
optimization problem. The derivatives of the objective function with
respect to the design variables were computed using a forward finite

difference method as shown below,

3013i __ Objdv + &lv — Objdv (4'5)
adv - (ilv

 

 

where “Obj” is the objective function and “dv” is the design variable. The
derivatives of the constraints with respect to the design variables were
computed analytically. The objective function was computed by solving
the governing equations for the effective properties using the finite

element method. The domain was discretized using approximately 6000

. . . . -5
triangular elements in an unstructured mesh wrth an element srze of 4c

m as shown in Figure (4.4).

43

. 21:31.5, . 11-51.
v A. \r‘t. 1 s
. 13].. 3, ti
.wu/y. . tum» ......i».

- T '1
meanwhnfi

. 14...).
15.54:»... _.
Erase...

 

element analysis.

Figure 4.4 Schematic of the unstructured triangular mesh used for finite

The Method of Moving Asymptotes (MMA) is used for optimization.
The results shown in this section are all for a limiting volume fraction of

0.3 and fixed value of 0.002 m for h. The physical properties of the fluid

are taken as dynamic viscosity it to be 0.001 Pas, density of fluid pf to be

1000 Kg/m3. The dependence of the step size and initial conditions in the

Optimization algorithm were verified by varying it for two different cases.
For case 1 the step size is set to 0.1 and an initial aspect ratio is 1, radius
(major axis, minor axis) set to 5x104m and 0 value to 90°. For case 2 the
step size is set to 0.5 and an initial aspect ratio is 1 (fixed height of
0.002m), radius (or major axis, minor axis) set to 7x104m, and 6 value to
90°. The iteration history for the different cases are shown in Figure (4.5 —

4.7) for circular, elliptic and rectangular solid regions.

45

lterration History

 

 

25~ """ l”

N
o
1

 

_A.
a.
1

—CASE 1 - - -CASE 2

 

 

 

..5
o
J

 

Objective Function

 

 

 

 

o 20 40 60 so
Iteration Number

Figure 4.5 Graph of the iteration history for circular solid region with
varied step size and initial conditions.

46

100

Iteration hlstory

 

70

60-

   

 

    

   

' - -CASE1 —-CASE 2

Objective Function

 

 

 

 

L

I

0 20 40 so so 100 120
Iteration Number

 

Figure 4.6 Graph of the iteration history for elliptic solid region with
varied step size and initial conditions.

47

Iteration History

 

 

 

 

 

 

 

 

 

 

35
304
C 5 :0.”
2 25" 0:"
ii ' =
I2 201 '
0
.5 15 . .
-‘-’- ; —CASE1 ---CASE2
3 10 'I I
° :
5 - .' J
o 3* . . . . .
o 20 40 60 so 100 120

Iteration Number

Figure 4.7 Graph of the iteration history for rectangular solid region with
varied step size and initial conditions.

48

The optimal microstructure of the porous medium for the desired
objective function to maximize isotropic effective permeability as
discussed in section 4.3.1 for volume fraction of 0.3 for solid region is
shown in Figure (4.8). The dimensions of the microstructures are shown

in Table (4.1).

49

 

 

 

 

 

 

 

 

 

(b)

 

 

 

.-____-___-______-_---J

 

(C)

Figure 4.8. a—c. Schematic of the optimal periodic cells obtained after
Optimization. Dimensions of the periodic cells are shown in Table (4.1)

50

Table 4.1 Dimensions of the microstructures shown in Figure (4.8)
Figure (4.8) l/h a/h b/h l9

 

 

 

 

 

(a) . l 0.31 89.6
(b) 0.65 0.19 0.33 89.6
(c) 1.07 0.25 0.28 89.6

 

 

 

 

 

 

 

51

The effective permeability computed for the microstructures shown in
Figure (4.8) a, b and c are;

1 e 1 e
K = 3.9(10‘5)[£ C]m2, K = 2.5(10‘5)[‘E e]m2 and
C

r

1 e
K=3.5(10‘5)[£ ']m2

where ac = 407(10‘3), 2e = 222(10‘3), and er =1.35(10‘ 3)

respectively. Figure (4.9) shows a contour of the distribution of

longitudinal permeability (uUjl). It is observed in Figure (4.8) that the

nricrostructure forms a square array for both the cases with circular and
rectangular microstructure. Among the three microstructures, the circular
microstructure gives the maximum effective permeability; but the
diagonal terms of the permeability tensor are considerably higher than the
other microstructures. Overall, the microstructure with rectangular solid
region shows the best result. This may be due to the fact that it can be
considered as a flow in a channel. It can be noted that for a given pressure
gradient the velocity will have similar profile as the longitudinal

permeability.

52

x105

I14
a l7
o

(a)

x105

. 10
i 5
I o

(b)

x105

I12
u .6
'0

(C)

Figure 4.9. a-c. Schematic of distribution of longitudinal permeability
(”U1 1) for microstructures shown in Figure (4.5)

53

4.3.2 Maximizing the effective dispersivity

Dispersion of a solute can be of importance for applications involving
reactors. The governing equations for solute transport are similar to the
equation for heat transfer.. Industrial examples of solute or heat transport
in porous media include static mixers (solute transport), heat exchanger
(heat transport), and material for microchip (heat transfer). For numerous
applications it is important how to disperse the solute or heat as fast as
possible. In this section, the optimization algorithm seeks a periodic cell
that maximizes the dispersion tensor. Maximizing the dispersion tensor
maximizes the solute or heat transfer. The objective function is posed in a
manner such that the algorithm also finds the required pressure gradient to
maximize the isotropic dispersivity, which can be defined as

d! = D (4.6)
The complete formulation of the optimization problem is stated as

follows:

Find the optimal value of6, l, a, b V (O) and V (0) for a

X” 1'”
given value of h, that will
maximize: Isotropic Effective Dispersivity

where isotropic effective dispersivity is expressed as d — 8d

Here d and 8d are computed from

54

1 2 (4.7)
d=— z 0..
2i=1 U
and
1 2 I 2 2 1 2 (4-8)
€d=(2) i§1(Dii‘D(i+1)(i+1)) +i§le2(Dij)

Dij is computed from solving equation (3.12).
subjected to: volumetric and geometric constraints

(2 (4.9)
Constraint 1: —i 2 V

m frac

Constraint 2: a < .1—
2

Constraint 3: b < M

where Os is the volume of the solid region, Om is the total volume of the

cell, and mec is value (0<Vf,ac<l) specifying the minimum solid volume.

The first constraint in equation (4.9) restricts the minimum volume of

solid region to mec. The next two constraints are imposed to insure the

continuity of each phase.

4.3.2.1 Example Problems

This section shows examples of shape optimization that will maximize the
dispersivity as descried in section 4.3.2. The solutions were obtained by

casting the problem in the fashion of a standard shape optimization

55

problem . The derivatives of the objective function with respect to the
design variables were computed using forward finite difference method as

shown below,

301))“ = Obj dv +a1v —Objdv (4.10)

 

adv (ilv

where “Obj” is the objective function and “dv” is the design variable. The
derivatives of the constraints with respect to the design variables were
computed analytically. The objective function was computed by solving
the governing equations for the effective properties using the finite

element method. The domain was discretized using an unstructured

. . . -5
triangular mesh wrth 6000elements of an average srze of 4c in as shown

in Figure (4.4) Thee Method of Moving Asymptotes (MMA) as proposed

by Krister Svanberg is used for optimization.

The results shown in this section are all for a limiting volume fraction of

0.3 and fixed value of 0.002 m for h. The physical properties of the fluid

are taken as dynamic viscosity [1 to be 0.001 Pa.s, density of fluid pf to be

1000 Kg/m3, molecular diffusivity D set to 10.9 m2/s. The dependence of

the step size and initial vale were in the optimization algorithm were
verified by varying it those for two different cases. For case 1 the step size
is set to 0.1 and initial aspect ratio is 1 (fixed height 0.002m), radius (or

major axis, minor axis) set to 5x104m and 0 value to 90°. For case 2 the

56

step size is set to 0.5 and initial aspect ratio is l, radius (or major axis,
minor axis) set to 7x104m and 0 value to 90°. The iteration history for the
different cases are shown in Figure (4.10 - 4.12) for circular, elliptic and

rectangular solid region.

57

Iteration History

 

.......

 

 

 

300 . - - -CASE1 —CASE 2

 

 

Objective Function
N
O
O

 

 

 

150 .
100 2
50 i
0 A i T . 1
0 20 40 60 80 100

Iteration Number

Figure 4.10 Graph of the iteration history for circular solid region with
varied step size and initial conditions.

58

Iteration History

 

450
400 i

8
8 o
L .l

M M (A)
O 01

O O

1

Objective Function
e

100 4

 

  

 

- - - CASE 1 -—CASE 2

   

 

 

.....
'O

 

...
_"

 

Iteration Number

Figure 4.11 Graph of the iteration history for elliptic solid region with

varied step size and initial conditions.

59

Iteration History

 

 

 

250

200 r

—L
(1"
o

 

d
o
o

 

 

 

Objective Function

50*

 

 

 

 

0 20 40 60 80 100
Iteration Number

Figure 4.12 Graph of the iteration history for rectangular solid region with
varied step size and initial condition.

60

The optimal microstructures of the porous medium for the desired
objective function to maximize isotropic effective dispersive as discussed
in section 4.3.2 for volume fraction of 0.3 are shown in Figure (4.13). The

dimensions of the microstructures are shown in Table (4.2).

61

"OQOOiOOOOI.

 

U

I

loo-.0000...
(a)

'I...‘.......'

 

0
l
O
l
I
I
O
I
O
I
I
J

(b)

-------------q

1

 

'0.-.-._.-..
‘-.-.D-.---

......O......
(C)

Figure 4.13. a-c. Schematic of the optimal periodic cells obtained after
optimization with dimensions tabulated in Table (4.2)

62

Table 4.2 Dimensions of the rrricrostructures as shown in Figure (4.13)

 

Figure (4.13) l/h a/h b/h a

 

 

 

 

(a) 0.94 0.30 79.3
(b) 1.07 0.33 0.31 73.4
(c) 1.18 0.33 0.27 87.3

 

 

 

 

 

 

 

63

0
. . V X p ( ) —1 3
The optimal pressure gradients computed are (0) = 0 9 N/m ,
V p T '
Y

(0) V (0)
- -
v p(0) —0.5 I "I H" -0.4

] N/m3 for Figure
Y

(4.13) a, b, and c respectively. The effective dispersion tensor computed

for the periodic cell shown in Figure (4.13) a, b, and c are

D =1.5(10 ) C m /s, D = 1.6700 ) e m /s and
ac 1 are 1

-3 1 a 2
D =1.33(10 ) r m /s
ar 1

where ac = 2000‘ 2), are =1.80(10" 2) and ar = 602(10‘ 2)

respectively.

One notable observation made from Figure (4.13): the angle between
the length and the height goes to a lower value for circular and elliptic
cross section. But for a rectangular cross section, the 0 value approaches
90. This complements the previous argument that it forms a channel flow;
and hence, it requires a comparative smaller pressure gradient to optimize
the dispersion tensor, which is of the same magnitude of the other micro

structures.

Figure (4.14) shows a contour of the distribution of the longitudinal
component of the dispersion tensor (Zn) in m2/s. It can be noted from
Figure (4.14) that the longitudinal term in the dispersion tensor is high in
the horizontal direction, hence higher solute or heat transfer in the given
direction. Figure (4.15) shows the velocity vectors which are periodic in

nature.

65

 

-10

(b)

x1043

-1
(C)

Figure 4.14 a-c. Schematic of distribution of longitudinal dispersivity
(21 1) for microstructures shown in Figure (4.7)

(a)

(C)

Figure 4.15 a-c. Schematic of velocity for microstructures shown in
Figure (4.7)

67

4.4 Discussion

In this chapter a methodology to find the optimal microstructure for
two different conditions; maximize effective permeability and effective
dispersivity using shape optimization.

Attempts are made to maximize the effective properties by modifying
the shape of the periodic cell. The magnitude of effective permeability
computed is comparable to the work done by Guest and Prévost (Guest
2007). It is observed in Figure (4.8) that the microstructure forms a square
array for both the cases with circular and rectangular microstructure.
Among the three microstructures, the circular microstructure gives the
maximum effective permeability; but the diagonal terms of the
permeability tensor are considerably higher than the other microstructures.
Overall, the rrricrostructure with rectangular solid region shows the best
result. This may be due to the fact that it can be considered as a flow in a
channeL

In the next optimization example, the effective dispersion tensor is
maximized. Here we try to maximize the effective property not only by
modifying the shape of the periodic cell but also by finding the
appropriate boundary condition. One notable observation made from
Figure (4.13) is that the angle between the length and the height goes to a
lower value for circular and elliptic cross section. But for a rectangular

cross section, the 0 value approaches 90. This complements the previous

68

argument that it forms a channel flow; and hence, it requires a
comparative smaller pressure gradient to optimize the dispersion tensor,
which is of the same magnitude of the other micro structures.

Finally, it can be concluded from the examples that a small change in

geometry changes the effective properties significantly.

69

CHAPTER 5

DESIGN OF POROUS MEDIA USING SHAPE OPTIMIZATION
TO MAXIMIZE “DISPERSIVE POWER” AND MINIMIZE
DISSIPATION POWER

5.1 Introduction

A porous medium may have to satisfy more than one desired criteria
and such problem is studied in this chapter. The desired criteria studied
are maximizing the “dispersive power” and minimizing the dissipation
power. The possible applications for such a medium include a static mixer
(fluid flow and solute transfer), a heat exchanger (fluid flow - heat
transfer), where the objective is to mix or disperse the solute (or heat)
with minimum energy loss.

The problem is cast in the fashion of a standard shape optimization
problem where various shapes described parametrically are studied. The

periodic cells are defined in similar manner as described in Section 4.1.

5.2 Optimization of “dispersive power” and dissipation power

The topic of interest of this section is to find a microscopic geometry
that will minimize the dissipation power while maximizing “dispersive

power”.

70

Two different approaches were taken to find the optimal
microstructure that minimize dissipation power and maximize dispersive
power. In the first approach the pressure gradient was fixed and in the
second approach the pressure gradient was kept as a design variable. The
motivation behind two separate approaches was to understand the effect of
the pressure gradient on microstructure in dispersion of fluid.

The complete formulation of the optimization problem for first approach
is stated as

Problem definition 1:

Find6, l, a and b for a given value of h, v X pm) and V pm) that

Y

will

minimize: 3dissipation — 3dispersive

The complete formulation of the optimization problem for second
approach is stated as

Problem definition 2:

Find 6, l, a, b V X pm) and VYp(O) for a given value of h, that will

mimmrze: SSdissipation — Sdispersive

71

 

Optimization algorithm MMA always minimize the objective function

and since we are interested in minimizing 3 d

dissipation an

“'8. .,wi ntivsi frisriv
maxrnrrzrng disperszve e ncludea ega e gnbeoed pes e

power term.

subjected to: volumetric and geometric constraints

Q
Constraint 1: —s 2 V
9 frac

"I

Constraint 2: a (é (5.1)

hsin6

 

Constraint 3: b <

where Os is the volume Of the solid region, Om is the total volume of the

cell, and mec is value (0<Vf,ac<l) specifying the minimum solid volume.
The first constraint in equation (5.1) restricts the minimum volume of
solid region to mec. Constraints 2 and 3 are imposed to insure the

continuity of each phase.

disipation is the normalized value of dissipation power (Borrvall

2003). An expression for the dissipation energy for the fluid in the unit

cell is derived from a mechanical energy balance. This is obtained by

multiplying equation (2.9) with 1"”), where Va” is any weighted periodic

72

velocity satisfying equation (2.10) and boundary condition as shown in

equation (2.12), and integrating over the entire domain

1: x -
Q n
m m
l Vx pm) .40)“,
Q
m
where the double dot operator (z) is defined here as a ; p = aij'Bji' The

first term in equation (5.2) is the viscous power and is responsible for

dissipating the energy during fluid flow. This term is also called the

dissipation power. After setting vm) = “(0)

then, (5.3)

=fl j qu(0):(qu(0))TdQ

Q
m

dissipation

The dissipation power is sealed with the square of the trace of the
effective permeability after the expression for the velocity field is
substituted

The dissipation power can be written in a dimensionless form as follows

—_— j thiw) :(Vxfi(0))TdQ (5'4)

9

m

dissipation

. . is the normalized value of “dispersive power”, which is
disperszve

computed by multiplying equation (3.10) with 0(0), and integrating over

73

the entire domain, the following expression for “dispersive power” is

obtained after using Green’s theorem:

T 5.5
(u) j vchWr-(OMO: [VXc(O)-D-[VXC(O)] an ( )

Q Q
m m

The right hand term is dispersive power term.

= l chw) -D-[ch(0)]ng, (5.6)

Q.
m

dispersive

where D is the effective dispersion tensor, which can be computed from

D 2 D<(—;-(va +vaTj+Ij gfiuwinjthmmjr»

The dispersive power term is scaled with the product of the trace of the
dispersion and the trace of the concentration gradient. Dispersive power

can be written in dimensionless terms as follows

T (5-7)
5(0)] do,

" _ ~(0) , " ,
SSdispersive _ Q! VX c D [V
m

X

5.3 Example Problems

This section shows few examples of shape optimization for the objective
functions descried in section 5.2.1 and 5.2.2. The solutions were obtained
by casting the problem in the fashion of a standard shape optimization

problem as before. The derivatives of the objective function with respect

74

to the design variables were computed using forward finite difference

 

method as
30W _ Objdv + adv — Obj dv (5:8)
adv - div 1

where “Obj” is the objective function and “dv” is the design variable. The
derivatives of the constraints with respect to the design variables were
computed analytically.

The objective function was computed by solving the governing equations
for the effective properties using finite element method. The domain was

. . . . . -5
discretized usrng unstructured triangular mesh of element srze of 4e m

as shown in Figure (4.1).

75

Avail-g -

é?

a
u

..,. }
31

V ’
' 1 UV

,. . 3‘

ix; 31$?“th

"fl
is

-‘

,4

ctr} t" 3‘
1...,M' .‘

“um...
10’“

W.

21-41:

:3“;

 

.15

Figure 5.1 Schematic of the unstructured triangular mesh used for finite
element analysis.

76

The Method of Moving Asymptotes (MMA) as proposed by Krister
Svanberg is used for optimization. The results shown in this section are all
for a limiting volume fraction of 0.3 and fixed value of 0.002m for h. The

physical properties of the fluid are taken as dynamic viscosity )1: to be

0.001 Pas , density of fluid pf to be 1000 Kg/m3, molecular diffusivity D

set to 10.9 mZ/s.

5.3.1 Optimization of disperslve power and dissipation power

for a given base dimension and pressure gradient

As mentioned earlier, the motivation of this work is not only to find
the optimal rrricrostructure for given boundary conditions but also to
understand the effect of pressure gradient on the objective function.
Computations have been performed for two different shapes of the solid
region (circular and elliptic solid region) for two given boundary
conditions,

1. The concentration gradient is parallel to the flow field.

2. The concentration gradient is normal to the flow field.

For all the above-mentioned cases, the pressure gradient is varied while
keeping the values of h.

Several optimal microstructures of the porous media for parallel and
cross flow are shown in Figure (5.2 a-h). The dimensions of the
microstructures as shown in Figure (5.2- a—h) are tabulated in Table (5.1)

where l, h, a, b and 6 are the length, height, major axis, minor axis and

77

angle between length and height for the periodic cell.. The macroscopic
equation is solved for a domain of 1x1 m2 with a horizontal concentration
gradient of 1 moi/m3 and the pressure gradient has a unit of N/m-3
applied along horizontal direction. Periodic cells obtained when the
concentration gradient is applied parallel to the flow field are shown in
Figure (5.2 a—d). Periodic cells Obtained when the concentration gradient
is applied normal to the flow field are shown in Figure (5.2 e-f). The
applied pressure gradients are -2 in Figure (5.2 a, c, e, g) and -0.5 in
Figure (5.2 b, d, f, h) N/m3. Figure (5.2 a, b, e, t) represents the circular
and Figure (5.2 c, d, g, h) represents the elliptic solid region. From figure
5.2 we see the optimized results shows a square array in most of the cases.

Further comparison of the magnitude of the objective functions versus

pressure gradient is provided in Figures (5.3) and (5.4).

78

 

......OCOCC‘

...--.C---....-.--.-
\

(a)

\
.

 

'C.......-.‘
......O...

.3--€---------------‘
(C)

0......
I

e
o
o
o
I
J

   

'cooooooooo
......o-.

I

I
O
l

----..

" (e)

--- _ o

(f)

I

I
I
I
I
I
I
I
I
I
I
d

 

I I
I I
I I
I I
I I
I I
I I
I I
I I
I I
hr:arc-uo.uc-c.uc.ccd
(b)

 

I
I
I
I
I
I
I
I
I
l
d

. . (d)

   

Figure 5.2 a-h Schematic of the optimal periodic cells obtained after
optimization for the dimensions tabulated in Table (5.1)

79

Table 5.1 Dimensions of the periodic cells shown in Figure (5.2)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure (5.2) I/h a/h b/h 0
(a) 4 1.939 0.231 90
(b) 1.864 0.297 90
(c) 1.840 0.233 0.285 90
(d) 1.852 0.304 0.226 90
(e) 0.516 0.172 89
(f) 0.573 0.191 80
(g) 0.498 0.165 0.276 82
(h) 0.564 0.188 0.331 83

 

80

 

35'

 
    

30-

  

+ Circular solid region
’ 25 ~
+ Elliptic solid Region

 

I:

.9; 20 i

fl

0i

I: i i

E1 15 4

0‘ .

.29

‘5; ‘0 7

0‘ f

O 5 i

l“ T "“ I r I r ‘r i 0

-2.1 -1.3 -1.5 -1.2 -o.9 -0.6 -o.3 0

Pressure Gradlent

Figure 5.3 Computed values of the objective functions for a parallel flow
configuration.

81

 

    
 

 

 

0 I l
g .
:3 -0.1 — i
c + Circular solid region ‘
E . . . .
g '0 2 i + Elliptic SOIId region
a:
8
3- -0.3 e ‘
o .
-0.4 1
’
'0.5 I r l l l r i
-2.1 -1.8 -1.5 —1.2 -0.9 -0.6 -O.3 0

Pressure Gradient

Figure 5.4 Computed values of the objective functions for a cross flow
configuration.

82

It is seen in Figures (5.3) and (5.4) that the magnitude of

- - . ‘3 - - decreases for a decrease in ressure
dzsszpatlon drsperszve p

gradient. This can be explained with reference to equation (5.3) and (5.4).

Since the dissipation power directly proportional to the pressure term, for

a higher pressure gradient a higher value of dissipation power is expected.

5.3.2 Optimization of the “dispersive power” and the

dissipation power for a glven base dimension

In the previous section the variation of objective function with
pressure gradient is studied. In this section, I am interested in finding the
optimal pressure gradient that satisfies both the desired criterion;
maximize “dispersive power” and minimize dissipation power. The
optimal microstructures of the porous medium for the desired objective
function as described in section 5.1 for volume fraction of 0.3 for solid
phase and prescribed molecular diffusivity of 10'9 m/s2 are shown in
Figure (5.5). The macroscopic equation is solved for a domain of 1x1 m2
with a horizontal concentration gradient of 1 mol/m3.

The dependence of the step size and initial values in the optimization
algorithm were verified by varying those for two different cases. For case
1 the step size is set to 0.1 and initial aspect ratio is 1, the radius (major
axis, minor axis) is set to 5x104m and 0 to a value of 90° . For case 2 the

step size is set to 0.5 and initial aspect ratio is 1, the radius (major axis,

83

minor axis) set to 7x10'4m and 0 to a value to 90° . The iteration history
for the different cases are shown in Figure (4.10 — 4.12) for circular,

elliptic and rectangular solid region.

84

Iteration History

 

 

Objective function

 

 

 

 

o 20 4o 60 80 100 120
Iteratlon Number

Figure 5.5 Graph of the iteration history for the circular solid region with
a varied step size and initial conditions.

85

Iteration History

 

 

     

-o.54' . i
' :

—1.5 .

 

l:-CASE1 —CASE 2]

 

  

-2.5 « 1

-3.5 « ‘.

Objective Function

.45 ~

 

 

 

'5.5 I I I I I I
o 20 4o 60 80 100 120
iteration Number

Figure 5.6 Graph of the iteration history for the elliptic solid region with a
varied step size and initial conditions.

86

iteration History

0 - _ ,_ -!___ “7__W__ __.

 

 

---CASE 1 —CASE 2

 

 

Objective function

 

 

 

o 20 4o 60 so 100 120
iteration Number

Figure 5.7 Graph of the iteration history for the rectangular solid region
with a varied step size and initial conditions.

87

The dimensions of the microstructures are shown in Table (5.2)
where l, h, a, b and 6 are the length, height, major axis, minor axis and
angle between length and height for the periodic cell. The optimal values
of normalized dissipation power and normalized dispersion power

evaluated are tabulated in Table (5.3).

88

 

rucu-I-C...
I-..---.

E)

 

.---....---

(b)

---------'-----'—-‘

 

'O.........‘

(C)
Figure 5.8. a-c. Schematic of the optimal periodic cells obtained after
optimization with the dimensions tabulated in Table (5.2)

89

Table 5.2 Dimensions of the microstructures as shown in Figure (5.8).

 

 

 

 

 

Figure (5.8) l/h a/h b/h 6
(a) , 1.07 0.32 89
(b) 1.44 0.47 0.33 89
(c) 1.74 0.42 0.33 89

 

 

 

 

 

 

 

90

Table 5.3 Optimized value of normalized dissipation power and
normalized dispersion power for the microstructures as shown in Figure
(5.8)

 

 

 

 

F1 gure (5.8) S disipation S dispersive
(a) 0.01427 4-63
(b) 0.007658 4-96
(C) 0.000772 4-93

 

 

 

 

 

91

The optimal pressure gradient computed are

(0) _ (0) _
VXp - 8(10 6) N/m3, VXp = 800 6) N/m3 and
-0.7 VYP(O) -0.4

0
V Xp( ) =[8(10_6)

] N/m3 for Figure (5.8) a, b, and c respectively.
- 0.1

The average effective permeability and the average dispersion tensor

computed for the periodic cells shown in Figure (5.8) a, b, and c are

K—341(10‘5) 1 EC 2
‘ ' e 1.38 m’
C
_ 1 1.31
D=5.4(10 9) mzls
1.31 3

K—307(10‘5)l 86 2
‘ ' 33 1.81 m’

_9 1 1.31 2
D = 5.20(10 ) m /s,
1.31 3

-5 1 5 2
andK=1.91(10 )8 r m,

I”

_ 1 1.56
D = 450(10 9) m2/s,
1.56 3

where ac =1.4(10“ 2), £8 = 212(10‘ 2) ,er =1.61(10“) respectively.

92

Figures (5.9) and (5.10) show a contours of the longitudinal component of

the effective permeability and dispersivity

93

 

(C)

Figure 5.9 a-c. Contours of the longitudinal component of the
permeability (pUl 1) for the microstructures shown in Figure (5 .8)

94

 

(b)

11109

19
n I10
-1

(c)
Figure 5.10 a-c. Contours of the longitudinal component of the
dispersivity (211) for the microstructures shown in Figure (5.8)

95

...............

........
..........

...........

(a)

........

Illlli Ii

(b)

 

 

 

 

 

 

 

II .I l ililiiln
n A I I. II I. 0.11 I. Jill. I. I ..I: I I
A. I'OI .Al I ITIII I ILI‘..II Iii.
a. i D 4 oil .1 I ’i
IIIIIIIIIIIIIII
J
u
_
IIIIIIIIIIIIIII
7 Al 4.10 i o
I II I A Ii I 11.. II I. III I. I
_
I I I iii 0. t I I: I... 1_1
ii I

(C)

Figure 5.11 a-c. Plots of the velocity vectors in microstructure

corresponding to Figure (5.8)

96

5.4 Discussion

In this chapter a methodology to find the optimal microstructure that
maximize mixing andminimize energy loss using shape optimization is
presented

A goal is also to find a periodic cell that provides a better solute
dispersion while minimizing the energy loss for a fixed or a variable

pressure gradient. It is seen in Figures (5.3) and (5.4) that the magnitude

of 3 decreases for a decrease in pressure

dissipation — 3 dispersive
gradient. For a higher pressure gradient, a higher value of dissipation
power is expected compared to the dispersive power. For a flow normal to
the concentration gradient, the optimization process does not converge to
a finite value for a very low pressure gradient (01 N/m'3). This may be
due to the fact that the applied pressure gradient is not sufficient to
disperse the solute though out the domain. Also since because of the
nature of the problem the final results are highly dependent on initial

value and sometimes the solution is not the global solution. Furthermore

the sensitivity of 8 are approximated using

dissipation — S dispersive
finite difference method, which may cause some additional errors.

In the second analysis it can be observed from Table (5.3) that for the
microstructure with a rectangular solid region, the energy loss due to

dissipation is the least and that the solute dispersion is comparable to the

97

other microstructures. Since the pressure gradient is also included as the

design variable, it takes the advantage of a higher transverse effective

permeability (K22) and adjusts itself to minimize the energy loss. In other

words it minimizes the pressure gradient in the horizontal direction and
diverts the flow in vertical direction. This is also because the height was
fixed; whereas, other dimensions could vary.

The results obtained using shape optimization for minimizing
dissipation power and maximizing “dispersion power” ,in which the
pressure gradient is a design variable, can be compared with the
microstructure formed by a square array of same porosity. The
comparison shows significant improvements in reducing energy losses in
the fluid flow and increased mixing. Figure (5.12) shows the macroscopic

P —6
layout with an applied pressure gradient X x = 800 )
pry 04

obtained for elliptic cross section (Figure 5.8 (b)). The macroscopic
velocity field computed from the Darcy’s equation is compared at the
cross section shown by dashed line in Figure (5 .13). When compared with
the microstructure form by square array of same porosity, the optimized

microstructure results in a higher velocity field as shown in Figure (5. 14).

98

 

I V

pr I

Figure 5.12 Schematic of the macroscopic layout with applied boundary
conditions

0.05 ~- —-—-—-—— -- -—.——- -- --——— ——— r——w~ *—

0045 1 — MinO structure wrth elliptic cros section

- - Square array with same porosity.

0.04 4

0.035 4

.°

0

on
L

0.025 1

Velocity field (m/s)

 

 

 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Location (m)

Figure 5.13 Comparison of the velocity fields for a domain with a
microstructure formed with the optimized elliptic cross section with the
velocity field for a square array of same porosity.

100

The objective function in this chapter is set in such a way that it also
maximizes the “dispersive power” or mixing. Figure (5.11) shows higher
solute dispersion at time t for the optimized microstructure. The
concentration are compared for both the cases and shown in Figure (5.11

c).

101

  
 

   

After timet

(a) (b)

0.75 - :

Location (m)
0
01

   
 
 

 

 

0.25 -
o l l
0 0.5 1
Concentration (c)
(C)

Figure 5.14 a—c. Comparison of solute dispersion after time t = 0.02 sec
for the macroscopic domain formed with a microstructure of elliptic cross
section cylinders (a) obtained from shape optimization with the
microstructure obtained from a square array (b) of same porosity. The
magnitude of c is shown in right (c).

102

CHAPTER 6

TOPOLOGY OPTIMIZATION OF FLUID FLOW IN POROUS
MEDIA TO MINIMIZE DISSIPATION POWER

6.1 Introduction

Topology optimization is a layout optimization technique that was
originally developed to design mechanical structures (Bendspe and
Kikuchi, 1988). Its scope was rapidly expended to the diverse field such
as optics and acoustics (Bendsde and Sigmund 2003, Eschenauer 2001,
Jensen 2003, 2004) and recently it has been applied in the field of fluid
flow by Borrvall and Petersson (Borrvall 2003). Most of the above
mentioned works were mainly applied to macroscopic layout problems.
With further development of the homogenization theory, topology
optimization was used to find the layout of the microstructure or the base
cell of the material. Few such noted works are: Sigmund (Sigmund and
Torquato 1996) to design microstructure for material that yield negative
thermal expansion, Larsen (Larsen 1996) to design material that yields
negative Poisson’s ratio and Diaz and Bénard (Diaz and Bénard 2003) to
design material that matches prescribed elastic properties.

The formulation of the problem is given in the following section. The

Method of Moving Asymptotes (MMA) as proposed by Krister Svanberg

103

(Svanberg 1978, 2002) is used for optimization. The derivatives of the
objective function and the constraints with respect to the design variables
were computed using adjoint method as proposed by Olesen, Okkel and
Bruus (Olesen 2006). Figure (6.1) shows the flowchart of step by step

approach to find the optimal microstructure using topology optimization.

104

 

 

 

Pre-processing

 

1

 

Optimization

 

 

 

i
1

initial Guess
i

Solve for dissipation power
i—" using finite element
‘ analysis

 

 

 

 

 

 

 

 

i 51

 

 

Constraints analysis
l i

 

 

 

.5 Sensitivity analysis

: l

i Filtering the sensitivities
1

Optimization with MMA

i

t

 

 

 

 

 

 

 

 

 

 

 

 

Design variables update

 

 

 

T

     
 

 

Termination condition
achieved 7

Yes

Y

 

Post-processing

 

 

 

Figure 6.1 Flowchart showing the steps of the algorithm based on
topology optimization for designing microstructures.

105

6.2 Formulation of the optimization problem

In the analysis, the macroscopic domain (QM) is assumed to be
homogeneous with periodic microscopic cells (Qm ). The microscopic

cell has two distinct regions: a solid region (525) and a fluid region ((2

f)

such that .Qm = $23 UQ andQs m!) =0. The interface between the

f f

solid-fluid in each periodic cell is 1",; and the interface between two

periodic cells is denoted by Fee".

106

 

1flcell

 

 

 

 

Figure 6.2. Schematic of a periodic microscopic cell (0,, ). The interface
between the solid-fluid in each periodic cell is F51 and the interface

between two periodic cells is denoted by Teen“

107

6.2.1 Minimizing Dissipation power

The objective in this section is to minimize the energy losses or the
dissipation power.

Find the optimal layout for the periodic cell that will

.8 . . .
“““imize' dzsszpatzon

subjected to: volumetric and geometric constraints

Q (6.1)

Constraint: 3 2 V
52

 

frac
m

where 05 is the volume of the solid region, Om is the total volume of

the cell, and mec is value (0<mec<1) specifying the minimum solid

volume. The constraint in equation (6.1) restricts the minimum volume of

solid region to mec.

The governing equations here are the homogenized Navier-Stokes
equation and the incompressibility constraint as derived in Chapter 2. To
generalize the above governing equation for both solid and fluid phase,
the fluid flow through the solid region is subjected to a friction force,
which is proportional to the fluid velocity. The governing equations for
the zeroth order problem are given by

x x m

XP

108

and

v 11(0) :0 in o (6-3)
X m

where V (mis the pressure gradient across the macroscopic domain

X P
and is applied as a source term (a known quantity). These governing
equations are valid for both solid and fluid regions. The penalty term a,
which is a function of the design variable p, allows to set the velocity
equal to zero in a solid domain. The design variable p controls the phase

of the medium as follows;

IVerf
Following reference Borrvall (Borrvall 2003) and Olesen (Olesen 2005),

{0 if x e as (6.4)
120) =

the penalty term 01 and design variable p are related by the convex

interpolation.

a( )=a +(a — ____q[1-p] (6'5)
'0 — min max min q+p

where q is a real and positive parameter used to tune the shape of the

penalty term and taken as 0.01. In this work “min is taken as zero and
amax is a large number (105). The dependence of a on q is shown in

Flgure (6.3). For q = 1, 1t gives a lmear plot from amax to amin over the

range of p. For lower value of q (0.01), 01 value is close to amin for p

varying from 0.2-1. This approach is done to penalize the grey material

109

formation, which is a mixture of solid and fluid and promote either the

fluid or liquid phases.

 

1.E+5

 

 

8.E+4 —

 

 

 

2.E+4 ~

Penalty function (a)

 

 

 

 

0.E+0

0 0.2 0.4 0.6 0.8 1

 

 

Design variable (p)

 

Figure 6.3 Plot showing the dependence of a on q.

110

The boundary condition are that of periodic flow through the boundaries

of the cell and no slip boundary condition at the solid fluid interface of the

inclusion.

“(0) =u(l) = ...... =0 onl" (solid-fluid interface), (6'6)

u(0),u(1)... and p(0),p(1).... are Qm-periodic.

The 3 . . . is the dissipation energy for the fluid in the unit cell
disszpatlon

and is derived from mechanical energy balance. Just as before, it is
. . . . . (0) (O) .
obtained by mult1ply1ng equat1on (6.2) With v , where v 18 any

weighted periodic velocity satisfying equation (6.3) and boundary

condition as shown in equation (6.6), and integrating over the entire

domain After setting vm) = “(0)
disszpation x x
Q
m
In dimensionless form equation 6.8 can be rewritten as,
. . . = j v 6(0) ;(v 13(0))ng (6.9)
dzsszpation x x
Q
m
where fim) is the dimensionless form of “(0) related by

l3(0) = .1(0) ,<u(0)>.

111

6.2.2 implementation Issues

To promote a discontinuous microstructure and avoid cases of a
continuous solid region, the layout is divided into two initial zones. Zone
1 is formed only by fluid material and forms the outer boundary of the
periodic cell, whereas Zone 2 is formed on both solid and fluid material as

shown in Figure (6.4). In Zone 1 the a value is taken as zero.

112

 

Zone 1
)1

 

 

 

 

Figure 6.4 Schematic of a periodic microscopic cells with a liquid zones
imposed on Zone 1 above to avoid a continuous solid region between the
cells.

113

6.3 Example Problems

The optimization problem was solved using commercial software's
Comsol and Matlab. “The governing equations were mentioned in earlier
sections were solved to compute the objective function and the sensitivity
were using Comsol. Matlab was used to run the optimization process.
Work done by Olesen (Olesen 2005) is used here to compute the
sensitivity. The domain of Zone 2 was discretized using an unstructured

triangular mesh of 50 elements along each boundary as shown in Figure

(6.5).

114

. - .1',
4.74)“ Ir

. ,; . a .
7+?i .. 39381
. ,5

‘ J

1%

"s

{i
O c'.
8’11
..u . . .u .
r:

r
11

.1- I q; p
N. .3. '
an .
‘3.“7.‘ V
#11521
M

'3}:

5"

5'31»; '
~22 ‘ gt, 1%.?” Ar.
... $13.8?fi
L $§§°3"I‘¢“0-Jafir-
& -Wli‘v"- 0.- fi’ '

t 0_ Jr ,3! .
.fifii‘éfifaiifim
’35,} ..'-.-'-'~'h' 'KiL 5““; - 3 It"
fiiflr‘tss‘vfi-‘f i1»-
.— ‘59 ”life?“ '3
.1 #3:?fécfir3.’ '
'12 . -

5

3
rise-fgwfifi‘
. “1.2mm?

“tiara-rm:
Rind; a.

3;,
. - -.....

{5

a

V

v

1111?- 4:1

are unit A;

I:

1
may
..- u .

:-‘.
IV
1:

‘1... ..-

I?"
‘3
.1 I r.
5 V

'vwv.
{l-

.. .
535,1 7

15):»;

J u"...
fin‘bz‘ll"
1‘ 9‘4

11. -

I
P .
rt“

1.

It

8

4...:

2:1):
_ \H

'3”

'-. 2

51

7
.—

IQ»

‘-

A!”
'l. .

_ 1:
1 REESE?

"J
s

4 .
'H
‘l

J_r

I

'3

L

-.,-_. a;
'K'
.

m

.36..»x
83' 1

7:11.!

”. 1?: ‘-'E
.t’e’.

finite element analysis.

115

r .
.‘pg.
IL

251..

1.
’4
"I"

;- I. 1-

F

. fl
433

83.1"

$1
«$6.16 If;

«a; .

.wxa.

.r—vgn

.i'lf‘}

a" .

j.
r on 'n a;
l.eA":u$.ib" .

r
1.
1%..

A
.v.
”III-n.1- ’1

:.V
41'; '

In

.- "‘-
M1.

I Is- QWVAJA'I

AL

1 5’1 "ya:
a. ailing- if):

 

Figure 6.5 Schematic of the unstructured triangular mesh used for the

6.3.1 Minimizing Dissipation Power

In this section, the optimized microstructures for the porous media for
different boundary condition are computed. Figure (6.6) and (6.7) shows
the microstructure obtained when a macroscopic pressure gradient of
1N/m2 is applied in the horizontal direction with volume fraction of solid
region is limited to 0.5 and 0.7. In Figure 6.6 we see that for lower volume
fraction of solid the microstructure forms a elliptic solid region, whereas
for higher volume fraction it forms a rectangular shape solid region with

curved edges.

116

Figure 6.6 Schematic of microstructure layout for volume fraction 0.5 and
with macroscopic pressure gradient of lN/m2 in horizontal direction

 
 

117

Figure 6.7 Schematic of microstructure layout for volume fraction 0.7
with macroscopic pressure gradient of lN/m2 in horizontal direction

118

Figure (6.8) — (6.14) shows the microstructure obtained when
macroscopic pressure gradient of lN/m2 is applied in both horizontal and

vertical direction with volume fraction of solid region is limited to 0.1, to

0.8.

119

 

Figure 6.8. Schematic of the microstructure layout for a volume fraction
0.2 with an equal macroscopic pressure gradient of 1N/m2 in the
horizontal and vertical directions.

120

 

Figure 6.9. Schematic of the microstructure layout obtained for a volume
fraction of 0.3 with an equal macroscopic pressure gradient of 1N/m2 in
the horizontal and the vertical directions

121

Figure 6.10. Schematic of a microstructure layout for volume fraction 0.4
with equal macroscopic pressure gradient of lN/m2 in the horizontal and
the vertical directions.

122

Figure 6.11 Schematic of the microstructure layout for a volume fraction
0.5 with equal macroscopic pressure gradients of lN/m2 in the horizontal
and vertical directions.

   

 

 
   

123

Figure 6.12. Schematic of the microstructure layout for a volume fraction
0.6 with equal macroscopic pressure gradients of lN/m2 in the horizontal
and vertical directions.

124

Figure 6.13 Schematic of the microstructure layout for volume fraction of
0.7 with equal macroscopic pressure gradients of IN/m2 in the horizontal
and venical directions.

125

Figure 6.14. Schematic of the microstructure layout for a volume fraction
0.8 with equal macroscopic pressure gradient of 1N/m2 in the horizontal
and vertical directions.

126

6.4 Discussion

In this chapter I tried to find an optimal microstructure to minimize
the energy losses using topology optimization.

The results show that the elliptic cross section for a lower volume
fraction has the least energy loss. This also agrees with the results
obtained from shape optimization. For a higher volume fraction, the
optimized layout obtained appears as a rectangle with smooth comers.
This validates the previous argument made in Chapter 5 that the
microstructure may form a channel flow. When the macroscopic pressure
gradient is applied along both the horizontal and vertical direction, the
optimized layout forms the shape of leaf. In this layout the specific
surface which is expressed as the ratio of the pore surface area to the total
volume of the periodic cell is less and hence increasing the permeability.

When optimized microstructures obtained from topology optimization
were compared with the square array of same porosity, it shows the same
trend. When the macroscopic flow is compared for optimized
microstructure and square array formed by circular solid region with same
porosity, a higher magnitude of velocity is obtained for the optimized
microstructure as shown in Figure (6.6). Figure (6.15) shows the

macroscopic domain and the applied boundary condition

Vpr 1
where = . Figure (6.16) shows the macroscopic velocity
V X py O

127

profile for both the optimized microstructure and the un-optimized

microstructure for Vfrac 0.5 at the dashed line shown in Figure (6.16).

Here the un-optimized microstructure is formed by square periodic cell
with circular solid region of identical porosity as the optimized

microstructure.

128

pr'i

 

I pr 'j I
Figure 6.15 Schematic of the macroscopic layout with applied boundary
conditions

129

0.3 .. ~
— Velocity profile for optimized microstructure

0.25 7 — - For Square array of same POIOSitY

.0
m
L

 

Velocity Field (m/s)
0
6‘.

 

 

 

0 0.2 0.4 0.6 0.8 1
Location (m)

Figure 6.16 Comparison of the velocity profiles in the macroscopic
domain for a microstructure formed form a cross section obtained from
topology optimization with the velocity field obtained for a square array

of same porosity

130

CHAPTER 7

TOPOLOGY OPTIMIZATION OF SOLUTE TRANSPORT IN
POROUS MEDIA TO MAXIMIZE “DISPERSIVE POWER”

7.1 introduction

Solute transport in porous media and finding the optimal layout using
topology optimization is the subject of this section. Topology
optimization is a layout optimization technique that was originally
developed to design mechanical structures (Bendspe and Kikuchi, 1988).
Topology optimization was then used to find the layout of the
microstructure or the base cell of the material. Examples of such works
include Sigmund (Sigmund and Torquato 1996) who propose to design
microstructures of materials with negative thermal expansion, Larsen
(Larsen 1996) to design material that yields negative Poisson’s ratio and
Diaz and Bénard (Diaz and Bénard 2003) to design material that matches
prescribed elastic properties.

The problem is formulated in a similar manner to a standard
optimization problem, in which the Method of Moving Asymptotes
(MMA) is used for optimization. The derivatives of the objective function
with respect to the design variables were computed using an adjoint

method as proposed by Olesen, Okkel and Bruus (Olesen 2006). Figure

131

(7.1) shows the flowchart of step by step approach to find the optimal

microstructure using topology optimization.

132

 

Pre-processing

 

 

 

!

Optimization

 

 

 

 

.'
Initial Guess

I

Solve for dispersive power
. using finite element
analysis

 

 

 

 

 

 

 

 

 

V

Constraints analysis

 

 

 

 

i
Sensitivity analysis
1 i
Filtering the sensitivities

i

1

Optimization with MMA

l

Design variables update

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

T
i

   

 

   

Termination condition
achieved ?

  

Yes
1

 

Post-processing

 

 

 

Figure 7.1 Flowchart of the topology optimization methodology
employed

133

7.2 Formulation of the optimization problem

In the analysis, the macroscopic domain (QM) is assumed to be
homogeneous with periodic microscopic cells (Qm ). The microscopic

cell has two distinct regions: a solid region (525) and a fluid region (I)

f)

such that Qm = as UK) and as n Q = 0. The interface between the

f f

solid-fluid in each periodic cell is F31 and the interface between two

periodic cells is denoted by FCC“.

134

 

1flcell

 

 

 

 

Figure 7.2 Schematic of a periodic microscopic cell.

135

7.2.1 Maximizing “Dispersive power”

The main objective here is to maximize mixing or “dispersive power”.

max. dispersive

subjected to: volumetric and geometric constraints

(23 (7.1)
a— 2 Vfrac
m

where 05 is the volume of the solid region, Qm is the total volume of

the cell, and mec is value (0<mec<l) specifying the minimum solid

volume. The constraint in equation (7.1) restricts the minimum volume of

solid region to mec.

The governing equation for the solute transport along with the Stokes
equation and the incompressibility constraints are for the volumetric

concentration c of the solute in the fluid region expressed as

<..>.v,.<o>.vx .[D.vx.<o>] (7.2)

where <u>=<u (0)>+<u (1)> and for weak inertia flow <u(l)>=0

0)

(Mei and Auriault 1991). u( is computed from the homogenized

N avier-Stokes equation and the incompressibility constraint as derived in
Chapter 2.

D in Equation (7.2) is the effective dispersion tensor (Mei 1992)

136

expressed as
D = 0(2) (73)

where Z is “characteristic dispersivity tensor” given as,

11111+120era-((1111118)

And D is the molecular diffusivity and N is any vector satisfying the

(7.4)

equations
V(D(1+v ...»-..«nv ...-..«n (7.5)
x x x
— <u(0)>/n — a(va +01an
(7.6)

where n is the porosity, and

Dll+V -N)-n=Oon F.
x
N is Qm -periodic and (N) = 0.

(0)

Additional governing equations to solve u are the homogenized
Navier-Stokes equation and the incompressibility constraint as derived in
Chapter 2. For generality the fluid flow through the solid region is
subjected to a friction force, which is proportional to the fluid velocity.

#szum) _pr(1) :17 pm) _ cm(0)in Om (7.7)

X

V -u(0) :0 in Q (7'8)
x m

137

where V (0)1s the pressure gradient across the macroscopic domain

X P
and is applied as an source term. The above governing equations are valid
in both solid and fluid phases since the penalty term 01, which is a function
of the design variable p, is used to extend the domain of validity to all the
cell domain. The design variable p controls the phase of the medium as

follows;

OiferS (7.9)
p(x)= lifxeflf

Following reference Borrvall (Borrvall 2003) and Olesen (Olesen 2005),
the penalty term a and design variable p are related by the convex

interpolation.

1— (7.10)
mp) E“min +(“max — min ‘11 P]

q + P
where q is a real and positive parameter used to tune the shape of the

penalty term and taken as 0.01. In this work “min is taken as zero and
. 5
“max rs a large number (10 ).

The boundary conditions are

“(0) = “(1) = ...... =0 onl" (solid fluid interface), (7'11)

(0) (1)

u ,u ...andp(0),p(1)

are Qm -periodic.

. . is the normalized value of “dispersive power”,
disperszve

which is computed by multiplying equation (7.2) with do), and integrating

138

over the entire domain, the following expression for the “dispersive

power” is obtained after using Green’s theorem:

T (7.12)
(It) I VXc(O)-c(0)dfl= [VXC(O)-D-[V C(O)] d9

9 Q X
m m

The right hand term is the dispersive power term.

T (7.13)
. . = iv C(O)'D- V cm) do,
dzspersrve Q X X
m
In dimensionless form, equation (7.13) can be written as,
A . T (7.14)
8. . = iv 6(O)-D-V 5“” d9,
disperszve Q X X
m
(0)

where E and D are dimensionless form of can and D related by

5“” = C(O) /C and 1‘): D/D.

7.2.2 Implementation issues

A factor that had to be considered in formulating the problem for the
optimization Is that the microstructure should be discontinuous as shown
in Figure (7.2) to model pores in a porous material. To avoid cases of
continuous solid regions, the layout is divided into two initial zones. Zone
1 is formed only by fluid material and forms the outer boundary of the
periodic cell, whereas Zone 2 is formed on both solid and fluid material as

shown in Figure (7.3). In Zone 1 a value is taken as zero.

139

 

Zone 1
i

 

 

 

 

Figure 7.3 Schematic of a periodic microscopic cells with the zones to
avoid continuous solid region

140

7.3 Example Problems

The optimization problem was solved using Comsol and Matlab. The
objective function and the sensitivity were computed using Comsol. Work
done by Olesen (Olesen 2005) is used here to compute the sensitivity.
Total nodes in Zone 2 were approximately 30 elements along each

boundary.

141

 

 

 

 

 

 

 

._ . , ,
. ... 1,.,1.- .~~' ..f-;
. - . , ‘ .wr‘.
. « be
e 1,, 7f. .
‘ 1'
-.. ._ .,. 1.
. J.
\ v-‘- .,lrv ,"(«
1-1‘ '
t - i
. -. 4‘.
. . _
717,-, e. I
,‘ > ' .
. 1
\ V}
7 ~. ‘ w
- _. 1 5 L-.
.1. .
1
’. ja'._ .
- .1 , 1 .
, . .
e \ ,e .7 ,4 x
l '-
. ,.
".4 s ’1 -.
‘ . g i
4 1,: ' r
e ,‘ ‘t
r ‘ r, .i
1 “
1 e '1 ‘ V"
.1. 4
‘_ > ‘ 1 1 1
. y
o . "‘
‘r'v,_ r. ,Y
, . , . .‘_~. ‘
\ 1
~,'A.' 1‘ ITt '
. J
v V ' t V v c

 

 

Figure 7.4 Schematic of the unstructured triangular mesh used in the finite
element analysis.

142

7.3.1 Maximizing “Dispersive power”

In this section, the optimized microstructures for the porous media that
minimize dispersion are shown. Figure (7.5) shows the microstructure
obtained when a macroscopic pressure gradient of 1N/m2 is applied in
both the horizontal and vertical directions with a volume fraction of the
solid region limited to 0.7. The concentration gradient of l moi/m3 is
applied in horizontal direction and the objective function is set in such a
way that the microstructure should maximize the mixing in the horizontal
direction.

A sharp change in the microstructure is seen when the volume fraction

is reduced to 0.3 as shown in Figure (7.6).

143

Figure 7.5 Schematic of microstructure layout for volume fraction 0.7
with equal macroscopic pressure gradient of lN/m2 in the horizontal and
vertical directions and a concentration gradient in the horizontal direction.

144

L...L.....L..L...

L...L...L.....L_..

L_.L..L..L..

Figure 7.6 Schematic of the microstructure layout for a volume fraction
0.3 with equal macroscopic pressure gradient of 1N/m2 in the horizontal
and vertical directions and a concentration gradient in the horizontal
direction

145

7.4 Discussion

In this chapter we tried to find an optimal microstructure to maximize
mixing using topology optimization. Due to the limitations of the solver
used, the results for maximizing dispersive power cannot be solved for a
fine-mesh . The results obtained and shown in this chapter provide a hint
of the appearance of the porous medium. In the results shown in Figures
(7.5) - (7.6) the inflow in at negative 45 degree to the horizontal and the
aim was to distribute the solute in the horizontal direction. In Figure (7.6)
the microstructure diverts the flow as much as possible to the horizontal
direction. For a lower volume fraction the solid materials form a wall
around the boundary of zone 1 and zone 2 to restrict the flow in vertical
direction. For a higher velocity in horizontal direction the longitudinal
dispersion tensor will be higher and hence maximize the solute flow in

that direction.

146

CHAPTER 8

THEORY FOR MULTIPHASE FLUID FLOW IN POROUS MEDIA

8.1 Introduction

Practical applications of porous materials often involve two or more
fluids. For example, use of porous material as a filter in the petroleum
industry often involve oil and water. As discussed in Chapter 1, naturally
occurring porous materials are heterogeneous there have been many
theoretical attempts to deduce the phenomenological equations by starting
from the micro-scale based on the idealized models of the microstructure.
In this chapter we review the existing theories for deriving the Stokes
equation for two phase flow in the micro-scale and derive the
phenomenological equations that describe the macroscopic behavior of the
porous media. The governing equations are developed using mixture
model. Shape optimization is used to determine the optimal
microstructure for porous media that will minimize the dissipation power,

for a given flow condition.

147

8.2 Derivation of Effective Permeability Tensor using Mixture
Model

In two phase flow problems, the average velocity of one phase is
typically distinctly different from the other phase (Kleinstreuer 2003).
Also, one or more physical properties such as density or viscosity of each

phase distinctly differ in magnitude. The macroscopic domain ((2 M ) is
assumed to be homogeneous with periodic microscopic cells (9m ). The
microscopic cell has two distinct regions: a solid region (528) and a

mixture fluid region (Q ) such that 9m =QSUQ

fmix fmix

and as m 9 = 0. The Mixture model computes the average behavior

finix

of a two phase flow field as a single phase flow that is rather general and

useful as shown in Figure (8.1), where the mean density (pf mix) can be

expressed as a function of the volume fraction Vf (Vf = V2/V). The

effective mixture density in terms of volume fraction is given as

8.]
meix =vfpf2+(1—vf)pf1 ( )

where the indices mark the individual phases i=1 (carrier fluid) and i=2

(dispersed phase). For example pf 1 is the density for fluid phase 1 and

f

,0 2 is the density for fluid phase 2.

148

Pf mix, if mix

 

 

Figure 8.1 Schematic of the process of representing a two phase flow
with an equivalent a single phase flow model.

149

 

The average steady-state mixture continuity and momentum equations

(Ishii, 2006) are given as

 

. f _ . (8.2)
V (pmixumix _0 1n Qfmix
f 2 _ _ (8.3)
#muV'flmx VKMx—
p .u .-Vu .+V- ii in!) .
mix mix mix f . 21 fmrx
(l—ijp mix

where “mix is defined as

 

plf (l-vf) pécvf (8'4)
11 . =———u + 11
mix f l f 2
pmix pmix

Gravity and mass transfer between two fluids are ignored. 112, is the slip

velocity between two fluid phases. It is also assumed that the solid phase
is chemically inert. For Newtonian fluids the momentum transfer due to
shear stresses within fluid is negligible as compared with the momentum

transfer to the solid matrix (Allen 1985). So

f f (8.5)
VfP1P2 ~ 2 =0

(1 —vf)p,{,,x “2'

and equation 8.3 reduces to

V .

 

f V2u —Vp 2pf 11 u ' (8'6)

. . . . . - . in Q .
mix mix mix mix mix mix fmrx

,u

The perturbation expansion for “mix and Pmix is as follows,

150

“(0) mm + 211+”) é.u3 (3) (8.7)

“mix: umix+ wmix umix+ umix
p(0) (1) £2 (2)+ £3 (3) (8-3)
p mix pmix +6pmix+ p mix p mix

The differential term related the two length scales as a function of a

V=Vx+eVX (8.9)

Substituting perturbation expansions for um,r and pm in equation 8.6 and

8.2
#f (V ”V “H (0) a1(1) +£2u(2) +5 3‘10) )_ (8.10)
mix mix+ mmix ix mix

(v ”Vx p+(0)£p(1) 2 p(2) ”316(9)—
pmix+ Epmix+ pmix Mix

p’iix(u£?;+ 811(1). +£2u(2.) u+e3ugzxj
(V +eV {u (0)+ sum +.92u(2 .) +8 3u(3) )
x X mix+ almix mix mix
and
f (0)+ a1+(1) 211+(2) 3(3)
(Vx+€VX)'(pmix(u mix+ anmix+ umix+ EBumixD: 0 (8°11)

Assuming the mixture density is not a function of macroscale X, from

Equation 8.10 and 8.11, at orders from 0(80) I get

f 2 (0) (1) (0) - (8.12)
’umixVx umix prmix= Xpmix 1n Qfmix
and
, (0) _ - (8.13)
VX umix — 0 1n Qfmix

151

(0) (1) (0) (1) _ 1 - (0)
and “mix ’umix and pmix ’pmix are 0m per1od1c. p

is only a function of the large scale X. The needed boundary conditions

are
“(0) = “(1) = ...... = 0 on I" (solid fluid interface), (8'14)
mix mix

In equation 8.10 and 8.1 1,8 = l/ L <<1 is the ratio of two well separated

length scales: the micro scale I=O(x) and macro scale 1,1000. Here the
macrostructure is assumed to be homogeneous consisting of periodic
micro cells.

Equation 8.12 relates the microscopic fluid flow with the macroscopic
pressure gradient. This is similar to the approach taken by Darcy to relate

the fluid flow in porous media for a given macroscopic boundary

condition. umixw) and pmixm can also be expressed in terms of Pmix (0)
from Darcy’s equation,

umixm) = —Umix . VX pmixm)’ (8.15)
and

pmixa) zamix' Xpmixw)’ (8.16)

where Umix is the characteristic mixture velocity and amix is the

characteristic mixture pressure can be obtained from the solution of

8.17
f.V2U.—Va.-I ()
mlx x mix xmlx

,u

152

and

V U , :0. (8.18)
x mlx

The boundary conditions applied are no slip and periodic boundary

conditions for Umix and amix,

The effective permeability of the porous media is computed from

K = ,uf (Umix> (8.19)

mix mix

Equations (8.17 - 8.18) are similar to the single phase equations (2.9 —
2.10), where the velocity, pressure of a single phase is replaced by a
mixture velocity and pressure. This equation is solved to the effective
properties of the mixture fluid.

Shape optimization problem is used to find the optimal microstructure
where various shapes described parametrically are studied. Since the
governing equations are very similar to the governing equation for single
phase fluid, this work can be taken as an extension of the shape
optimization work done in earlier chapters. The Method of Moving
Asymptotes as proposed by Krister Svanberg is used for optimization. The
derivatives of the objective function with respect to the design variables

were computed using finite differences.

153

8.3 Identifying the Periodic Cell

. In the analysis, the macroscopic domain (QM) is assumed to be
homogeneous. The microscopic cell (am) has two distinct regions: a

solid region ((23) and a mixture fluid region (Q ix) such that

fm

Qm=QSUQ andflan =0. P, l and L are the

fmix fmix

characteristic pressure, micro length and macro length scales. Two scales

are introduced: a small scale (x 2: 0(1)) and a large scale (X z 0(L)). The

multiscale coordinates are related by X = 8x in the asymptotic
expansion.

The homogeneous macroscopic domain is formed of periodic
microscopic cells as shown in Figure (8.2). Each cell is identified by the
geometric parameters such as l, h, a, b, and 9 . L is the cell length, h is the
cell height, a is the major axis or length of the solid region, I) is the minor
axis or height of the solid region, and 6 is the angle between I and h in
each cell. Three different shapes for the solid region such as circle (a = b),

ellipse (Figure (8.3)), and rectangle (Figure (8.4)) are considered below.

154

 

 

 

 

 

 

Figure 8.2 Schematic of a homogeneous macroscopic domain studied
consisting of periodic microscopic cells.

155

 

Figure 8.3 Schematic of periodic microscopic cells with an elliptic solid
region

156

 

I I
I I
I I
I I
I . I
I g . I
I I
I 0 I
I I
L - - - - - --------- - - - - - J l

(C)

Figure 8.4 Schematic of periodic microscopic cells with a rectangular
solid region

157

8.4 Formulation of the optimization problem

In this section the aim is to find a microstructure of the periodic cell that
minimizes the dissipation power for multiphase fluid flow with the added

constraint of volume fraction using shape optimization.

8.4.1 Minimize dissipation power for multiphase fluid flow

The complete formulation of the optimization problem is stated as

follows:

(0) and V p(0) that

Find6, h, a and b for a given value of l, VXp Y

will

minimize: 3 . . .
disszpanon

subjected to: volumetric and geometric constraints

52

Constraint 1: ——S— 2 V
Q frac

m (8.20)

Constraint 1: a < i
2
h sin 0

Constraint 1: b <

where 05 is the volume of the solid region, 0m is the total volume of the

cell, and mec is value (0<mec<1) specifying the minimum solid volume.

The first constraint in equation (8.20) restricts the minimum volume of

158

solid region to mec. The next two constraints are imposed to insure the

continuity of each phase.

. . . is the normalized value of dissipation power for the
dzssrpatzon

mixture fluid, which is computed by multiplying equation (8.12) with

(0) (0) . . . . . . .
vmix , where vmix 18 any periodic mixture veloc1ty satisfying

equation (8.15) and the periodic boundary condition, and integrating over

the entire domain

d. . . =,uf, IV u “(0):(V v .(0))TdQ (8.21)
13571761th)? "'1le x MIX x mlx
m

where the double dot operator (z) is defined here as a : B = aij'Bji'

8.5 Example Problems

The results shown in this section are all for a limiting volume fraction of
0.3 and fixed value of 0.002m for l. Mixture of water and crude oil are
solved here for volume fraction of 0.2. The physical properties of water
are density 1000 Kg/m3 and dynamic viscosity of 0.001 Pa.s. Physical
properties of crude oil are density 915 Kg/m3 and dynamic viscosity of
0.02 Pa.s. The properties of water and crude oil are taken at 60° F. The
macroscopic equation is solved for a domain of 1x1 m2. The applied

pressure gradient is 0.1 N/m'3 applied along X axis.

159

8.5.1 Minimize dissipation power for multiphase fluid flow

The optimal microstructure of the porous medium for the desired
objective function as discussed in section 8.3 for volume fraction of 0.3
for solid region is shown in Figure (8.5). The dimensions of the
microstructures are shown in Table (8.1) where l, h, a, b and 6 are the
length, height, major axis, minor axis and angle between length and height
for the periodic cell. The dependence of the step size and initial vale were
in the optimization algorithm were verified by varying it those for two
different cases. For case 1 the step size is set to 0.1 and initial aspect ratio
is 1, radius (major axis, minor axis) set to 5x104m and 0 value to 90°. For
case 2 the step size is set to 0.5 and initial aspect ratio is 1, radius (major
axis, minor axis) set to 7x10'4m and 0 value to 90°. The iteration history
for the different cases are shown in Figure (8.5 — 8.6) for circular, elliptic

and rectangular solid region.

160

Iteration hlstory

 

\l
O
4

O)
O
l

8

 

.5
O

'1 F —CASE1 ---CASE2

Obiectiye Function
0

N
O

—b
O
l

 

 

 

 

 

O

o 20 4o 60 80 100 120
Iteration Number

Figure 8.5 Schematic of the iteration history for a circular solid region
with varied step size and initial conditions.

161

01 O) \l
O O O
1 1 J

3s
0

Objective function

10*

 

iteration History

 

 

l—CASE1 ---CASE2

 

 

 

 

5 10 15 20 25 30
Iteration Number

35

40

Figure 8.6 Schematic of the iteration history for an elliptic solid region

with varied step size and initial conditions.

162

Iteration history

 

 

 

 

 

 

350 1 — ~——.—- —— -— ———— —— ——— — —— —— .—-—————~ — ——
300
I:
.2 250
up
0
5 200
IL '1
2
0§ 150 ‘
3100 «
o
50
0 T l r l r ——_i
0 1o 20 30 4o 50 60

iteration Number

Figure 8.7 Schematic of the iteration history for a rectangular solid region
with varied step size and initial conditions.

163

 

-

I I
I I
I I
I I
I I
I I
I I
fiounoocoucol

(a)

'C---------.
0 t
I t
I I
I

(b)
too-0.0.0...

 

 

L a a n a a u a .- u o
(C)
Figure 8.8 a-c. Schematic of the optimal periodic cells obtained after
optimization

164

Table 8.]. Dimensions of the microstructures shown in Figure (8.5)
Figure (9.2) M a/l b/l t9

 

 

 

 

 

(a) 0.73 0.21 84.4
(b) i 0.39 0.26 0.12 82.7
(c) 0.96 0.26 0.14 81.9

 

 

 

 

 

 

 

165

The average effective permeability computed for the microstructures

.. 0.1358 0
shown in Figure (8.8) a, b and c are<K . >=1(10 6) m2,
mix 0.0303 0
—6 0.4045 0 2
<K . >=1(10 ) m and
mix 8.5826e—6 0

-6 0.7727 0 2 . .
<K . >=l(10 ) m. Figure (8.9) shows the veloc1ty
mzx 0.0072 0

distribution for the mixture fluid.

166

 

 

x10
i575
4.3
o
(a)
x10—4
33.4
1.9
i‘f‘“ “' /
i... _ J
0
x104 0”
H30
(C)

Figure 8.9 a-c. Schematic of velocity distribution of fluid mixture for
microstructures shown in Figure (8.8)

167

x10'2
:1 3.45

 

(a)

x10'2
H 3.18

 

 

-8.24

(C)

Figure 8.10 a-c. Schematic of pressure distribution of fluid mixture for
microstructures shown in Figure (8.8)

168

The optimal values of normalized dissipation power evaluated are
tabulated in Table (8.2). The dissipation power tabulated is normalized
with the dissipation power for elliptic solid region. A possible reason for
the high dissipation power observe in a rectangular solid region is
because of the presence of a sharp edge. Also due to the flat surface the

stagnation point is spread over a area which adds to the energy loss.

169

Table 8.2 Optimized value of normalized dissipation power and

normalized dispersion power for the microstructures as shown in Figure
(8.8)

 

 

 

 

Flguré (8'8) 8disipation
(a) 2.85

(b) l

(c) 9.83

 

 

 

 

170

 

 

 
 

8.6 Discussion

In the last two chapters a theory that governs the multiphase fluid flow
in a porous media is presented and the partial differential equations are
used in an algorithm to design porous media. The governing equation
developed using mixture model gives a very approximate behavior of the
multiphase fluid flow in a porous media and do not allow to show phase
separation since the divergence of the slip velocity was set to zero.
Commercial software Fluent was used to solve the differential equations
using the finite volume method. After optimization, the microstructure
with an elliptic microstructure shows the least energy loss. From Figure
(8.8) it can also be added that for microstructure with elliptic solid region
the velocity profile shows less variance as compared with other two
microstructures. For a microstructure with rectangular solid region, the
velocity profile is more that of a channel flow.

When the macroscopic flow is compared (figure 8.12) for optimized
microstructure obtained from shape optimization with square array formed
by circular solid region with same porosity, a higher magnitude of
velocity is obtained for the optimized microstructure as shown in Figure
(8.8). The mixture velocity is obtained from the Darcy’s equation solved
at the macroscopic domain shown in Figure (8.11). The applied boundary

V

. . X p x 0.1
condition where V = .
X p y 0

171

 

Vpr

Figure 8.11 Schematic of the macroscopic layout with applied boundary
conditions

172

4.005051 - ,,,, ”MM— ___, — , M—.— i

 

 

 

 

 

 

 

— Velocity profile for optimized microstructure
- - -Velocity profile for square array of same porosity
3.00E-05 «
t
I
Q 2.00E-05 -
o
>
1.00E-05 4
0.00E+oo . . . .
0 0.2 0.4 0.6 0.8
Location (m)

Figure 8.12. Comparison of velocity field in a macroscopic domain for
microstructure form with the cross section obtained from shape
optimization to minimize energy loss for multiphase flow with that
formed by square arrays of the same porosity

173

CHAPTER 9

SUMMARY AND CONCLUSION

In this work, two methodologies based on the theory of
homogenization, combined with two material design are presented, one is
based on shape optimization, the other on topology optimization. The
theory of homogenization is used to obtain the macroscopic equations
and the microscopic equations. The effective properties such as
permeability and dispersion of the porous media are computed using the
finite element method and compared with experimental results. The
method of moving asymptotes is used to find the optimal periodic cell that
will satisfy a given criteria such as the maximization of some function of
the effective properties, minimization of energy loss, or maximization of
solute mixing using both shape optimization and topology optimization.

The expressions derived for the effective dispersion using theory of
homogenization gives reasonable argument with experimental results.
Figure (3.2) and (3.3) shows that for low Peclet number the computed
values obtained from solving equation 3.11 matches closely with the
experimental results.

Further work was done to find the optimal microstructure using shape
optimization that maximizes the effective properties such as permeability

and dispersivity. The magnitude of the effective permeability computed is

174

and dispersivity. The magnitude of the effective permeability computed is
comparable to the work done by Guest and Prévost (Guest 2007). It is
observed in Figure (4.8) that the microstructure forms a square array for
both the cases with circular and rectangular microstructure. Among the
three microstructures, the circular microstructure gives the maximum
effective permeability; but the diagonal terms of the permeability tensor
are considerably higher than the other microstructures. Overall, the
microstructure with rectangular solid region shows the best result. This
may be due to the fact that it can be considered as a flow in a channel.
One notable observation made from the results obtained from shape
optimization for maximizing dispersivity is from Figure (4.13): the angle
between the length and the height goes to a lower value for circular and
elliptic cross section. But for a rectangular cross section, the 0 value
approaches 90. This complements the previous argument that it forms a
channel flow; and hence, it requires a comparative smaller pressure
gradient to optimize the dispersion tensor, which is of the same magnitude
of the other micro structures.

In the subsequent chapters I focused on finding optimal microstructure
that will minimize the energy loss and maximize mixing. Figure (5.3) and
(5.4) exhibits the dependence of pressure gradient on the objective
function (minimize the energy loss and maximize mixing). When the
results obtained using shape optimization for minimizing dissipation

power and maximizing “dispersion power” where the pressure gradient is

175

a design variable, were compared with microstructure formed by square
array of same porosity it showed significant improvement in reducing
energy loss while fluid flow and increasing mixing as shown in Figure
(5.11).

Similarly, when optimized microstructures obtained from topology
optimization were compared with the square array of same porosity shows
the same trend. When the macroscopic flow is compared for optimized
microstructure and square array formed by circular solid region with same
porosity, a higher magnitude of velocity is obtained for the optimized
microstructure as shown in Figure (6.6 a).

The work done in the following chapter to find the optimized layout of
the microstructure that will maximize mixing gave a mixed idea of how
the optimized layout should look. In Figure (7.5) the microstructure layout
is such that the flow is deviated toward the vertical direction. In all the
results Figure (7.5) - (7.6) the inflow in at negative 45 degree to the
horizontal and the aim here was to distribute the solute in the horizontal
direction. In Figure (7.5) the microstructure diverts the flow as much as
possible to the horizontal direction. For a lower volume fraction the solid
materials form a wall around the boundary of zone 1 and zone 2 to restrict
the flow in vertical direction. For a higher velocity in horizontal direction
the longitudinal dispersion tensor will be higher and hence maximize the

solute flow in that direction.

176

All the above mentioned work was done for single phase fluid flow.
But in most of the application of porous media the fluid is not of single
phase. In the last two chapters we developed theory that governs the
multiphase fluid flow in a porous media and then solve those partial
differential equations to design porous media. The governing equation
developed using mixture model gives the approximate behavior of the
multiphase fluid flow in a porous media. These governing equations were
solved using finite difference method. After optimization microstructure
with elliptic microstructure shows the least energy loss. From Figure (8.8)
it can also be added that for microstructure with elliptic solid region the
velocity profile shows less variance as compared with other two
microstructures. When the macroscopic flow is compared for optimized
microstructure obtained from shape optimization with square array formed
by circular solid region with same porosity, a higher magnitude of
velocity is obtained for the optimized microstructure as shown in figure
(8.8). The methodology developed here can also be used to understand
the effect of multiphase flow with variable volume fraction. This also
gives an opportunity to study further in the field of multiphase flow to
develop time depended expression for effective properties.

Further future work should be directed towards solving more problems
of similar nature. In addition to the solved optimization problem, such
methodology can be further extended to modeling complex materials such

as poroelastic material for artificial bones or teeth. Also experiments with

177

prototype optimized material should be pursued to validate the results
obtained from this work. From simulation perspective, further algorithms
and techniques should be developed to solve for transient cases since for
multiphase fluid the transient example will be more realistic as compared
to the results shown in this work. Volume of fluid method should be
further investigated using theory of homogenization to derive governing
equations for different scales that will identify each phase separately at
each time step. Finally for better understanding of the results further work
should also be done in solving the same examples in three dimensional

domain.

178

APPENDIX A

179

APPENDIX A

Appendix A shows the Matlab program used to solve the shape

optimization problem to maximize isotropic permeability in Chapter 4.

function main

clear all

clc

This part of the code defines the number of variables, constraints and
other parameters associated with it.

fnamel = ['G&P'];

fptl = fopen( fnamel, 'wt');

volfrac = 0.3;

mcons = 3;%number of constraints

numdesvar = 4;%number of variables

xmin = [001,001, 0.01, 0.01]';%lower area limit on area
xmax = [1, l, l, l]';%upper area limit on area
ITERMAX = 100;%maximum number of iteration
STEPSIZE = 0.1 ;%stepsize

A = 2;

R = 9.4413-4;

T = pi/2;

xval = [1,0.5,0.5,1]';%initial value of area
iter=1;

converged = 0;

30:01;

xold2 = xval;

xoldl = xval;

low = xmin - 30 *(xmax - xmin);
upp = xmax + 50 *(xmax - xmin);

df0dx2 =zeros(numdesvar, 1);
dfdx2=zeros(mcons,numdesvar);
a0mma = l;

cmma = 1000*ones(mcons, 1 );
dmma = 0*ones(mcons,1);
amma = 0*ones(mcons,l);

obj3 = 0;
0ij = 0;
obj] = 0;

while iter < ITERMAX & converged == 0
alpha=zeros(numdesvar,l);
beta=zeros(numdesvar, l );

180

alpha = max(xmin, (1 - STEPSIZE)*xval);
beta = min(xmax, (l + STEPSIZE)*xval);
This command calls a function to compute the objective function and
its sensitivity with respect to the design variables
[Obj, dObjda, dObjdb, dObjdc, dObjdd, KK, k, error, fem] =
sensitivitycalculation(xval(1)*A,xval(2)*R, xval(3)*R, xval(4)*T);
f0val = Obj;
df0dx = [dObjda*A, dObjdb*R, dObjdc*R, dObjdd*T]';
This command calls a function to compute the constraints and its
sensitivity with respect to the design variables
[con, dcon] = constraint(xval(1)*A,xval(2)*R,xval(3)*R, xval(4)*T,...
volfrac,A,R,T);
fval 2 con;
dfdx = dcon;
This command calls the MMA function to compute the optimization
process
[xmma,ymma,zmma,lam,xsi,eta,mu,zet,s,low,upp] =
mmasub(mcons,numdesvar,iter,xval,alpha,beta,xold1,xold2,
f0val,ddex,df0dx2,fval,dfdx,dfdx2,low,upp,...
a0mma,amma,cmma,dmma);
disp(['i=',sprintf('%3.0f, iter) ,
' Objective function=', sprintf('%9.9f‘, f0val ),...
' kiso=', sprintf('%9.9f‘, k ),...
' error=', sprintf('%9.9f', error ),...
' Constraint=', sprintf('%9.3f‘, fval ),...
' Aspect Ratio=', sprintf('%9.4f‘, xval(l) ),...
' major a =', sprintf('%9.4f', xval(2) ),...
minor a =', sprintf('%9.4f', xval(3) ),...
' theta=', sprintf('%9.4f', xval(4)*90 )]);
fprintf(fptl, '%f %f %f %f %f %f %f %f %f \n', iter,
f0val, xval(l), xval(2), xval(3), xval(4), 1000*KK(1,1),...
1000*KK(1,2), 1000*KK(2,1), 1000*KK(2,2) );
obj = f0val;
converged] = abs(obj-obj l );
converged2 = abs(obj l ~obj2);
converged3 = abs(obj2-obj3);
if iter > (ITERMAX-6)
set(gcf,'outerposition',[10,400,400,400]);
fname2 = ['K11-',int2str(iter)];
postplot(fem,
'tridata',{ 'Kl 1','cont','intemal' } ,
'trimap','jet(1024)',
'refine',3);
saveas(gcf, fname2,'jpg')
set(gcf,'outerposition',[410,400,400,400]);

181

fname3 = ['K22-',int2str(iter)];
postplot(fem,
'tridata', { 'K22','cont','intemal' } ,
'trimap','jet(1024)',
'refine',3);
saveas(gcf, fname3,'jpg')
set(gcf,'outerposition',[10,10,400,400]);
fname4 = ['K12-',int23tr(iter)];
postplot(fem,
'tridata', { 'Kl2','cont','intemal' } ,
'trimap','jet(1024)',
'refine',3);
saveas(gcf, fname4,'jpg')
set(gcf,'outerposition’,[410,10,400,400]);
fname5 = ['K21-',int2str(iter)];
postplot(fem,
'tridata', { 'K2 1 ','cont','internal' } ,
'trimap','jet(1024)',
'refine',3);
saveas(gcf, fname5,'jpg')
end
if converged] <= 0.0001 & converged2 <= 0.0001 & converged3...
<= 0.0001 & fval <= 0
converged = 0;
else
converged = 0;
end
iter = iter + l;
xold2 = xoldl;
xoldl = xval;
xval = xmma ;

obj3 = obj2;

0ij = objl;

objl = obj;

end
fclose(fptl );
clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Function to compute the constraints and its sensitivity with respect to
the design variables

function [con, dcon] = constraint(AR, major_axis,...

minor_axis, theta, volfrac,A,R,T)
B = 0.002;
conl = l-(pi*major_axis*minor__axis)/((volfrac)*AR*B"2);

182

dconlda = (pi*major_axis*minor_axis)/((volfrac)*AR"2*B"2)*A;
dconldb = -(pi*minor_axis)/((volfrac)*AR*B"2)*R;
dconldc = -(pi*major_axis)/((volfrac)*AR*B"2)*R;
dconldd = 0;
con2 = 2*minor_axis*l.5/(B*sin(theta))-1 ;
dcon2da = 0;
dcon2db = 0;
dcon2dc = 3/(B*sin(theta))*R;
dcon2dd = -3*minor_axis*cos(theta)/(B*sin(theta)"2)*T;
con3 = 3*major_axis/(AR*B)- l;
dcon3da = —3*major_axis/(AR"2*B)*A;
dcon3db = 3/(AR*B)*R;
dcon3dc = 0;
dcon3dd = 0;
con = [conl ;con2;con3];
dcon = [ dconlda, dconldb, dconldc, dconldd; dcon2da, dcon2db,
dcon2dc, dcon2dd; dcon3da, dcon3db, dcon3dc, dcon3dd];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Function to compute the objective function and its sensitivity with
respect to the design variables
function [Obj, dObjda, dObjdb, dObjdc, dObjdd, K, k, error, fern] =...
sensitivitycalculation(AR, major_axis, minor_axis, theta)
[OB], K, k, error, fem] = permeabilityanddispersionv2(AR,...
major_axis, minor_axis, theta);
Obj = OBJ;
Obj__old = Obj;
AR_new = AR*1.01;
[OBJ] = perrneabilityanddispersionv2(AR_new, major_axis,...
minor_axis, theta);
OBJa= OBJ;
dObjda = (OBJa-Obj_old)/(AR_new-AR);
major_axis_new = major_axis* 1.01 ;
[OBJ] = perrneabilityanddispersionv2(AR, major_axis_new,...
minor_axis, theta);
OBJma= OB];
dObjdb = (OBJma-Obj_old)/(major_axis_new-major_axis);
minor_axis_new = minor_axis* l .01 ;
[OBJ] = perrneabilityanddispersionv2(AR, major_axis,...
minor_axis_new, theta);
OBJmi= OBJ;
dObjdc = (OBJmi-Obj_old)/(minor_axis_new-minor_axis);
theta_new = theta* 1.01;
[OBJ] = perrneabilityanddispersionv2(AR, major_axis,...
minor_axis, theta_new);

183

OBJth= OBJ;

dObjdd = (OBJth-Obj_old)/(theta_new-theta);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Function to solve the governing equations and compute the objective
functions. .
function [031, Kavg, k, error, fem] =

permeabilityanddispersionv2(AR, major_axis, minor_axis, theta);

flclear fern

clear vrsn

vrsn.name = 'COMSOL 3.2';

vrsn.ext = ";

vrsn.major = 0;

vrsn.build = 222;

vrsn.rcs = '$Name: 3';

vrsn.date = '$Date: 2005/09/01 18:02:30 $';

fem.version = vrsn;

%geometry

B=0.002;

L: AR*B;

B1 = B*cos(theta);

B2 = B*sin(theta);

cl={curve2([0,L],[0,0],[1,l])};

c2={curve2([0,Bl],[0,B2],[l,1])};

c3={curve2([Bl,(Bl+L)],[B2,B2],[l,1])};

c4={curve2([(B1+L),L],[B2,0],[1,1])};

g1l=geomcoerce('curve',cl);

g2 l =geomcoerce('curve',c2);

g3 l =geomcoerce('curve',c3);

g41=geomcoerce('curve',c4);

gl=geomcoerce('solid',{gl l, g2], g3], g4] });

Pll = (Bl+L)/2;

P22 = B2/2;

g2=ellip2(major_axis,minor_axis,'base','center',‘pos',[Pl 1,P22]);

g3=geomcomp( { g1 , g2 } ,'ns',{ 'Rl’,'El ' } ,'sf,‘R 1 -E1 ','edge','none');

g4=geomcomp({ g3 } ,'ns', { 'CO 1 ' } ,'sf‘,’COl ','edge','none');

% geomplot(g4)

clear s

s.objs={g4};

s.name={ 'CO2'};

s.tags={'g4'};

fem.draw=struct('s',s);

fem.geom=geomcsg(fem);

%meshing

fem.mesh=meshinit(fem,'Hmaxsub',[ l, 4e-5]);

184

%defining the partial differential equations
clear appl
appl.mode.class = 'FlPDEC';
appl.dim = {'Kll','K12','K21','K22','A1','A2','K11_t','K12_t',...
'K2l_t','K22__t','A1_t','A2_t' } ;
appl.gporder = 4;
appl.cporder = 2;
appl.assignsuffix = '_c';
clear pnt
pnt.constr = {0,[0;0;0;0;'A1';'A2'} };
pnt.ind =[1,2,l,1,l,1,1,1];
appl.pnt = pnt;
clear bnd
bnd.type = {'neu','dir' };
bnd.h = {1,{l;0;0;0;0;1;0;0;0;0;l;0;0;0;
0:1;0;0;0;0;O;0;0;0} };
bnd.ind = [1,1,l,1,2,2,2,2];
appl.bnd = bnd;
clear equ
CQU-be = {{{0;Ol.{0;0}.{0;0}.{0;0}.{0;0}.{0;0};
{0:0},{0:0}.{0:0}.{0;0}.{0;0}.{0;0};{0;0}.{0;
0}.{0;0}.{0;0}.{-1.0;0}.{0;0};{0;0}.{0;0}.{0;
0}.{0;0},{0;0}.{-l.0;0};{0;0}.{0;0}.{0;0}.{0;
0},{0;-1-0},{0;0};{0:0}.{0;0l.{0;0}.{0;0},{0;
0}.{0;-1.0} } };
equ.c = {{0,0,0,0,0,0;0,0,0,0,0,0;'-eta',0,0,
0,0,0;0,'-eta',0,0,0,0;0,0,'-eta',0,0,0;0,0,
0,'-eta',0,0} } ;
equal = {{i1:0}.{0;0}.{0;1}.{0;0l,{0;0}.{0;0};
{0;0l.{ 1:0}.{0;0}.{0;l }.{0;0},{0;0};{0;0}.{0;
0}.{0;0}.{0;0l.{0;0},{0;0};{0;0}.{0;0}.{0;0},
{0;0l,{0;0},{0;0};{0:0},{0;0}.{0;0}.{0;0},{0;
0}.{0;0};{0;0}.{0;0}.{0;0},{0;0}.{0;0}.{0;0} } };
equ.f = {{0;0;-l.0;0;0;-1.0} };
equ.ind = [l];
appl.cqu = equ;
fem.appl{ l} = appl;
% Shape functions
fem.shape = {'shlag(2,"K1 1")','shlag(2,"Kl 2")',...
'shlag(2,"K2l ")','shlag(2,"K22")',...
'shlag(2,"A l ")','shlag(2,"A2")'};
fem.border = l;
fem.units = 'SI';
% Subdomain settings
- clear equ

185

equ.be= {{{0;0,} {0;,0} {0;,0} {00} { l.;,00} {0;0};.
{0;,0} {0;,0} {0;,0} {0;,0} {00} { 1.;00}; {00} .
{00} {00} {00} {0;- 10,} {0;0}; {00} {00}
{0;0},{0;O}.{0;0};{0;0}.{0;0}.{0;0}.{0;0}.{0;
0},{0;0} l};

equ.c = {{' -eta';' -eta'; '-eta'; '-eta' ;0 ,;0}}

equ.da = 1;

6911-31 = l{{0;0},{0;0},{0;0},{0;Ol.{0;0}.i0;0};
{0;0}.{0;0}.{0;0}.{0;0},{0;0}.{0;0};{0;Ol.{0; ..
0}.{0;0}.{0;0},{0;0}.{0;0};{0;0}.{0;0}.{0;0}.
{0;0},{0;0},{0;0};{-l.0;0},{0;0},{0;-1.0},{0;
0}.{0;0},{0;0};{0;0}.{-1.0;0}.{0;0}.{0;-1.0}.
{0;0}.{0;0} } };

equ.f = {{-l.0;0;0;-1.0;0;0} };

equ.ind = [1];

equ.dim 2 {'K1 1','K12','K21','K22','A1','A2'};

equ.var = {'absKl lx_c','sqrt(Kl 1x"2+K1 ly"2)',
'absculx_c','sqrt(cu1x"2+cu1y"2)',
’absK12x_c','sqrt(Kl2x"2+K12y"2)',
'abscu2x_c','sqrt(cu2x"2+cu2y"2)',
'absK21x_c','sqrt(K21x"2+K21y"2)',
'abscu3x_c','sqrt(cu3x"2+cu3y"2)',
'absK22x_c','sqrt(K22x"2+K22y"2)',
'abscu4x_c','sqrt(cu4x"2+cu4y"2)',
'absAlx_c','sqrt(A1x"2+Al y"2)',
'abscu5x_c','sqrt(cu5x"2+cu5y"2)',
'absA2x_c','sqrt(A2x"2+A2y"2)',
'abscu6x_c',’sqrt(cu6x"2+cu6y"2)' } ;

equ.cxpr = {'U','-(Kl 1*P1+K12*P2)', . .
'V',’-(K21*P1+K22*P2)'};

fem.equ = equ;

% Global expressions

fem.expr = {'Pl',0,

'P2','0',
'eta','0.001'} ;

% Coupling variable elements

clear elemcpl

% Extrusion coupling variables

clear elem

elem.elem = 'elcplextr';

elem.g = {'1'};

src = cell(l,l);

clear bnd

' bnd-exp” {{{},{}.'K11'},{{},{},'A2'}.{'K22'.{},{}},{'A2',-~

186

{}.{}l.{'K12',{}, {l}.{'K21'.{},{l}.{'K11'.{},{l}.i{}.{},--.
'K12'}.{'A1'.{ },i } },ll },l },'K21'}.{ { }.{ }.'K22'},{{ },i },'A1'}};

bnd.map = {'0','0','0','0','0','0','0','0','0','0','0','0' };

bnd.ind = { {'1'},{ '2','4','5','6','7','8'},{'3'} };

srC{1} = {{ }.bnd,{ } };

elem.src = src;

geomdim = cell(l,l);

clear bnd

bnd-map = {{{1.'1'.il}.{{l.'1'.{l},{{l.{}.'2'}.{{},{}.'2'l.--.
{l}.{},'2'}. {{H },'2'}.{{}.{},'2'},{{}.'1'.{}}.{{}.{}.'2'},-~
{{I.'1'.{ } M l },'1'.{ } H { },'1',{ l } };

bnd.ind = {{'1','3','5','6','7','8'},{ '2'},{'4'} };

geomdim{1}={{},bnd.{}};

elem. geomdim = geomdim;

elem.var = {'pconstr9','pconstr8','pconstr6','pconstr2',...
'pconstr ','pconstrS','pconstr3','pconstr10','pconstr1',...
'pconstr] l','pconstr12','pconstr7'} ;

map = cell(l,2);

clear submap

submap.type = 'linear';

submap.sg = '1';

submap.sv = {'l','7'};

submap.dg = '1';

submap.dv = {'2','8'};

map{ 1} = submap;

clear submap

submap.type = 'linear';

submap.sg = '1';

submap.sv = {'7','8'};

submap.dg = '1';

submap.dv = {'l','2'};

map{2} = submap;

elem.map = map;

elemcpl{ l} = elem;

% Point constraint variables (used for periodic conditions)

clear elem

elem.elem = 'elpconstr';

elem.g = {'1'};

clear bnd

bnd.constr = { {'pconstr9-(Kl1)','pconstr8-(A2)','0','0','0',...

'0','0','pconstr10-(K12)','0','pconstrl 1-(K21)’,...

'pconstr] 2-(K22)','pconstr7—(Al )' } , { '0','0','pconstr6-(K22)',
'pconstr2-(A2)’,'pconstr4-(K12)','pconstr5-(K21)',...
'pconstr3-(Kl 1)','0','pconstr1-(A1)','0','0','0'} };

bnd.cpoints = {{'2','2','2','2','2','2','2','2','2','2','2','2'},...

187

{'2','2','2','2','2','2','2','2','2','2','2','2'} };

bnd.ind = {{'2'},{'4'}};

elem.geomdim = { { { },bnd,{ } } };

elemcpl{2} = elem;

fem.elemcpl = elemcpl;

% Solution form

fem.solform = 'general';

% Multiphysics

fem=multiphysics(fem,

'sdl'.[]);

% Extend mesh

fern.xmesh=meshextend(fem);

% Solve problem

fem.sol=femlin(fem,
'conjugate','on',
'solcomp',{'K21','Kl1','K22','A1','K12','A2'},
'outcomp',{ 'K21','Kl 1','K22','Al ','K12','A2'});

% Save current fem structure for restart purposes

fem0=fem;

% Integrate

K1 la = 0.001*postint(fem,'Kl1','dl',[l])/(L*B);

K12a = 0.001*postint(fem,'Kl2','dl',[1])/(L*B);

K21a = 0.001*postint(fem,'K21','dl',[l])/(L*B);

K22a = 0.001*postint(fem,'K22’,'dl',[1])/(L*B);

Kavg = [K] la K12a; K21a K22a];

k: (K1 1a+K22a)/2;

error=((Kl 1a-K22a)"2+((Kl2a+K21a)/2)"2)/k"2;

OBJ=-k+error;

188

APPENDIX B

189

APPENDIX 8

Appendix b shows the Matlab program used to solve the topology
optimization problem to minimize dissipation power in Chapter 5.

flclear fern

clc

clear all

warning off

flclear fem

fnamel = ['minimisedissipation’];

fptl = fopen( fnamel, 'wt');

Da = le-3;

Define the Geometry

g]=rect2(],1,'base','corner','pos‘,[0,0])+rect2(1.2,l.2,'base',...
'comer','pos',[-0. 1,0 1]);

s.objs={gl };

s.name={ 'Rl’};

s.tags={'gl'};

fem.draw=struct('s',s);

fem.geom=geomcsg(fem);

Define Mesh

fem.mesh = meshinit(fem,'Hmaxsub',[l, 0.02, 2, 0.02]);

meshplot(fem.mesh)

fem.const.alphamin = 0.001;

fem.const.alphamax = l/Da;

fem.const.q = 0.01;

fem.const.eta = 0.01;

Phi0 = 1;

Define the variables and the governing equations

fem.sdim = {'x' 'y'};

fem.dim = {'u' 'v' 'p’ 'rho'};

fem.shape = [2 2 1 l];

fem.bnd.shape = {[1:3]};

fem.bnd.type = 'neu';

fem.bnd.g = {0,{0.001;0.001;0;0} };

fem.bnd.ind = [2,2,1,l,l,l,1,1];

fem.equ.shape = {[1:3] {1:41};

fem.equ.init = { {0;0;0;0.7} };

190

 

 

fem.equ.f = {{'alpha*u-l';'alpha*v+1';'ux+vy';1 },{'alpha*u-l ';...
'alpha*v+l';'ux+vy'; l } };

fem.equ.ga = { { {'-p+2*eta*ux';'eta*(uy+vx)' } ; { 'eta*(uy+vx)';...
'-p+2*eta*vy' }; {0;0};{0;0} } };

fem.equ.ind = [1,2];

% Coupling variable elements

clear elemcpl

Enforcing periodic boundary condition

clear elem

elem.elem = 'elcplextr';

elem.g = {'1'};

src = cell(l,l);

clear bnd

bnd-eXPr= {{{}JV',{}}.{{}.'U',{}}.{'U',{ },{}}.{'V',{ },{}}};

bnd.map = {'0','0','0','0'};

bnd.ind = { {'1’},{ '2'},{ '3','4','5','6','7','8'} };

SfCl1}= “},bnd,{”;

elem.src = src;

geomdim = cell(l,l);

clear bnd

bnd-map = {ll },'1',{}}.{{}.'1'.{}},{{ }.{}.'2'}.{{l,{}.'2'}l;

bnd.ind = {{'1','2','4','5','6','7'},{'3'},{'8'} };

geomdim{1}={{}.bnd.{}};

elem. geomdim = geomdim;

elem.var = { 'pconstr4','pconstr3','pconstrl ','pconstr2' } ;

map = cell(l,2);

clear submap

submap.type = 'linear';

submap.sg = '1';

submap.sv = {'2','8'};

submap.dg = '1';

submap.dv = {'l','7'};

map{ 1} = submap;

clear submap

submap.type = 'linear';

submap.sg = 'l';

submap.sv = {'7','8'};

submap.dg = '1 ';

submap.dv = {'l','2'};

map{2} = submap;

elem.map = map;

elemcpl{ 1} = elem;

Point constraint variables (used for periodic conditions)

clear elem

elem.elem = 'elpconstr';

191

elem.g = {'1'};
clear bnd

bnd.constr = {{'pconstr4-(v)','pconstr3-(u)','0','0'},{'0','0','pconstr1-(u)',

'pconstr2-(v)'} };

bnd.cpoints = { {'2','2','2','2'},{'2','2','2','2'} };

bnd.ind = {{'3'},{'8'}};

elem.geomdim = {{{},bnd,{ } } };

elemcpl{2} = elem;

fem.elemcpl = elemcpl;

Define the objective function

fem.equ.expr = {'A' 'eta*(2*ux*ux+2*vy*vy+(uy+vx)*(uy+vx))'...
'alpha', { 0,'alphamin+(alphamax-alphamin)*q*( l -
rho)/(q+rho)' }, } ;

fem.bnd.expr = {'B' '0' };

fem=multiphysics(fem);

Define the adjoint problem

fem = femdiff(fem);

fem.xmesh = meshextend(fem);

fem.sol = asseminit(fem);

femadj = fem;

femadj.equ.ga = {{{'diff(A,ux)’ 'diff(A,uy)'} {'diff(A,vx)' 'diff(A,vy)'}

{’diff(A,px)’ 'diff(A,py)'} {'diff(A,rhox)' 'diff(A,rhoy)'} } };

femadj.equ.f = {{'diff(A,u)' 'diff(A,v)' 'diff(A,p)' 'diff(A,rho)'} };

femadj.bnd. g = {{'diff(B,u)' 'diff(B,v)' 'diff(B,p)' 'diff(B,rho)'} };

femadjxmesh = meshextend(femadj);

flngdof(fem);

i4 = find(asseminit(fem,'Init',{'rho' 1},'Out','U'));

L = assemble(fem,'Out',{ 'L' });

Vgamma = L(i4);

Vdomain = sum(Vgamma);

i123 = find(asseminit(fem,'Init',{'u' 1 'v' 1 'p' 1},'Out','U'));

a0 = 1;

a=0;
c=20;
d=0;
xmin=0.01;
xmax: l;

xold = fem.sol.u(i4);
xolder = xold;

low = 0;
upp = 1;
penal = 3;

for iter = 1:1000
fem.sol = femnlin(fem,'Solcomp',{'u' 'v' 'p'},'U',fem.sol.u);
[K N] = assemble(fem,'Out', { 'K' 'N' } ,'U',fem.sol.u);

192

[L M] = assemble(femadj,'0ut', { 'L' 'M' } ,'U',fem.sol.u);

femadj.sol = femlin('In',{ 'K’ K(i123,i123)' 'L' L(il23) 'M'...
zeros(size(M)) 'N' N(:,i123)});

rho = fem.sol.u(i4);

Phi = postint(fem,'A','Edim',2) + postint(fem,'B','Edim',1);

dPhidgamma = K(i123,i4)'*femadj.sol.u-L(i4);

x = rho;

f = Phi/PhiO;

g = rho'*Vgamma/Vdomain-0.5;

dfdx = dPhidgamma/PhiO; dgdx = Vgamma'N domain;

d2fdx2 = zeros(size(rho)); d2gdx2 = zeros(size(rho'));

nel = size(x); nely = nel(l); nelx = nel(2); rrnin = 0.1;

[dfdx] = check(nelx,nely,rmin,x,dfdx);

[xnew,y,z,lambda,ksi,eta,mu,zeta,s,low,upp] = mmasub( 1 ,length(rho),iter,

x,xmin,xmax,xold,xolder,f,dfdx,d2fdx2,g,dgdx,d2gdx2,low,upp,a0,a,c,d);
xolder = xold; xold = x;

rho = xnew;

if iter >= 100

break

end

u0 = fem.sol.u; u0(i4) = rho;

fem.sol = femsol(u0);

disp(sprintf('lter.:%3d Obj.: %8.4f Vol.: %6.3f Change: %6.3f,
iter,f,rho'*Vgamma,max(abs(xnew-xold))))

fprintf(fptl, '%f %f \n',iter,f);

set(gcf,'outerposition',[ 10,400,400,400] );

fnameS = ['topo-',int23tr(iter)];
postplot(fem,'tridata','rho','trimap','gray')

axis equal; shg; pause(0.1)

saveas(gcf, fname5,'jpg')

end

193

APPENDIX C

194

APPENDIX C

In the topology optimization work shown in Chapter 6 and 7 the
derivative of objective function with respect to the design variable is
obtained using the adjoint problem as described by Olesen et al.
According to their work the derivative of the objective function can be
expressed using the chain rule,

30m+ i 30%ng (AC.1)

d

—Ob'=

dpl J] 8p 9 an 3p
m

where Obj is the objective function, p is the design variable and u is the

velocity-pressure vector. The starting point of the finite element analysis

is to approximate the solution component u,- on a set of finite element

basis functions or also known as shape function {goi’n}

(AC2)

where “Ln are the expansion coefficients. Similarly, the design variable

field p is expressed as

p = Loni?“ (AC3)
n

For a homogenized incompressible Stokes problem the standard Taylor—

Hood element pair with quadratic velocity and pressure approximation is

195

commonly used. For the design variable linear Lagrange elements are
chosen.

The governing equations and the boundary conditions are discretized by
the Galerkin method as I

3 T (AC.4)
L.(U,p)— z N..A.=0
I j=1 11 J

Mi(U. p) = 0 (AC.5)
where Ui’ Ai and p are column vectors with the expansion coefficients

for the solution “in, the Lagrange multipliers yin, and the design variable
field ,0", respectively. The column vector L,- contains the projection of

Neumann boundary conditions onto (an. The partial integration is given

by

Li n = I (¢i nFi +V¢i n fig“)..- I ¢i nGidan (AC6)
9 Q 9 9 an I
m m
The column vector M,- enforces the Dirichlet constraint in the boundary

conditions

Mi,n =Ri (AC.7)

The matrix Ni]: —6Mi/6Uj describes the coupling to the Lagrange

multipliers in Neumann boundary condition. The sensitivity analysis as

shown in equation (AC.1) can be also written as

196

i[Obj]=—aObj+ )2 LOW -a—U- (AC8)
dp 8p i=1 8U 8p

Equation (AC.8) is computed using the standard adjoint method. From

3
equation (AC.4-AC.5) we have for any p that Li — Z NSAJ. =0and

i=1

Ml.(U,p)=0. Therefore also the derivative of those quantities with

respect to p is zero. Adding any multiple such as O, and A, to equation

(AC8) does not change the result. Equation (AC8) can be written as

aObj+ % BObjBELJr (AC9)

d
—Ob'=
dp[ J] 8p 1&1an 8p

Reorganizing the terms,

- 3 ,_ al.. ,_ aM. (AC.10)
—d—[Obj]=——aObJ+ z UTi—'——A.T—‘ +
dp 3;) i=1 81) l 3p

U

'_-

i:

3 ' - 3 ~ aL. ~ 36.
)3 a—O—bl+ z UTi—i—A.TN.. -—'—
an wt. 1 ap

 

3
Z 2 U iN
i=lj=1 '1

The derivatives (BU/63p and aft/0p of the implicit functions can be

eliminated by choosing O, and A, such that

197

.. j: 3)ij (AC.11)

(AC. 12)

where Kij = —6L,-/6Uj . This problem is the adjoint of equation (AC.4) and

0i and Al. are the corresponding Lagrange multipliers.

198

REFERENCES

199

REFERENCES

Allaire, G (2002). Shape Optimization by the Homogenization Method.
New York: Springer-Verlag.

Allen, M. III (1985). Numerical Modeling of Multiphase Flow in Porous
Media. Advmes in Water Resources: 8, 162-187.

Amaziane, B, Bourgeat, A. and Jurak, M. (2005). Effective Macro
diffusion in Solute Transport through Heterogeneous Porous Media.
Society for Industrial and Applied Mathematics- Multiscale Modeling &
Simulation: 5, 1, 184-204.

Auriault, J.-L. (1980). Dynamic Behavior of a Porous Medium Saturated

by a Newtonian Fluid. International Journal of Engineering Science: 18,
775-785.

Auriault, J .-L. and Lewandowska, J. (1993). Homogenization analysis of
Diffusion and Adsorption Macrotransport in Porous Media:
Macrotransport in absence of advection. Geotechnigue: 43, 457—469.

Auriault, J.-L. and Lewandowska, J. (1997 May). Effective Diffusion
Coefficient: From Homogenization to Experiment. Transport in Porous
Media: 27, 205-223.

Auriault, J.-L. (2002, August). Up scaling Heterogeneous Media by

Asymptotic Expansions. Journal of Engineering Mechanics: 128, 817-
822.

Auriault, J.-L, Geindreau, C. and Boutin, C. (2005 July). Filtration Law in
Porous Media with Poor Separation of Scales. Transport in Porous Media:
60, 1, 89-108.

Baltean, D, Levy, T. and Balint, S. (2003 April). Diffusion—Convection in
a Porous Medium with Impervious Inclusions at Low Flow Rates.

minsport in Porous Media: 51, 19-39.

Bandyopadhyay, D, Benard, A and Diaz, A (2008). Design of Porous
materials using topology optimization. Submitted to AIAA

200

Bandyopadhyay, D and Benard, A (2008). Optimization of effective
properties and solute transfer in porous media using shape optimization.

In preparation.

Bandyopadhyay, D and Benard, A (2007) Optimization of solute transfer
in porous media using homogenization. Proceeding of the Inverse
Problems, Design and Optimization Symposium: p 1-8

Bang, B. and Dag Lukkassen, N. (1999). Application of Homogenization
theory related to Stokes flow in porous media. Applications of
Mathematics: 4, 309-319.

Bendspe, M. and Kickuchi, N. (1998). Generating optimal Topologies in
structural Design using a homogenization Method. Computer Methods in

Applied Mechanics and Engineering: 71, 197-224.
Bendspe, MP. and Sigmund, O. (2003). Tomlogy thimization-Theog,

Method_s_ and Applications. Berlin: Springer.

Bensoussan, A. (1978). Asymptotic Analysis for Periodic Structures, New
York: North-Holland Pub. Co.

Beran MJ (1968). Statistical Continuum Theories, New York: Interscience
Publishers

Bird, R. B, Stewart, W. E, and Lightfoot E. N (2001). Transport
Phenomena, Wiley.

Blackwell, R. J, Rayne, J. R. and Terry, W. M. (1959). Factors
influencing the efficiency of miscible displacement. Transactions AIME:
216, 1-8.

Borrvall, T. and Petersson, J. (2003 January). Topology optimization of
fluids in stokes flow, International Journal for Numericg Methods in
Fluids: 41, 77-107.

Buyuktas, D. and Wallender, W. W. (2003 October). Dispersion in
spatially periodic porous media. Heat and Mass flansfer: 40, 261-270.

Chella, R, Lasseux, D and Quintard, M ( 1998). Multiphase,
Multicomponent Fluid Flow in Homogeneous and Heterogeneous Porous
Media. Oil & Gas Science and Technology - Revue De L’Institut Francais
Du Petrole: 53, 3, 335-346.

 

 

201

COMSOL version 3.2 and 3.3 Reference Manual. Comsol Inc.

Delgado, J. M. P. Q. (2005 September). A critical review of dispersion in
packed beds. Heat Mass Transfer: 42, 279—3 10.

Diaz, A. R. and Sigmund, O. (1995 August). Checkerboard patterns in
layout optimization. Structural and Multidisciplinary Optimization: 10, 1,
40-45.

Diaz, A. R. and Benard, A. (2001 July). On the discretization of problems
involving periodic planar tilings. Communications in Numerical Methods
in Engineering: 17, 543—549

Diaz, A. R. and Bénard, A. (2003 March). Designing materials with
prescribed elastic properties using polygonal cells. International Journal
for Numerical Methods in Engineering: 57, 301—314.

Diaz, A. R. and Bénard, A. (2003 September). Topology Optimization of
Heat-Resistant Structures. Proceedings of DETC’03 ASME 2003 Design
Engineering Technical Conferences: Chicago, IL.

Edwards, M. F. and Richardson, J. F. (1968). Gas dispersion in packed
beds. Chemical Engineering Science: 23(2), 109-123.

Edwards, D. A, Shapiro, M, Brenner, H. and Shapira, M. (1991).
Dispersion of Inert Solutes in Spatially Periodic, Two-Dimensional Model
Porous Media. Transport in Porous Media: 6, 337—358.

Eschenauer, H. A. and Olhoff, N. (2001). Topology optimization of

continuum structures: a review. Applied Mechanics Reviews: 54, 4, 331—
389.

Eshelby, JD. (1961) Elastic inclusions and inhomogeneities. Proggess in
Solid Mechanics, vol. II. North-Holland, Amsterdam, pp. 87—140.

Evgrafov, A. (2004). Approximation of topology optimization problems

using sizing optimization problems. PhD Dissertation. Goteborg
University, Sweden.

Evgrafov, A (2005 July). Topology optimization of slightly compressible
fluids. Journal of Applied Matherrinics and Meganics: 86, 1, 46 — 62.

Evgrafov, A (2005 October). The Limits of Porous Materials in the

Topology Optimization of Stokes Flows. Journal Applied Mathematics
and Optimization: 52, 3, 263—277

202

Fluent version 6.2. Reference Manual. Ansys Inc.

Gersborg-Hansen, A. (2003). Tomlogy optimization of incompressible
Newtonian flows Moderate Reynolds Numbers. MS Thesis, Technical

University of Denmark. .

Gersborg-Hansen, A, Sigmund, O. and Haber, R. B. (2005 September).
Topology optimization of channel flow problems, Structural and
Multidisciplinary Optimization: 30 181—192.

Gersborg-Hansen, A, Bendspe, M. P. and Sigmund, O. (2006 March)
Topology optimization of heat conduction problems using the finite

volume method. Structural and Multidisciplinary Optimization: 30 181—
192.

Guest, J. K. (2005). Design of Optimal Porous Material Structures for
Maximized Stiffness and Permeability using Topology Optimization and
Finite Element Methods. PhD Dissertation. Princeton University.

Guest, J. K. and Prevost, J. H. (2005 December) Topology optimization of
creeping fluid flows using a Darcy—Stokes finite element. International
Journal for Numerical Methods in Engineering: 66, 461-484.

Guest, J. K. and Prevost, J. H. (2006 March) Optimizing multifunctional
materials: design of microstructures for maximized stiffness and fluid
permeability. lntemational Journal of Solids and Structures: 43, 22-23,
7028-7047.

Guest, J. K., Prevost, J. H. (2007 January). Design of maximum
permeability material structures. Computer methods in applied mechanics
and engineering: 196, 1006-1017.

Gunn, D. J. and Pryce, C. (1969). Dispersion in packed beds. Trans. Inst.
Chem. Eng: 47, T34l-T350.

Ishii, M. and Hibiki, T. (2006) Thermo-Fluid Dyflmics of Two-Phfi
Flow. New York: Springer

Jensen, J. S. (2003). Phononic band gaps and vibrations in one- and two-

dimensional mass-spring structures. Journal of SouncL and Vibration: 266,
5, 1053 -1078.

203

Jensen, J. S. and Sigmund. O. (2004). Systematic design of photonic
crystal structures using topology optimization: low-loss waveguide bends.
Applied Physics Letters: 84, 12, 2022-2024.

Jorgens, A . (2003) Topology Optimization of heat resistant composite
structures, PhD Dissertation Institut fiir Kunststoffverarbeitung.

Kaviany, A. (1991). Principle of Heat Trapsfer in Porous Media, New
York, Springer - Verlag.

Kitandis, K, and Dykaar, B. B. (1997 January). Stokes flow in a slowly
varying two dimensional pore. I_r_ansport in Porous Media: 26, 89, 89-98.

Klarbring, A, Petersson, J, Torstenfelt, B. and Karlsson, M. (2003
August). Topology optimization of flow networks. Computer Methods in
Applied Mechanics and Engineering: 192, 35, 3909-3932.

Kleinstreuer, C. (2003). Two-Phase Flow: Theory and Applications. New
York: Taylor and Franscis.

Kroner, E. (1978) Self-consistent scheme and graded disorder in
polycrystal elasticity. Physics of Fluid: 8, 2261-2267.

Kroner, E. (1986) Statistical Modelling. Netherlands: Elsevier.

Kurochkina, E. P, Soboleva, O. N and Epov, M. I (2004). Numerical
modeling of two- phase liquid flow in fractural porous medium. Journal of
Mixing Science: 40, 5 474-481

Larsen, U. D, Sigmund, O. and Bouwstra, S. (1996). Design and
fabrication of compliant mechanisms and material structures with negative

Poisson's ratio. IEEE, International Workshop on Micro Electro
Mechanical Systems, MEMS-96.

Lee, C. K, Sun, C-C. and Mei, C. C. (1996 March). Computation of
permeability and dispersivities of solute or heat in periodic porous media.
International Journal of Heat_and Mass Transfer: 39,4, 661-676.

Lomen, D.O, Islas, A. L, Fan, X. and Warrick A.W. (1991 October). A
perturbation solution for nonlinear solute transport in porous media.
Ipansport in Porous Media: 5, 5-6

Man, W, Donev, A, Stillinger, F. H, Sullivan, M. T, Russel, W. B,
Heeger, D, Inati, S, Torquato, S and Chaikin, P. M (2005). Experiments

204

on Random Packings of Ellipsoids. Physical Review Letter: 94, 198001,
1-4

Martys, N. S, Torquato, S and Bentz, D. P (1994). Universal scaling of
fluid permeability for sphere packings. msigal Review E: 50, l, 403 -
408 »

Mei, C. C. and Auriault, J.-L. (1989 December). Mechanics of
heterogeneous porous media with several spatial scales. Proceedings of
the Royal Society of London. Series A. Mathematical and Physical
Sciences: 426(1871), 391-423.

Mei, C. C. and Auriault, J. L. (1991). The effect of Weak inertia of flow
through a porous medium. burn}; of Fluid Mechanics: 222, 647-663.

Mei, C. C. (1991). Heat Dispersion in Periodic Porous Media by
Homogenization Method. Multiphase Transport in Porous Media: 122,
11-16.

Mei, C. C. (1992). Method of Homogenization Applied to Dispersion in
Porous Media. mnsport in Porous Media: 9, 261-274.

Mei, C. C. (1995). Mathematical Analysis in Engineering. New York:
Cambridge University Press.

Olesen, L. H, Okkels, F. and Bruus, H. (2005 September). A high-level
programming-language Implementation of topology optimization applied
to steady-state Navier—Stokes flow. International Journal for Numerical
Methods in Engineering: 65, 975-1001.

Pan, Y. and Home, R. N. (2001 October). Generalized Macroscopic
Models for Fluid Flow in Deformable Porous Media 1: Theories.
T_ransport in Porous Media: 45, 1, 1-27.

Pan, Y. and Home, R. N. (2001 November). Generalized Macroscopic
Models for Fluid Flow in Deformable Porous Media H: Numerical Results
and applications. unsport in Porous Mega: 45, 2, 195-213.

Pfannkuch, H. O. (1963). Contribution a l'étude des déplacements de
fluides miscibles dans un milieu poreux." Rev. I Fr. Petrol: XVIII(2),
170—215.

Pharm, D. C and Torquato, S. (2004 December). Exactly realizable

bounds on the trapping constant and permeability of porous media.
Journal of Applied Physics: 97, 013535, 1-8

205

Quintard, M and Whitaker, S (1988). Two-phase flow in heterogeneous
porous media: The method of large-scale averaging. lransport in Porous
Media: 3, 357-413

Quintard, M and Whitaker, S (1990). Two-phase flow in heterogeneous
porous media I: The influence of large spatial and temporal gradients.
Ignsport in Porous Mga: 5, 341-379

Rifai, M. N. E, Kaufman, W. J., and Todd, D. N. (1956). Dispersion
Phenomena in Laminar Flow through Porous Media. University of
glifomiaaBerkelei Berkeley, California.

Rubinstein, J. and Torquato, S. (1988). Diffusion-controlled reaction:
Mathematical formulation, variational principles and rigorous bounds.
Journal of Chemical Physics: 88, 10, 6372-6380.

Rubinstein, J. and Torquato, S. (1989). Flow in Random Porous Media:
Mathematical Formulation, Variational Principles, and Rigorous Bounds.
Journal of Fluid Mechanics: 206, 25-46.

Saez, A. E, Otero, C . J and Rusinek, I. (1989). The effective
Homogeneous Behaviour of Heterogeneous Porous Media. Transmrt in
Porous Media: 4, 213-238

 

Sanchez-Palencia, E. (1980) Non-Homogeneous Media and Vibration
Theory. Berlin: Springer. 1980.

Sigmund, O. (1995 June). Tailoring materials with prescribed elastic
properties. Mechanics of Materials: 20, 351-368.

Sigmund, O. and Torquato, S. (1996). Composites with extremal thermal
expansion coefficients. Applied Physics Letters: 69, 21, 3203-3205.

Sigmund, O. (2000). A new class of extremal composites. Journal of the
Mechanics and Physics of Solids: 48, 397-428.

Sigmund, O. (2001). A 99 line topology optimization code written in
Matlab. Structural and Multidisciplinary Optimization: 21 120—127.

Svanberg, K. (1987). The method of Moving Asymptotes-A new method
for structural optimization. International Journal for Numerical Methods
in Engineering: 24, 359-373.

206

Svanberg, K. (2002 January). A Class Of Globally Convergent
Optimization Methods Based On Conservative Convex Separable
Approximations. Journal on Optimization: 12, 2, 555—573.

Torquato, S and Pharn, D. C. (2004). Optimal Bounds on the Trapping
Constant and Permeability of Porous Media. Physical Review Letter: 92,
255505, 1-4

Wang, Y. (2003). A study on Microstructures of Homogenization for
Topology Optimization. PhD Dissertation. Victoria University of
Technology

Zijl, W and Trykozko, A (2002). Numerical homogenization of two-phase
flow in porous media. Computational Geosciences: 6, 49-71.

N .N Industrial Quick Search URL:
http://www.heatexchangcrs.org/info/heatexchangers, February 2007

N.N: H-IP—O URL: www.hipo-online.de, February 2007

N .N : Filtration and Separation Buyers Guide URL:
www.filtsepbuyersguide.com, February 2007

N.N: Cleveland Eastern Mixers URL: www.clcvelandmixcrcom,
February 2007

N .N: Corn Burner URL: www.thccornburncr.com, February 2007

207

   

u(ilijjjijijjjuiji1)u