 

 

 

 

:_ 7 x .v.
.1. 5: a: :5.» 5

This is to certify that the
thesis entitled

DESIGN AND ANALYSIS OF MESO-SCALE COOLING AND
STATIC MIXING DEVICES

presented by

RAMKUMAR SUBRAMANIAN

has been accepted towards fulﬁllment
of the requirements for the

 

 

 

 

 

Master of degree in Mechanical Engineering
Science
3 7 Major Professor’s Sig?ature
57/1/4954
Date

MSU is an Affirmative Action/Equal Opportunity Institution

 

 

LIBRARY
Michigan State
University

 

 

 

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

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 c:/ClRC/DateDue.p65-p. 15

DESIGN AND ANALYSIS OF MESO-SCALE COOLING AND
STATIC MIXING DEVICES

BY

RAMKUMAR SUBRAMANIAN

A THESIS
Subnﬁﬂedto
Michigan State University
in partial fulﬁllment of the requirements
for the degree of
MASTER OF SCIENCE
Department of Mechanical Engineering

2004

ABSTRACT

DESIGN AND ANALYSIS OF MESO-SCALE COOLING AND
STATIC MIXING DEVICES

By

Ramkumar Subramanian

In recent years, there has been a growing emphasis in manufacturing research on the
design and fabrication of Microtechnology-based Energy and Chemical Systems
(MECS). MECS are mesoscopic devices, which use embedded microstructures to
enhance heat & mass transfer and chemical reactions. Mesomachines are expected to
provide a number of important functions where a premium is placed on mobility,
compactness, or point application. The first part of this thesis investigates the design
aspects for fabricating a two-ﬂuid counter flow meso—scale heat exchanger. For this
investigation, a finite element analysis was carried out using Abaqus CAE to determine
the limits on ﬁn aspect ratios in two-ﬂuid microchannel arrays produced by diffusion
bonding. The second part of this thesis addresses the issue of using Functionally Gradient
Materials (FGM), achieved by spatially varying the compositional ratio between two
powders, to fabricate an active cooled and efficient meso-scale heat exchanger. Coupled
fluid-thermaI-structural finite element analysis was done using Ansys Multiphysics 7.0 to
minimize residual thermal stresses in the heat exchanger. The final part of this thesis
deals with the design of ribbed ducts to facilitate mixing of two liquids in microfluidic
devices. Fluent 6.0, a computational fluid dynamics package, was used to determine the
extent of mixing in the ribbed ducts for various combinations of fluid properties and inlet

conditions.

Man is made by his belief. As he believes, so he is.
Bhagavad-Gita

iii

To my loving parents and sister

ACKNOWLEDGMENTS

At the end of my thesis I would like to thank all those people who made this thesis
possible and an enjoyable experience for me. First of all I wish to express my sincere
gratitude to Dr. Patrick Kwon, who guided this work and helped whenever I was in need.
I think his presence at MSU was the best thing that could have happened to me as it
allowed me to hone my skills in multi-disciplinary fields of mechanical engineering. I
also wish to thank Dr. Kwon for allowing me to go on a semester long internship during

which I gained valuable experience in the field.

I am grateful to the members of my thesis committee, Dr. Dahsin Liu and Dr. Hyungson
Ki, for their valuable time and support. I also wish to express my deepest gratitude to Dr.
Tom Shih who helped me immensely with grid generation of the ribbed duct by allowing
me to modify his code. I acknowledge the financial support extended to me by the
Michigan State University for three semesters in the form of a Teaching Assistantship. I
am also grateful to Dr. Kwon for giving me the opportunity to work as a Research

Assistant under him in the summer semester.

Finally, I would like to express my deepest gratitude for the constant support,
understanding and love that I received from my parents and sister Ramya during the past
years. I am also thankful to my friends Neeraj Kohli, Vaibhav Ekbote, Nimisha Patel,
Sivom Manchiraju, Praveen Halepatali, Abhyudai Singh, Parmeet Batra, Monojit

Bhattacharya and Mohammad Hyderuddin for their constant encouragement and support.

TABLE OF CONTENTS

LIST OF TABLES viii
LIST OF FIGURES ix
1. INTRODUCTION 1

1.1 LIMIT ON ASPECT RATIO IN TWO-FLUID MICRO-SCALE HEAT EXCHANGER 2

1.2 DESIGN AN ACTIVE COOLED AND EFFICIENT MESO-SCALE HEAT EXCHANGER USING
A FUNCTIONALLY GRADIENT MATERIAL 7

1.3 DESIGN OF RIBBED DUCTS TO FACILITATE MIXING OF TWO LIQUIDS IN MICRO-

FLUIDIC DEVICES 12
2. FINITE ELEMENT ANALYSIS OF THE MICRO-LAMINATED HEAT EXCHANGER 18
2.1 PROBLEM DEFINITION 18
2.2 FINITE ELEMENT MODEL 21
2.3 EXPERIMENTAL APPROACH 27
2.4 RESULTS 28
2.5 DISCUSSION 30
2.6 CONCLUSIONS 32

3. FINITE ELEMENT ANALYSIS OF A OF FUNCTIONALLY GRADIENT MESO-SCALE

HEAT EXCHANGER 33
3.1 THEORETICAL BACKGROUND 34
3.2 FUNCTION ALLY GRADIENT MATERIALS 35

3.3 MORI-TAN AKA SCHEME 36

vi

3.4 DESCRIPTION OF DEVICE

3.5 FINITE ELEMENT MODEL

3.6 RESULTS AND DISCUSSION

3.7 VALIDATION

3.8 CONCLUSIONS

3.9 RECOMMENDATIONS FOR FUTURE WORK

4. COMPUTATIONAL STUDY OF MIXING MICROCHANNEL FLOWS

4.1 CHANNEL DESIGN AND SPECIFICATIONS

4.2 MESH GENERATION

4.2 MODELLING DETAILS, INITIAL AND BOUNDARY CONITIONS

4.3 SIMULATION RESULTS

4.4 ANALYSIS OF RESULTS

4.5 RESULTS AND DISCUSSION

4.6 FABRICATION OF MICROCHANNELS

4.7 FUTURE WORK

4.8 CONCLUSIONS AND EXTENSIONS

APPENDIX A

APPENDIX B

APPENDIX C

REFERENCES

vii

62

63

64

65

77

79

8O

80

81

89

90

91

LIST OF TABLES

Table l: The Geometrical Specifications of All Test Coupons and the Results of FE
Simulations for Average Critical Stresses

Table 2. Material property values of alumina and yttria-partially stabilized zirconia

Table 3: Comparing nodal temperatures from the hand calculated FD model with the
Ansys FE model

Table 4: Simulation results tabulated for the ribbed duct

Table 5: Simulation results tabulated for the ribless duct

viii

29

41

61

76

77

LIST OF FIGURES

(Images in this thesis are presented in color)

Figure l: Micro-lamination scheme used to fabricate a dual micro-channel array with
arrows indicating direction of flow 4

Figure 2: Exploded view of the FEA model of counter-ﬂow micro—channel array
investigated in this study 6

Figure 3: Photo of the Alumina-Zirconia FGM — one with internal channels and the other

with external channels 8
Figure 4: How FGMs work 10
Figure 5: The heat exchanger with cooling channels (side view) 11
Figure 6: Front view of the heat exchanger showing 8 cooling channels 11
Figure 7: Photo of the ribbed duct 15
Figure 8: Cut-out photo of the internal ribbed ducts 15
Figure 9: Photo of internal ribbed ducts 16
Figure 10: CFD model of the ribbed duct 16

Figure 11: a) The plan view of the two-ﬂuid interleaved, counter-ﬂow micro-channel
array comprised of micro-channel and fin laminae; and b) Cross-section of the device at
cross-section A-A 20

Figure 12: Cross-section of a diffusion-bonded counter-flow micro-channel heat
exchanger showing sections through which the bonding pressure: a) was transmitted; and
b) was not transmitted. Cross-section b resulted in leakage between the two ﬂuids

(Courtesy Dr. Brian Paul, Oregon State University) 21
Figure 13 (a) Geometry of the test coupon studied in this paper and verified by 23
Figure 14: S33 stress contours on test coupon with h = 152.4 um and a =l.3mm 25
Figure 15: S33 stress contours on test coupon with h = 152.4 ,um and a =l.7mm 25

Figure 16: (a) S33 stress contours from FEA on counterflow microchannel model (h =
152.4 mm and a =l.7mm). (b) Zoomed in view 26

Figure 17: Micrographs of test specimens consisting of 508 pm thick laminae: A) 5.08

mm SPAN (50X); B) 3.175 mm SPAN (50X), AND C) 2.54 mm SPAN (50X). (Courtesy
Dr. Brian Paul, Oregon State University)

Figure 18:

A plot of the experimental results showing the conditions under which

bonding and leakage occurred. Bonding conditions for all of these samples were

28

31

Figure 19: A plot of the experimental results showing lamina thickness Vs aspect ratio 31

Figure 20:

Figure 21

Figure 22:

Figure 23:
Figure 24:
Figure 25:
Figure 26:
Figure 27:
Figure 28:
Figure 29:
Figure 30:
Figure 31:
Figure 32:
Figure 33:

Figure 34:

Different layer thickness ratios for a three-layered FGM model

: Different layer thickness ratios for a six-layered FGM model

Temperature contours of the 8 -channeled alumina model

Maximum principal stress contours on the 8 -channeled alumina model
Maximum principal stresses at x = 0, y= 0

Alumina model - maximum principal stress contours at x=0

Three equi-layered FGM model - maximum principal stress contours
Maximum principal stresses at x = 0, y = 3mm

Maximum principal stresses for the 3 layered FGM model
Maximum principal stresses for the 6 layered FGM model
Maximum principal stresses for the 3 layered FGM model
Maximum principal stresses for the 6 layered FGM model

The two dimensional model for validation

Meshing the symmetric model for finite-difference calculations

Finite element model with one cooling channel, used to compare results

obtained from the finite difference calculations

Figure 35:
Figure 36:

Figure 37:

Meshing of ribbed-duct
Velocity vectors of the ribbed duct for water-water mixing at Re 2 800

Contour plots of volume fraction of liquid 2 in water-water mixing

42

43

46

46

47

48

48

49

50

50

51

55

56

59

67

69

70

Figure 38: Contour plots of volume fraction of liquid 2 in water-water mixing 71
Figure 39: Contour plots of volume fraction of liquid 2 in water-water mixing 71

Figure 40: Contour plots of volume fraction of liquid 100 times more viscous than water
when mixed with water in a ribbed duct with inlet velocity = 0.5m/s 72

Figure 41: Contour plots of volume fraction of liquid 300 times more viscous than water
when mixed with water in a ribbed duct with inlet velocity = 0.5m/s 72

Figure 42: Contour plots of volume fraction of water—glycerin mixing in a ribbed duct
with inlet velocity = 0.5m/s (top view) 73

Figure 43: Velocity vectors in a cross-section in the ribbed area 79

xi

1. INTRODUCTION

Over the last ten years, there has been a growing interest in the use of low—cost,
engineering materials to fabricate Micro Energy and Chemical Systems (MECS). MECS
are multi-scale fluidic devices, which rely on embedded micro-scale features to process
bulk ﬂuids for portable and distributed energy and chemical systems applications. MECS
are expected to provide a number of important functions where a premium is placed on
mobility, compactness, or point application. Examples include miniature refrigerators for
high-speed electronics portable power packs based on combustion [l] and miniature

chemical reactors for on-site waste processing [2].

One advantage of MECS devices over same-sized conventional devices is that they have
greater surface area available for heat transfer, reaction, separations, etc. The overall size
of MECS devices places them in the mesoscopic regime, i.e. in a size range between
macro objects such as automobile engines and laboratory vacuum pumps, and the
intricate MEMS based sensors that reside on a silicon chip. Thus MECS, although
having microfeatures, are large by MEMS standards straddling the size range between
the macro- and micro— worlds. This large surface area to volume ratio accounts for the

high rates of heat and mass transfer within MECS devices.

The internal processes of these devices rely on length scales that are much smaller than
traditional systems. For thermal and chemical applications, a small characteristic size
provides the benefits of high rates of heat and mass transfer, large surface-to-volume

ratios, and the opportunity of operating at elevated pressures. For other more

mechanically operated meso machines such as generators and motors, small dimensions
imply rapid response and compact design. Furthermore, these systems can often be

volume produced resulting in substantial cost reduction of each device.

In the energy area, MECS will find increasingly important uses where small scale heat
engines, heat pumps and refrigerators are needed. For example, the development of
miniature refrigerators could provide point cooling of high speed electronics and
communication equipment for enhancing performance. Also, power packs based on
combustion rather than electrochemistry could extend operating times of electronic
devices by a factor of ten. In the area of chemical processing, miniaturized chemical
reactors could provide on-site neutralization of toxic chemicals thereby eliminating the

need for transport and burial.

Because many MECS devices rely on fluidic processes, the same technology can be
applied to biological applications. Miniatun‘zed bioreactors could provide precisely
regulated environments for small groups of cells to enhance their production of
therapeutic drugs, or the detection of toxic compounds. Such bio-applications could
range from benchtop research to large scale production facilities. This thesis is an

attempt to design MECS devices which can be used for a variety of purposes.

1.1 LIMIT ON ASPECT RATIO IN TWO-FLUID MICRO-SCALE HEAT
EXCHANGER

The first part of the thesis focuses on the theoretical limits on fin aspect ratios within two-

ﬂuid micro-channel arrays produced by diffusion bonding. The counter-ﬂow device is

comprised of alternating layers of micro-channels, which allow the two fluids to flow in
opposite directions separated, by fins. The objective was to interpret the span of the fin as

a function of the fin thickness with the help of a finite element model.

Finite element analysis was performed on a model which assumes the micro-lamination
fabrication method. Micro-lamination involves the patterning and bonding of thin layers
of material called laminae [3]. Laminae could be metal shims or polymer ﬁlms having
the desired mechanical, thermal and chemical properties important to the device function.
The process begins by surface machining, or through cutting of a single laminate in the
desired pattern. Once the pattern is cut, the laminates are surface treated to enhance the
diffusion bonding and stacked in a prearranged order. The stack is then bonded together
forming a single block of material. For the method to have utility, a machining method
capable of fabricating structures in the laminating material is needed. The method must
be versatile, easy to use, and capable of rapidly machining (with through-cuts and surface
texturing) a wide variety of materials. Laser numerically controlled micromachining
satisfies these requirements. Other techniques such as through-mask electrochemical

machining and stamping are applicable to laminate forming.

Micro-lamination techniques have been used to fabricate MECS devices for advanced
climate control [4], solvent separation [5], micro combustion [6], fuel processing [7],
ﬂuid compression [8], microdialysis [9] and dechlorination among others. Several
research institutions and industrial companies have developed micro-lamination

approaches for producing MECS devices [3, 10]. In all, systems of micro-channel arrays

have been constructed in a wide array of materials including copper, stainless steel,

titanium, nickel, intermetallics and various polymers with features as small as five um.

 

 

 

 

 

 

 

 

 

 

 

.3,
.__-__l__
End Cap [Egljlwl
E D i—fD‘L‘i
Header and Channel ii ﬂ
Header and Fin __ [:1 D J]
a .1—
-.. I? [I H

 

 

 

 

 

Cl

 

 

Ir

 

 

\,\

End Cap ®\\ ®\\

 

 

 

 

Figure 1: Micro-lamination scheme used to fabricate a dual micro-channel array
with arrows indicating direction of ﬂow

In the case of heat transfer, the function of a typical micro-channel heat exchanger is to
transmit the heat from one fluid stream into a second ﬂuid stream. Several different
approaches can be used to accomplish this task. The most efﬁcient heat transfer methods
tend to interleave the two ﬂuids in an alternating succession of cross-ﬂow or counter-
ﬂow channels. In particular, counter-ﬂow channels are known to provide better
temperature distributions along the length of heat exchangers when compared to co-ﬂow
heat exchangers [11]. Whereas, in the case of single ﬂuid ﬂow, micro-channels with

relatively high aspect ratio’s have been fabricated successfully even in hard-to-produce

materials [12], the aspect ratio is generally constrained in two-ﬂuid micro-scale heat
exchangers because of the more complex geometry. This research investigates the
theoretical leakage limit on fin aspect ratios within two-ﬂuid MECS devices with the help

of finite element analysis.

High aspect ratio (AR) micro-channels can have very high surface area densities in
excess of 5,000 mZ/m3 where AR is defined as the ratio of channel width to channel
height. System miniaturization is provided by the reduction in residence time (and
consequent channel length) needed by ﬂuid molecules ﬂowing through micro-channel
reactors and heat exchangers. Brooks [13], Wang [14], Peng [15], Friedrich and Kang
[16] used the term surface area density to quantify this attribute of MECS devices defined

as the amount of surface area per unit volume of device.

Further, the high AR micro-channels within MECS devices tend to have low pressure
drops when compared with other high surface area density structures (e.g. sponges). In
fact, the characteristic reduction in channel length causes MECS devices to have no
theoretical increase in pressure drop over conventional systems. Overall, MECS devices
provide the opportunity to greatly reduce the size and weight of energy and chemical
system applications while providing the opportunity for enhanced process control for

novel material synthesis.

Results from metallography and leakage testing performed by Dr. Brian Paul at Oregon

State University revealed that for a given set of bonding conditions there exists a

 

 

 

 

maximum permissible value of aspect ratio (width/height of the cross-section of a
channel) during bonding beyond which leakage is bound to occur. Furthermore, the
results implied that higher aspect ratios can be achieved at smaller scales. Hence, the
FEA model was used to determine these limiting values of aspect ratios for various

laminae using Abaqus CAE to simulate the counter-ﬂow device.

 

Figure 2: Exploded view of the FEA model of counter-ﬂow micro-channel array
investigated in this study

1.2 DESIGN AN ACTIVE COOLED AND EFFICIENT MESO-SCALE HEAT
EXCHANGER USING A F UNCTIONALLY GRADIENT MATERIAL

Many researchers have attempted to remove heat from microelectronic devices. One of
the more popular methods to dissipate the heat is to spread rapidly by employing highly
conductive and high heat capacity materials such as diamond, silicon nitride,
molybdenum, etc. [21, 22, 23]. Heat is dissipated to the environment either by forcing air
through pin arrays or fins or by cooling naturally. Recently MEMS-sized heat exchangers
have emerged which rely on convective heat transfer via air. In such cases, heat

inevitably spreads to its surrounding and other peripheral devices.

More efficient cooling results when coolants circulated through networked channels and
manifolds in a closed system. Bower and Mudawar [24] and Gillot [25] developed a two-
phase heat exchanger with a higher efficiency. Bower and Mudawar [24] made mini and
micro-channels on blocks made of copper and nickel plates. This work is an attempt to
design an efficient active cooled heat exchanger by minimizing residual thermal stresses
in the heat exchanger. To achieve this, we propose to the body of the heat exchanger to

be micro-textured.

Micro-textured or Functionally Gradient Material (FGM) can be processed by spatially
modulating the compositional ratio between two powders. Further, the body of the heat
exchanger consists of a network of micro-scale channels and manifolds through which
the coolant is transported. Applications for the heat exchanger include a variety of other
industrial or medical uses, such as conduits for ﬂuids in engines or to carry medicines or

biological ﬂuids.

 

Figure 3: Photo of the Alumina-Zirconia FGM - one with internal channels and the
other with external channels

Structural and thermal properties of FGMs can be tailored by spatially varying the
composition of the second-phase particles in the matrix material. This is achieved by
modulating the concentration level of the second phase during fabrication. Therefore, the
FGM fabrication processes should possess the ability to produce various distributions of

particles to withstand assigned loading.

Powder metallurgy is the most versatile and most popular method of fabricating a bulk
FGM. The state of powder metallurgy technologies in fabricating FGM is reviewed by

Watanabe and Kawasaki (1990). This process involves sintering of powder mixtures with

varying volume fractions of the constituent materials. One drawback behind powder
metallurgy is that controlling the concentration levels of the constituent materials can be

quite cumbersome.

An FGM ideally has continuous gradiency. One can approximate such gradient properties
by preparing subcomponents in the form of layers. Each layer is prepared by
homogeneously mixing two powders in various ratios, which can be stacked to form the

‘discrete’ gradiency in thermo-mechanical properties.

One area of micromechanics is concerned with calculating the effective or bulk physical
properties of locally inhomogeneous, but statistically homogeneous, materials such as
composite materials. Composite materials consist of at least two distinct phases — the
matrix material and the reinforcing material, whose properties are known. Rather than
microscopically modeling each inclusion in a composite medium, a given

inhomogeneous material can be replaced by an equivalent “effective medium”.

For our finite element analysis, alumina and zirconia were chosen for their inherent
interfacial strength and difference in their coefﬁcient of thermal expansion (CTE) values.
The Mori-Tanaka (1973) approach was used to solve for the properties of the FGM as a

function of the volume fraction of the particles or fibers.

Because of the almost similar sizes of the alumina and zirconia particles, the Mori-

Tanaka mathematical model was solved twice-once for each of the two materials as the

matrix. It was found that the properties of the intermediate discrete layers were almost
similar by using either of the two models. Hence, the average values were taken as the

final material properties of the FGM.

In the micro—textured body of our heat exchanger with cooling channels, heat is

conducted and dissipated via coolant circulating through the manifold channels (Fig 5).

The novel processing techniques used to fabricate the heat exchanger involve:

(1) Processing techniques to fabricate internal channels into an FGM by using a
fugitive phase.

(2) Powder processing techniques to produce the micro-textured heat exchanger

undamaged by the process-induced stress.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

  

 

 

 

 

 

 

 

 

 

 

 

 

 

(a) (b) (C)
{D
0 0 ‘-
5 .‘é ii
2 S B
:- 0 D.
e e E
O 0 I—
l- :
l‘.’ “’ '-
o '5 0
m
«n m m
30 eo___ ___ er--- -__
‘7, U)
l
Homogeneous Uniform Temperature Thermal Gradient Thermal Gradient FGM
Material Field Loading Loading

 

 

 

 

Figure 4: How FGMs work

Th-——~

q),

__—_.

 

Q...

Figure 5: The heat exchanger with cooling channels (side view)

 

Figure 6: Front view of the heat exchanger showing 8 cooling channels

The key issue in producing micro-textured medium (FGMs) is to avoid damage due to
process-induced residual stress. Residual stresses in FGMs [26, 27] occur due to spatial
variations in Coefficient of Thermal Expansion (CTE) and densification behavior during
processing. However, it is also this variation in the CTE values which is one of our
design parameters used to minimize the thermal stresses because of a large thermal
gradient. Hence, with the help of FEA an attempt is made to design the heat exchanger
such that the effect of thermal gradient loading is offset by the variation of the CTE to

give us minimum thermal stress.

1.3 DESIGN OF RIBBED DUCTS TO F ACILITATE MIXING OF TWO
LIQUIDS IN MICRO-FLUIDIC DEVICES

Rapid mixing is essential in many of the micro-ﬂuidic systems targeted for use in the
testing of biochip analysis, drug delivery, and analysis or synthesis of blood, among
others. Mixing initiates many of these processes by rapidly increasing rates of enzyme
reactions, cell activation or protein folding. Effective mixing is characterized by an
increase in the contact area between the ﬂuids and a decrease in the distance over which
diffusion needs to act such that diffusion can complete the mixing process rapidly [28]. In
macroscopic devices, turbulent three dimensional ﬂow structures or mechanical actuators

can be used to accomplish this.

However in micro electro mechanical systems mixing of sample and reagent ﬂuids
quickly is a major issue since the physics of ﬂuids at these length scales precludes the
normal available mechanisms for mixing. Mixing is dominated by molecular diffusion

process in which the diffusion time is proportional to the square of the length scale.

12

However for a channel with dimensions of about 100 pm, the diffusion times are too slow

to rely upon for effective mixing [29].

The difﬁculty in mixing ﬂuids on the microscale lies in the size of the devices. Because
of the small scales, these ﬂows are predominantly laminar, so the efficient mixing
obtained in turbulent ﬂows is not practically attainable. While it is possible to design
active micromixers that utilize macroscale mixing techniques such as mechanical stining,
the inclusion of active elements often leads to very complex and costly fabrication
schemes. Turbulence occurs in ﬂows characterized by high Reynolds number Re > 2000.
For a characteristic velocity and length scale U, D and kinematic viscosityv, Re is

defined as: Re = UD/ v

For typical dimensions of MEMS devices around 200 um and assuming highest ﬂow
velocities to be around 100 mm/s, we have Re = 20. Typically Re is much lower in the l-
10 range and for creeping ﬂows Re < 1. This ﬂow regime severely precludes turbulence

as a possible mechanism.

For planar devices, diffusion is limited by the in-plane dimension of the ﬂuid channel.
Using Fick’s equation, the diffusion mixing time is deﬁned by, Ta = L2 / K , where L is
the characteristic mixing length, KiS the Fick’s diffusion coefficient. For salt in water, K
= 1000 um2 / s, and using L = 200 pm, we find Ta = 40 seconds. Thus mixing times by

diffusion are too slow for the viability of MEMS mixes in most applications [28].

13

Most designs of micromixers can be broadly classified as either active mixers or passive
mixers. Active mixers are those that exert some form of active control over the ﬂow ﬁeld.
They can usually do an excellent job of mixing. However they normally use
electrokinetic or thermally actuated ﬂow and actuators, valves and pumps which make
them difficult to fabricate, operate and integrate into microﬂuidic systems. Passive
micromixers utilize no energy input except for a pressure head to drive the ﬂow. The key
advantages that they offer are the simplicity in fabrication and operation and reduction in

COSI.

Partly due to the challenge of achieving high mixing rates in passive mixers, the literature
contains signiﬁcantly lesser number of designs of passive devices compared to active
ones. Our goal has been to design a passive micromixer that allows fast mixing of small
amounts of liquids. The ducts have inclined ribs or patterned grooves on all four walls

which assist the mixing process by the vortex phenomenon of ﬂuid stirring.

Most of the earlier computational studies on internal coolant passages have been two-
dimensional. In recent years a number of three-dimensional studies have been reported.
Three dimensional studies are needed on ducts with ribs, 180 degree bends, and/or
rotation. Very few investigators performed computational studies on duets with inclined
ribs, which are used in advanced designs. When computing duets with ribs, it is important
for the geometry of the ribs to be captured correctly. This is because they exert

considerable inﬂuence on the ﬂow and heat transfer.

14

 

Figure 7: Photo of the ribbed duct

 

Figure 8: Cut-out photo of the internal ribbed ducts

 

 
  

INIET 1

Figure 9: Photo of internal ribbed ducts

Figure 10: CFD model of the ribbed duct

The model was solved using Fluent 6.0, a commercial computational ﬂuid dynamics
code. Mixing conditions were simulated for a few combinations of liquids that are known
to be almost immiscible. The results of mixing between these liquids in the ribbed duct

were then compared to those obtained from a straight simple duct.

Shih, et al [3] studied the effect of a ribbed U-duct with or without rotation. Their model
consisted of a U-duct with inclined ribs on just two opposite faced walls. They developed
a mesh capturing the effects of the side ribs to perform computational ﬂuid dynamics
simulations. We used the code written by Shih et al [30] and modiﬁed it so that the mesh

now represented our geometry.

Our geometry was a straight duct which consisted of two pairs inclined ribs, both being at
an angle of 90 degrees to each other. Near the inlet the ribs were on the top and bottom
sides of the duct, while near the outlet the inclined ribs were placed on the two opposite

side walls. We came up with this design to ensure a more thorough mixing.

The code by Shih et a1 [30], which was in NASA plot3D format, was first converted to
ASCII dat format readable by Gambit a CFD meshing tool. Then a code in C++ was
written to modify the grid according to the required geometry. Later, the modiﬁed grid
was imported into Gambit where the boundary conditions were applied. The same mesh
was then used to simulate the ﬂow and mixing in Fluent 6.0 for a few combinations of

almost immiscible liquids.

l7

2. F INITE ELEMENT ANALYSIS OF THE MICRO-LAMINATED
HEAT EXCHANGER

This part of the thesis investigates the limits of fin aspect ratios within two-ﬂuid counter-
ﬂow micro-channel arrays based on the minimum required stress state needed to laminate
arrays by diffusion bonding. In the paper by Paul [20], it has been established that the
diffusion bonding of two-ﬂuid systems by micro-lamination can result in regions of the
device that do not directly transmit bonding pressure and, consequently, result in
unbonded regions leading to device leakage. A finite element model for analyzing the
stress state in these unbonded regions was used to determine the minimum stress state
necessary to adequately bond a given geometry. Model results were then compared with
experimental results for a range of different counter-ﬂow geometries. It has been found
that under a given set of bonding conditions, a minimum stress state exists which must be
met in every part of the geometry in order to produce leak-free bonds. Implications of
this ﬁnding on the design of two-ﬂuid micro-channel arrays were discussed. The
experimental work was carried out by Dr. Brian Paul at the Oregon state university. The
results of these ﬁndings were presented in Paul et al 2004 [39]. This chapter of the thesis

has been adapted from the above research paper.

2.1 PROBLEM DEFINITION

For this investigation, the case of a two-ﬂuid, counter-ﬂow, micro-channel array was
investigated though it is expected that the conclusions of this work will apply to other

types of two-ﬂuid micro-channel arrays. The counter-ﬂow device is comprised of

18

alternating layers of micro-channels, which allow the two ﬂuids to ﬂow in opposite
directions separated by fins. A concept of the lamina design for a micro-laminated two-
ﬂuid, interleaved, counter-ﬂow micro-channel array is shown in Figure 2. It is assumed
that once the laminae are stacked they are diffusion bonded according to Figure 2. The
plan view of the resulting monolithic device comprised of micro-channels and fins is

shown in Figure 11.

In order to diffusion bond the laminae, a uniform bonding pressure is applied on the stack
at an elevated temperature. The pressure is transmitted uniformly through the device
except at the neck of the micro-channel (where the cross-section AA is in Figure 11). The
lack of bonding pressure causes the laminae to remain unbonded at points A & B

resulting in leakage within the device and mixing of the two ﬂuids.

Figure 11 shows a region where pressure was directly transmitted through all of the
laminae and exhibits excellent bonding. Figure 12 shows a region where pressure was
not transferred directly between adjacent laminae because of the counter-ﬂow design

[20]. Figure 12 exhibits unbonded regions.

By experience, the width of the channel has to be constrained to a dimension where
bonding will occur. The limit on the width of the channel, at which bonding takes place,
depends on the thickness of the lamina. In other words, as the thickness of the ﬁn
laminae above and below the neck of the micro-channels increases, the neck can be made

wider.

19

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

III: ..................:::l
A (— '
C12: ................... 'El
A
(a)
P
/ / \/ / / / / .
Huidttl CHANNEL //7/] l// / 2
, 3
Fluid #2 (‘HANIEIIEE ”I,” 4
FIN / / 7//_5
Fluid #1 CHANNEL 6
\7
P
(b)

Figure 11: a) The plan view of the two-ﬂuid interleaved, counter-ﬂow micro-channel
array comprised of micro-channel and tin laminae; and b) Cross-section of the
device at cross-section A-A

20

-._..~‘._- ‘7.

$00.92*.

.~
l I

"l l l l l
[ I

 

a) b)
Figure 12: Cross-section of a diffusion-bonded counter-ﬂow micro-channel heat
exchanger showing sections through which the bonding pressure: a) was

transmitted; and b) was not transmitted. Cross-section b resulted in leakage
between the two fluids (Courtesy Dr. Brian Paul, Oregon State University)

The ratio of the channel span in the neck region to the fin lamina thickness is the fin
aspect ratio. It is expected that for a given set of bonding conditions (time, temperature,
bonding pressure), there will exist a critical minimum pressure that all areas of the
laminae must receive in order to bond. This minimum pressure will ultimately constrain

the fin aspect ratio for any given two-ﬂuid micro-channel device.

2.2 FINITE ELEMENT MODEL

An analytical approach to modeling this phenomenon has been attempted based on
simple elastic beam theory [20]. However, no generalizable conclusions were gained.
Here a finite element model is used to analyze the effect of changing lamina thicknesses
and channel (neck) widths.

To simplify the modeling of this phenomenon, a test coupon was developed (Figure 13)

to approximate the conditions of the necked region shown in Figures 11 (a) and (b). To

21

validate the geometry of the test coupon, FEA was performed on both geometries (Figure
2 and 13) to verify that the values of the out of plane stresses are comparable with those
in the test coupon model. To do this, a simple linear elastic model was found adequate
since it is recognized that the channels in the lamina become the stress concentration.
The two geometries were found be comparable across a wide range of channel widths and

fin thicknesses.

For the analysis, certain assumptions were made. The 6—layered test coupon used was
assumed to be 25.4 X 25.4 mm in width and breadth with varying thickness depending on
lamina thickness, h. Critical dimensions of the first and third laminae in the test coupon
are shown in Figure 13 (b) including the span (width) of the channel, a. As the model
was always under compression, to reduce contact non-linearity problems, the analysis
was performed on a 612 thick 3D deformable model instead of 6 separate layers each with
thickness h. Two rigid three-dimensional discrete plates were used on the top and bottom

to apply the bonding pressure.
The model was a linearly elastic isotropic model with the following material properties:

E = 124,800 N/mmz, u = 0.2235. These values were chosen as the respective elastic

properties at the temperature where diffusion bonding of the lamina took place.

22

   

6.34mm

(b)
Figure 13 (3) Geometry of the test coupon studied in this paper and verified by

leakage testing. (b) Specific critical dimensions of the test article used in this study

23

Contact was defined between the rigid plates and the set of the deformable plates
representing the heat exchanger. Load was applied in 3 steps to prevent numerical
singularity:
(1) The deformable body was constrained in the Z direction. Then the top &
bottom rigid plates were moved towards the body of the heat exchanger by a

very small distance (1 e-08 mm) to establish the contact between the nodes.

(2) After removing the Z direction constraint on the deformable body the bottom

rigid plate was constrained in all 6 degrees of freedom.

(3) After removing the displacement BC on the top rigid plate, a pressure of 5.86
MPa was applied on the top rigid plate. The load was applied in steps of 100

to prevent numerical singularity and negative eigenvalue.

An 8-noded brick element with reduced integration (C3D8R) was used for all the models.
For the rigid plates, 3D rigid elements were used. FEA was performed for the following
values of thickness of the lamina plates: 76.2 pm, 152.4 mm, 203.2 pm, 254 um, and 508
pm. For each of these values of thickness, a range of channel spans were chosen and FEA
was used to predict the limiting value of the channel span to avoid leakage. As shown in
figure 12, the leakage problem occurs in the neck region of the channels of the heat
exchanger model and on the third layer from the top of the test coupons. Thus, the critical
stress is deﬁned as the out-of-plane nodal stresses (833) on the “bridge region” of the
layer 3, i.e. the region on layer 3 between the channel and the slot on the right. The stress

contour images [Figure 14, 15, 16(a) and 16(B)] in this thesis are presented in color.

24

 

Figure 14: S33 stress contours on test coupon with h = 152.4 run and a =1.3mm

 

Figure 15: S33 stress contours on test coupon with h = 152.4 um and a =1.7mm

25

A}, -.
”M" -IWIW/l -irm A
_ " 'r u.
p x. ,‘l'lm Eh!!!”

.k ["71” Eb! \ p
I \,
‘ 94,. E4,

\ ‘Ar 1 i:
\\ “ M‘IIIIII .1

 

(b)

Figure 16: (a) S33 stress contours from FEA on counterﬂow nricrochannel model (h
= 152.4 um and a =1.7mm). (b) Zoomed in view

26

2.3 EXPERIMENTAL APPROACH

All experimental work was done by Dr. Brian Paul’s research group at Oregon State
University. The fabrication procedure for making the test specimens (Figure 13 (a))
consisted of a micro-lamination procedure involving laser machining and diffusion
bonding. First, laminae were patterned from blanks of 304 stainless steel shim stock of
various thicknesses (76.2, 152.4, 203.2, 254 and 508 am) using a 532 nm Nd: YAG laser
installed on a laser micromachining system. Overall size of the laminae was 25.4 x 25.4
mm. Critical dimensions of the ﬁrst and third laminae are shown in Figure 13 (b).
Efforts were made to remove laser burrs by abrasion after laser machining. Laminae
were cleaned with an acetone, methanol and DI water for degrease and rinse and put into
an alignment ﬁxture to be bonded in a vacuum hot press. The base of the hot press ﬁxture
was made of graphite with tungsten pins for facilitating the edge alignment of the
laminae. All laminae were diffusion bonded in the vacuum hot press at 900° C and 5.86

MPa for 2 hours.

Preliminarily, a set of test coupons were developed to verify the existence of the bonding
problem under study. Laminae thicknesses for preliminary coupons were 254 um and
508 pm. Metallography was conducted on these coupons to verify bonding conditions.
Once the process was veriﬁed, additional test specimen were produced in the remaining
lamina thicknesses. To speed up testing, the bonds of these devices were evaluated by
means of a leakage test in which the open header was pressurized by air to several

atmospheres and held under water to look for air leakage.

27

2.4 RESULTS

To begin the experimentation, three specimens of 508 pm thick SS with channel spans of
2.54 mm, 3.175 mm, and 5.08 mm were diffusion bonded at the conditions described
previously. Metallography of the specimens revealed that the specimens with 2.54 mm
and 3.175 mm spans were bonded below and above the channel, whereas small gaps
between laminae were detected for the specimen with the 5.08 mm span. Figure 17

represents the micrographs of the 508 um thick specimens.

 

 

 

 

 

 

A) B) C)

Figure 17: Micrographs of test specimens consisting of 508 um thick laminae: A)
5.08 mm SPAN (50X); B) 3.175 mm SPAN (50X), AND C) 2.54 mm SPAN (50X).
(Courtesy Dr. Brian Paul, Oregon State University)

Similarly, three specimens of 254 um device with a channel span of 1.524 mm, 2.54 mm,
and 5.08 mm were diffusion bonded. Metallography of the specimens showed that

specimens comprised of 5.08 mm and 2.54 mm were not bonded above and below the

channel span, whereas the specimen with 1.524 mm thick laminae did not have any gaps.

28

- ... [.154 Hill: It...“ .413:

 

Micrographs of these specimens are shown in Figure 17. Further experiments were
performed by fabricating and leakage testing specimens consisting of 76.2 um, 152.4 um,
and 203.2 pm thick laminae. A summary of results from the FEA model is shown in
Tablel. These results also compared well both with the metallography and leakage tests.
The values displayed in Table l are the averaged critical stress values that were obtained

after averaging all the nodal out-of—plane stress in this bridge region.

Table 1: The Geometrical Specifications of All Test Coupons and the Results of FE
Simulations for Average Critical Stresses

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

FE Lamina Channel Aspect Leakage/ Average Critical
simulation Thickness, h Span, a Ratio, a/h Bonded Stress (MPa)
Case (mm) (mm)
1 0.0762 0.75 9.84 Bonded - 0.17
2 0.1524 0.50 3.28 Bonded - 0.45
3 0.1524 1.30 8.53 Bonded - 0.05
4 0.1524 1.70 11.15 Leakage +0.43
5 0.1524 2.50 16.4 Leakage +0.44
6 0.2032 1.0 4.92 Bonded - 0.30
7 0.2032 1.5 7.38 Bonded - 0.07 to + 0.1
8 0.2032 2.2 10.83 Leakage + 0.45
9 0.2032 3.5 17.22 Leakage + 0.39
10 0.254 1.8 7.09 Bonded - 0.05 to + 0.15
11 0.254 2.6 10.24 Leakage + 0.36
12 0.508 3.1 6.1 Bonded - 0.026 to + 0.15
13 0.508 5.1 10.04 Leakage + 0.30

 

 

 

 

 

 

29

 

For most cases, the value of critical stress is almost constant throughout the bridge
region. However for a few cases (Cases 7, 10 and 12 in Table l), we see that this value is
not constant unlike other cases. There is a slight variation in the stresses showing some
nodes being in tension and other nodes being in compression. For example, the variation

in stress for Case 10 ranges from -0.05 to 0.15MPa.

The stress contours in Figures 16 (a) show the out-of-plane stress (S33) distribution on
the heat exchanger model for h = 0.1524 mm and a = 1.70 mm. Figure 16 (b) is a zoomed
in view of one of the channels of the same model. Figures 14 and 15 are S33 stress
contour plots on the test coupon for the same h = 0.1524 mm but different channel spans.
We see that the stresses for the same model (h = 0.1524 mm and a = 1.70 mm) on the

heat exchanger model and the test coupon are comparable.

2.5 DISCUSSION

From the FEA model we can predict that when the value of the S33 critical stress reaches
+ 0.2 MPa, then leakage takes place. Also when there are a few nodes in the bridge
region in tension and a few in compression, and when the value of the tensile stress on
these nodes does not exceed + 0.2 MPa then the layers are just bonded. However, we can
safely say that the layers will be bonded together when the stresses present in the bridge
region are purely compressive. The results of the FE simulation show that simple elastic
FE models with a set of lamina thickness and channel span without modeling the
complex creep behavior of the lamina allow you to determine the interfacial stress
between the lamina. Based on the sign of the stress, we can evaluate if diffusion bonding

can take place or not.

30

Aspect Ratio (a/h)

  

LEAKAGE

   
 

Grannel Span a (m)

BONDING

 

 

ﬂ ﬁr

04H—~ r . f w 1 v v
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Lamina Thickness h (m)

Figure 18: A plot of the experimental results showing the conditions under which
bonding and leakage occurred. Bonding conditions for all of these samples were
5.86 MPa at 900°c for two hours

20.
mi
16
w

12

  
 

 

 

 

 

 

10 LEAKAGE

8.

6‘ _

4* A BONDING

2i

0*“— ~- ——— . . —_7—.__--w.- . . . . . — . .
o 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Lamina Thickness h (m)

Figure 19: A plot of the experimental results showing lamina thickness Vs aspect
ratio

31

2.6 CONCLUSIONS

An investigation on limits of aspect ratio for counter-ﬂow heat exchangers was carried
out in this study. The simulation results were verified by carrying out the leak tests on
test coupons. Following are the conclusions of this study:

1. The fins in two-ﬂuid micro channel heat exchangers have a certain aspect ratio limit
for a certain given thickness of the lamina, beyond which poor bonding will occur within
the device resulting in leakage and mixing of the two ﬂuids in the device.

2. For a given set of bonding conditions, the width of the channel can be interpreted as a
function of the thickness of the lamina with the key bonding characteristic being that only
compressive out of plane stresses exist in the critical neck region. It is expected that the
laminae may just bond even if a tensile stress up to a maximum of 0.15 MPa exists in few
places in the critical neck region, provided that on all the other places of the neck region
the out-of plane stress is purely compressive.

3. It is expected that higher aspect ratios can be accommodated at smaller scales with
high aspect ratio (greater than 30) ﬁns possible at lamina thicknesses less than 1.36 um

for this set of bonding conditions.

32

3. F INITE ELEMENT ANALYSIS OF A OF F UNCTIONALLY
GRADIENT MESO-SCALE HEAT EXCHANGER

The miniaturization of electronic components and the rapid increase in power density of
advanced microprocessors and electronic components have created a need for improved
cooling technologies to achieve high heat dissipation rates. Some of future
microprocessors and power-electronic components have been projected to dissipate over
1000 W/cm2. The traditional methods of electronic cooling such as conduction and
natural/forced convection have proved to be insufficient for such high heat ﬂuxes. The
need for new technologies capable of dissipating high heat ﬂuxes and requiring low input
power is of critical importance to many industries. Coupled with these are requirements
of low volume/weight addition, high portability, high reliability and large-scale

fabricability of such future cooling devices.

The cooling system explored in this project was the ﬂuid ﬂow through an array of
parallel micro-channels. Among the advantages of this cooling method is the closeness of
the micro-channels to the heat source and therefore the thermal resistance is reduced.
Among the disadvantages is the non-uniform variation of heat ﬂux throughout the body
of the heat exchanger. This in turn causes uneven expansion which results in high
induced thermal strains, warping and even cracks. We propose to use Functionally

Gradient Materials to offset this effect of thermal gradient loading.

33

3.1 BACKGROUND

A key issue we face in our design is to distribute these FGM layers such that the thermal
strains are minimized. Considering the case of a long rod with a linear CTE, the change

in length due to temperature change is given by:
AI = laAt

where,

l = length of the rod

AI = change in length due to change in temperature (heating or cooling)
or = coefﬁcient of thermal expansion of the material of the rod

At 2 change in temperature of the rod

Hence the induced thermal strain along the length of the rod is given by:

82—Al—lzaAt

In our case, the bottom surface of the device is subjected to a high heat ﬂux [Fig 5]. Now
if the whole device is made up of the same material, then the part very near the high heat
ﬂux part tends to expand more than the rest of the body. Hence the objective is to build
the device using FGMs such that the layer closest to the heat ﬂux has the least coefficient
of expansion. Alumina and Zirconia powders were used to make the FGM. Because of
the almost similar sizes of the alumina and zirconia particles, the Mori-Tanaka
mathematical scheme was solved once for each of the two as the matrix material. The

average values were taken as the final material properties of the FGM.

34

3.2 F UNCTIONALLY GRADIENT MATERIALS

The unique feature of the heat exchanger is network of channels and manifolds in the
micro-textured medium achieved by spatially varying the compositional ratio between
two powders during conventional and microwave powder processing. The micro-textured
medium is necessary to compensate thermal gradient loading as a result of heat transfer
into a cooling ﬂuid circulating in the network of channels and manifolds. The residual
stress, induced by spatial variations in properties and densification, is minimized by using
different powder characteristics. The fugitive phase in the shape of desired channels is

used during compaction, which decomposes during sintering leaving cooling channels.

A key issue in producing functionally gradient materials is to avoid damage due to
process-induced residual stress. Residual stresses in FGMs [27] occur due to:
(i) Differential shrinkage in between the co-sintering layers due to their different
sintering kinetics during processing (densification behavior)

(ii) Thermal expansion mismatch

In order to resolve this problem the difference in the two shrinkage mechanisms existing
in powder processing must be utilized, namely the effect of CTE and densiﬁcation
behavior. During densification, the interstitial spaces among the powder are ﬁlled due to
mass transport mechanisms. The shrinkage differences due to CTE mismatch are intrinsic
for a given combination of materials after consolidation. However, the shrinkage from the
densiﬁcation of powder can be altered by designing the powder characteristics and
optimizing the sintering conditions such as time, temperature, heating/cooling cycles, etc

[31].

35

The starting powder characteristics and sintering conditions can compensate for the CTE
difference by exploiting the differences in densification behavior among the layers.
Powder shrinkage during densification must be understood as a function of various
processing variables such as the powder characteristics as well as sintering conditions.
The powder characteristics inﬂuence the initial powder packing density which in turn
affects the densiﬁcation behavior during sintering. Methods to control the initial powder
packing include mixing a few batches of powders with different average particle sizes

[32].

3.3 MORI-TANAKA SCHEME

Rather than microscopically modeling each inclusion in a composite medium, a given
inhomogeneous material can be replaced by an equivalent “effective medium”. Its
effective properties can be determined by applying the appropriate homogeneous
boundary condition with physical properties of the constituent materials. In particular, the
Mori-Tanaka theory [33] may be used for this purpose. Under the assumption of Mori-

Tanaka, the effective stiffness tensor is,

C“ =CMT =Cm +v<(c, —cm)AMT >

where,
v = represents the volume fraction of particles or fibers

C m = stiffness tensor of the matrix material
Cf = stiffness tensor of the particle (fiber) material

AMT = strain concentrator

36

The bracket < > represents orientation distribution function (ODF) — weighted
orientational averaging. ODE — weighted averaging for anisotropic and ellipsoidal
inclusions was extensively discussed by Ferrari and Johnson (1989). Limitations on the
applicability of the Mori-Tanaka theory have been discussed by Ferrari (1991). For the
FGM body made of isotropic matrix and isotropic spherical inclusions, both the Mori-

Tanaka and the poly-inclusion method predict the effective stiffness to be

Ct :CMT =Cm + V(Cf _Cm )AMT

VI T

Here, the strain concentrator, A‘ , is defined as

AMT =T[(l — v)I + m“

where,
I = fourth-rank identity tensor

T = Wu’s tensor

Wu’s tensor (1966), is defined as

T=[I+ESm(C, —Cm)]“

where,
E = Eshelby tensor

Sm = Compliance tensor of the matrix material

The Eshelby’s tensor for spherical inclusions written in 6 X 6 matrix form is

37

 

 

 

 

 

_ _ 5 _ _ 7
7 -v 51/ 1 5v 1 O 0 0
15(1—v) 15(1—v) 15(1—v)
5v —1 7 — SV SV —1 0 0
15(1-v) 15(1- v) 15(1 —v)
5v — 1 5v —1 7 — 5v 0 0
15(l—v) 15(1—v) 15(1—v)
[E]: O 0 0 ﬂ 0 0
15(1 —v)
0 0 0 0 ——4 7 5"
15(1 —v)
0 0 0 0 :11;
15(1—v)
where,
v = Poisson’s ratio of the matrix material
The stiffness tensors (Cm , Cf ) are defined as,
—Cll C12 C12 0 O O T
C.. C.. C. 0 0 0
C12 C12 Cl! 0 0 0
1C, 1 =
0 0 0 C11 — Clz/2 0 0
0 0 0 0 CH — C,,/2 0
L 0 0 0 0 0 C11 — Cu/2_

 

 

where CH and C12 are related to Lame’s constants defined as,

Cll = X + 2n

C ,2 = it

and where xi and u are related to Young’s modulus (Y) and Poisson’s ratio (v ),

xi = Ev/{(1+ v)(1—- 2v)} , and,

38

221 = E/(l + v)

In the Mori-Tanaka scheme, the effective compliance tensor, SMT , may be obtained as the

reciprocal of CMT [34].

For thermal gradient applications, the effective thermal expansion coefficient (CTE)

tensor, a/+ , is required. The CTE corresponding to the Mori-Tanaka approach is

a!“ zam + (a, —am)(S, -S )“(SMT —S )

II] III

In the above equation,

am = CTE of the matrix material

a, = CTE of the particle material

Sm = Compliance tensor of the matrix material
S, 2 Compliance tensor of the particle material

SMT = The effective compliance tensor

For steady—state thermal conduction, the expression for the effective thermal conductivity

for the composite spheres assemblage (CSA) model was obtained by

K" =Km + v(K, — Km )AMT

where,

K m = Thermal conductivity of the matrix material

K, = Thermal conductivity of the particle material

39

v = Volume fraction of particles or fibers

AMT = Strain concentrator
The specific heat of a specimen at temperatures below ambient temperature can be
measured by a dynamic heating method. It has been found that except for the specific

heat at 88K, those measured are close to the value predicted by the law of mixtures, i.e.
C = VmCm + V_,C,

where,

Vm 2 Volume fraction of the matrix phase

V, = Volume fraction of the particle phase

Cm = Specific heat of the matrix material

C, = Specific heat of the particle material

Density of a material at a constant temperature can again be predicted by the law of

mixtures, i.e.
p: Vmpm +prf

where,

pm = Specific heat of the matrix material
,0, = Specific heat of the particle material

The powders used to model our FGM, were alumina and Partially Stabilized Zirconia
(PSZ). Two-layer (alumina and P82) have been successfully made [35] in the past. The

ﬁrst layer was made with the alumina powder designated TM-DAR, manufactured by

40

Taimai Co. Japan. The second layer was formed with a mixture of 30% CERAC, a kind

of Fully Stabilized Zirconia (FSZ), and 70% P82 (TZ-3YS).

Using the above relations, the material properties of the composites with varying volume
fractions of TM-DAR alumina and TZ-3YS partially stabilized zirconia were determined.
Matlab was used to solve the matrices in the above equations [Appendix A]. Because of
the similar sizes of the two particles, the equations were solved once taking almunia as

the matrix phase and then taking PSZ as the matrix. Table 2 gives the material properties

 

of the two samples.

Table 2. Material property values of alumina and yttria-partially stabilized zirconia

 

 

Property Alumina Y-PSZ

Density (g/cc) 3.97 6.02

CTE (/°C) 8.1 x 1076 10.3 x 10‘6
Young’s modulus (GPa) 350 200
Poisson’s ratio 0.23 0.3

Thermal conductivity (W/m/K) 20 3

Speciﬁc Heat (J/Kg/K) 750 450
Flexural Strength (MPa) 355 660

 

 

3.4 DESCRIPTION OF DEVICE

The geometry of the proposed device is as shown in Figures 5 and 6. This device is
primarily meant to cool electronic chips or cutting tools. It is subjected to a high value of
heat ﬂux on its bottom surface. The primary means of cooling this heat exchanger is

through liquid coolants ﬂowing through internal channels in the device. The device also

41

looses heat to the ambient air through convection from its side walls and its top surface.
A ﬂux of 250,000 W/m2 enters through its bottom surface. The coolant used was water

and step-wise varying linear properties with temperature were assumed.

Water ﬂows in all four channels with an inlet velocity of 0.5 m/s along the length of the
channels. The water is discharged at the outlet to the ambient at atmospheric pressure.
The ambient air is assumed to be at 25°C and its convective heat transfer coefﬁcient is
approximated to be 60 W/ m72 /K. The objective is to prevent solid temperatures from
raising over 50°C. At the same time we also want to minimize the induced thermal strain
in the body. As discussed earlier, the material of the lowest layer and the one nearest to
the heat source should have the least CTE value. An attempt was made to study the
values of thermal strains for different ratios of thicknesses of the FGM layers for a three

layered model [Figure 20].

 

   

100°/
90% .
80°/- '
70% ,
60% ,
50% 7

40% ‘_
30% :
20%

10% I
0% -

1 2 3

Figure 20: Different layer thickness ratios for a three-layered FGM model
(Equal thickness layer model, Model 1 and Model 2)

42

For the three layered model, the bottom layer (layer 1) was alumina, followed by a
mixture of 50% alumina and 50% zirconia (layer 2), and finally a layer of pure zirconia
(layer 3). The combination with the least induced thermal stress was chosen. Then, the
FEA simulation was repeated for a six-layered FGM model with again the design
objective being to choose the one with the least induced thermal stress. The bottom layer
is named layer one and the top layer is layer six. Again the six-layered model was built
proportionately as follows:

Layer 1: 100% alumina

Layer 2: 80% alumina + 20% zirconia

Layer 3: 60% alumina + 40% zirconia

Layer 4: 40% alumina + 60% zirconia

Layer 5: 20% alumina + 80% zirconia

Layer 6: 100% zirconia

 

   

1 00%
90

 

 

 

 

 

0%

Figure 21: Different layer thickness ratios for a six-layered F GM model
(Equal thickness layer model, Model 1 and Model 2)

43

 

The properties of the three and six layered FGM, solved from the Mori-Tanaka model,

are tabulated in Appendix C.

3.5 FINITE ELEMENT MODEL

Our objective here is two-fold:
i. Solve the ﬂuid-thermal interactions for temperature distributions in the FGM
model as well as in the cooling ﬂuid.
ii. Solve the thermal-structure interactions in the solid part for induced thermal

strains.

We made use of Ansys Multiphysics to solve these coupled ﬂuid-thermal-structure
interactions. Care was taken to use only three-dimensional brick elements to mesh our
models especially when computing ﬂow through ducts. This is because it is important to
align the grids along the direction of ﬂow as they exert considerable inﬂuence on the

ﬂow and heat transfer.

For the ﬂuid-thermal analysis, 3-D CFD ﬂotran elements were used. The segregated
solver was used for ﬂow and TDMA solver was used for the velocity calculations. Here
the thermal properties of the non-ﬂuid material are several orders of magnitude different
from those of the ﬂuid, and this is called an ill-conditioned conjugate heat transfer
problem. In this situation, the TDMA method probably will not yield useful results, no
matter how many sweeps we specify. The most robust but most memory-intensive

method for solving conjugate heat transfer problems, the Preconditioned Generalized

 

Minimum Residual method, was used as the temperature solver. Preconditioned
Conjugate Residual method was chosen was the pressure solver because of its robustness
for solving ill-conditioned heat transfer problems. Appropriate boundary conditions were
applied (Figure 5) and the model was solved till desired convergence was achieved. Now
the 3-D CFD ﬂotran elements were replaced with 3-D structure elements in the solid
regions, while null element was used in the ﬂuid part. Then temperature boundary
conditions were applied by importing the results file from the previous ﬂuid-thermal
analysis. Now the system was solved for thermal strains only in the solid regions. The
analyses were performed on all the models of the three layered as well as six layered

specimen.

3.6 RESULTS AND DISCUSSION

The results that were obtained from the finite element analysis were in agreement with
our expected design objectives. First, the analysis was carried out on an alumina model

and then the results were compared with the three layered and the six layered model.

3.6.1 Alumina model

Without the cooling channels, it was observed that the steady state temperatures were in
the range of 280 °C to 250 °C. With the cooling channels we could bring down the steady
state solid temperature range to 118 °C to 59 °C [Figure 22]. After the thermal-ﬂuid
analysis, a structural analysis was performed to determine the maximum principal
stresses in the model [Figure 25 and 26]. It was observed that in the regions of high
temperature gradients, the stress was the highest. These images in this thesis are

presented in color.

45

 

 

 

Figure 22: Temperature contours of the 8 -channeled alumina model

Figure 23: Maximum principal stress contours on the 8 -channeled alumina model

46

3.6.2 Comparing Results between - pure alumina, 3 layered FGM and 6 layered
FGM

On comparing the results from the three models, it was found that the stress was least in
the six layered FGM model, followed by the three layered model and it was maximum in
the pure alumina model. Also for the models with the same number of layers, the ones

with the smallest ﬁrst layers produced the minimum stresses.

3.6.3 Maximum principal stress contours

1.00E+O7 —
8.00E+06 ~
6.00E+06 .
4.00E+06 .

2.00E+06 a

 

..--.----'- Alumina
__ 3 layered FGM
_ . . - 5 layered FGM

 

0.00E+00

 

 

 

-2.00E+06 ~

-4.00E+06 .

Maximum Principal Stress (N/rn2)

'6.00E+06 ‘ :

 

-8.00E+06 ~

 

-1.00E+07 -

Height along Zaxls (mm)

Figure 24: Maximum principal stresses at x = 0, y= 0
(F GMs have equal layer thickness)

47

 

 

 

Figure 25: Alumina model - maximum principal stress contours at x=0

Figure 26: Three equi-layered FGM model - maximum principal stress contours

at x=0

48

5.00E+07 T'WTT'""‘“ ,, a» - . . , -.__--. m.-.
4.00E+07 .

3.00E+07 4 r

 

2.00E+07 ~ : Alumina

—— 3 layered FGM
- ~ ~ - 6 layered FGM

 

 

1.00E+07 -

0.00E+OO -——/—,4— -— -- “ ' a
0 ,3 1 2 3 4 \0

-1.00E+07 *

 

 

Maximum principal stress (N/m2)

 

 

 

-2.00E+07
Height along 2 axis (mm)

Figure 27: Maximum principal stresses at x = 0, y 2 3mm
(FGMs have equal layer thickness)

From the above results it was seen that by building models using functionally gradient
materials, we succeeded in reducing the tensile stresses by almost half. Now, each of the
two FGM models — the three layered and the six layered, were analyzed to see the effect
in stresses by varying the thickness of each layer. The effect was not as pronounced as
seen between an alumina model and a FGM model. However the effect was seen more in
the case of the three layered FGM than the six layered FGM. For both the cases, the

models with the smallest ﬁrst layers (model 1) produced the minimum principal stresses.

49

 

6.00E+06

4.00E+06 .

2.00E+06 ~

 

 

0.00E+OO
0

-2.00E+06 *

-4.00E+06 a

Maximum Principal Stress (N/m2)

-6.00E+06 i

-8.00E+06 ~

 

 

 

 

 

-1.00E+07

Height along Z axis (mm)

 

 

__ Equally thick layers
_ Thin botom layer
_ . . - Thick bottom layer

 

 

Figure 28: Maximum principal stresses for the 3 layered F GM model

6.00E+06 -

4.00E+06 .

2.00E+06 i

atx=0,y=0mm

 

 
 

 

0.00E+00
(O

-2.00E+06 ,

-4.00E+06 4

Maximum Principal Stress (Nlm2)

 

-8.00E+06

 
 
  

 

 

-1 .OOE+O7 ~

Height along Z axis (mm)

 

 

Equally thick layers
— Thin bottom layer
- - - - Thick bottom layer

 

 

 

Figure 29: Maximum principal stresses for the 6 layered FGM model

atx=0,y=0mm

50

 

Maximum principal stress (N/m2)

Maximum principal stress (N/m2)

1.00E+07 -. .

8.00E+06 ~

6.00E+06 *

4.00E+06 <

2.00E+06 .

0.00E+00

 

-2.00E+06 .

-4.00E+06 *

-6.00E+06 *

-8.00E+06 —

(I)

 

   

 

-1.00E+07

Height along Z axis (mm)

  

r.— .._

I l
f I—Thin bottom layer
61m.

1

 

 

Equally thick layers

- Thick bottom layer

Figure 30: Maximum principal stresses for the 3 layered F GM model

8.00E+06

6.00E+06

4.00E+06 .

2.00E+06 i

0.00E+00

-2.00E+06 -

-4.00E+06 ~

-6.00E+06 ~

-8.00E+06

-1.00E+07 .

 

  

atx=0,y=3mm

 

 

(I)

 

    
 
  
  

 

5

 

 

Height along Z axis (mm)

 

 

 

Equally thick layers
-— Thin bottom layer
— - - - Thick bottom layer

 

Figure 31: Maximum principal stresses for the 6 layered FGM model

atx=0,y=3mm

51

 

 

 

 

Many studies in the past have looked at the effects of confining the brittle material [36].
In our case, the functionally gradient materials were modeled to fail at a tensile stress of
355 MPa which is also the tensile strength of alumina (lower of the two between zirconia
and alumina). At this value of stress the material underwent brittle failure. Elastic perfect
behavior was assumed for compressive behavior of the functionally gradient material and
the Rankine failure criterion was used in the models. In all the cases, it was seen that the

maximum principal stress in the models was well below the tensile limit.

For obtaining mesh convergence, the mesh was refined in all the areas until the required
convergence limit was obtained. It was verified that the ratio of stress values between two

successive reﬁnements was less than 0.1%.

3.7 VALIDATION

For validating our results, the finite difference approach was used. A simplified case of a
two dimensional two layered FGM model of alumina and P82 zirconia with only one
cooling channel was considered. Starting from the simpliﬁed momentum and energy
differential equations for fully developed (velocity and temperature proﬁles) ﬂow in the
two-dimensional parallel wall channel, the convective heat transfer coefﬁcient was
calculated when there is a constant ﬂux through the lower wall but the upper wall is
insulated. The general steps that were followed are as follows:

i. The laminar-ﬂow solution to fully-developed channel ﬂow was found and it was

expressed in terms of the mean velocity V", . For simpliﬁcation, y=0 was chosen to

be the mid-point of the channel.

52

 

From the continuity equation, we have:

(it By
The X- momentum equation can be simplified to:

1810 all!

___ _=0

,0 81‘ By 3

And the Y-momentum equation to:

“I..- 41.1.3:ch “I!”

__Li’izo

p 3y

 

On applying the boundary condition, u = 0 at y = -_i- h, and g}: = 0 at y = 0, we get:
)1

u =——L-a—I:(h2 — yz)
2,11 8x

Now, the area weighted mean velocity is defined as,

If}, udy

"I

 

h
dy
«h
Solving the above equations we get,
3 2

y
“__v 1__
2 m( hz)

The mixed-mean ﬂuid temperature Tm for channel ﬂow was defined and an

overall energy balance (integrating over y the advection-diffusion equation) was

calculated to relate q.” to dTm /dx

53

7

j for low speed

0

a:

for the long pipe and also the term ”[8
)1

 

Neglecting the term ,

H

ﬂows, the thermal energy equation gives us:

#0_T_ad2T
8x 8313

 

Now let’s define the mixed-mean ﬂuid temperature Tm as,

F quy
Tm(.r)= ',',
I udy

—h

 

On solving the above equations and applying the boundary condition q”, = constant on
the lower wall and ‘10... = 0 on the upper wall, we get:

37 . ”
V 212 m = A in
pCp .. ax q

 

iii. The dimensionless temperature profile is now defined as,

 

. T .
6: “ T ,where Tw rs the wall temperature

For a fully developed temperature profile, by putting %¢9_ = 0 and solving, we get:
x

2-8T

M'

__ 8T“, _ 8T
ax dx

"1

33c— 8x

and

 

iv. By using all of the above results, and on integrating the advection diffusion

equation twice we can ﬁnd T(y). Then using the deﬁnition of Tm and substituting

the above result, we get:

54

 

 

V dT ,
H1 ”(2)
a dx 35

. l/
v. Combining all the results and using q = H (T, — Tm) , where H is the convective

heat transfer coefﬁcient, we get:

_ﬂ
26h

In our case, we consider a channel with cross-section of 4mm X 2mm. Further, for the
sake of simplicity in calculations we consider our block to be 16mm X 16mm X 6mm in
dimension, and simplify it further to a 2D model with 16mm X 6mm dimensions.
Adopting a grid space of Ax = Ay =1 and identifying the line of symmetry, the foregoing
nodal network is constructed. The corresponding ﬁnite-difference equations may be

obtained by applying the energy balance method to nodes 1 to 61.

Th

 

H

Figure 32: The two dimensional model for validation

Q in

55

 

 

Figure 33: Meshing the symmetric model for finite-difference calculations

For all interior nodes, not heat ﬂux ﬂowing in is equal to zero. Referring ﬁgure 33,

4

; q(i)———)tm.n) — 0
I:

So by using simpliﬁed form of the Fourier’s law, the rate at which energy is transferred

by conduction from node m-l, n to m, n may be expressed as,

T — T
= k(Ay. 1) —”'“"—’“‘ , where k = conductivit of the material
Ax y

q(m—l,n)—r(m,n)
Applying similar approximations on all nodes in a region with the same material
properties, we obtain the following equation for all those interior nodes:

+ Tuna—1

+T

m+l.n

+T

m—l.n

T

m.n+l

— 4va,, = 0

56

For nodes from 1 to 7, the equations cannot be simplified as above because of different k
values. For these nodes we use two different values of k, k1 and k2, and use the basic

finite difference equation to solve for them.

Heat transfer to node 13 occurs by conduction from nodes 20 and 22, as well as by
convection from the outer ﬂuid. Application of an energy balance to this section yields a

finite-difference equation of the form,

 

l
(T20 +733) + 2%1, — 2P?”- +1)T,, = 0, where,

l' ‘lm I Jul . . .

h0 = convective heat transfer coefﬁcient of air

To, = bulk air temperature

k = conductivity of the material — either kl or k2

Nodes 14 to 20, l, 21, 22, 48 and 49 have the same physics (3 conduction, l convection)
and application of an energy balance to node 14, as an example, yields a finite-difference
equation of the form:

2110M

 

Ax
T11+ T15 + 2T24 — 2['hOT+ ZJTM = — Two

For nodes 2, 10 and 47 h0 is replaced with H of water, and T,” with Tm".
For nodes 12 and 45 (3 conduction, l adiabatic), taking node 12 as an example, it follows

that,

T9+2T24+Tn —4Tl2 =0

57

Nodes 8 and 46 are internal comers for which, taking node 8 as example,

T

”Ax jr. = 451:3:

2T32 + 2T34 + T2 +T,0—2(3 + T

On the bottom line, a heat ﬂux Q is applied. So, each element will have a heat ﬂux of
Q/8. So, nodes 37 to 43 have a heat ﬂux of Q/8 and nodes 35 and 36 have a heat ﬂux of

Q/l6. Again using the principle of heat balance on these nodes we get,

4

= 0
Zq(ii—)(m.n)

i=0

Applying the above relations we can solve for nodes 35 to 43.

 

a"

In case of nodes 9 and 44, heat transfer is characterized by conduction from two
adjoining nodes and convection from the internal ﬂow, with no heat transfer occurring
from an adjoining adiabat. Performing an energy balance for nodal region 18, it follows

that:

 

H Ax HAx
(T12 + T17 ) _[ 2)7718 2 _ TTwnter

For node 1 1, H of water is replaced with ho, and T with T0,.

water

For validation purpose a low value of heat ﬂux, 10,000 W/mz, was used to minimize the

error due to a coarse mesh. Also, from the relation,

H = :33]:- , by plugging in the value of h = 1mm and k = 0.6 W/m/K,

26h water

we get H = 807.7 W/mZ/K.

58

Using hm: 15 W/mZ/K, and T0,: 25°C, we solve the above finite-difference equations.
The nodal temperatures obtained from the above calculations are tabulated in Table 3.
These results are then compared with those obtained from a ﬁnite element model solved
in Ansys 7.0. The FE model used is a 3-D cylinder with radius of 8 mm and a thickness of
6mm. An inlet velocity of 0.3 m/s (Re = 800) and inlet temperature of 25°C for water was

used with outlet pressure set to atmospheric. Convection was modeled all external

surfaces except the bottom surface where a heat ﬂux of 10,000 W/m2 was applied.

 

Figure 34: Finite element model with one cooling channel, used to compare results
obtained from the finite difference calculations

It can be seen that considering the very rough estimate of the model, the nodal
temperatures are within reasonable comparable limits. The Ansys ﬁnite element model
was also solved for thermal strains by using a 2D coupled element. The results were then

compared with the hand calculated results.

59

 

For hand calculation of stress along the X axis, first the temperature was averaged in the
area with material 1, from elements I to 22, and the steady state averaged temperature
denoted by T,. Then similarly the temperature was averaged in the area with material 2,
from elements 23 to 44, and denoted by T2 . The materials would tend to elongate along
the X axis, and their final elongated lengths if allowed to expand can be defined as:

L.- =1(1+ a, (T, — 25)), where: T.
I
l

a, =Coefficient of thermal expansion of the materials

 

l = Initial length = 16 mm, and,
T, = Final steady state averaged temperature of the materials i
Depending on the value of 0:, one material would tend to expand more than the other.
Now, this is assumed to be similar to a bimetallic strip, and the final elongated length of
both areas is assumed to be L", . To achieve equilibrium, the compressive force exerted
by one area should equal the tensile force of the other. So, in other words considering the
2D layers to be also [mm wide (normal to the paper), we have:
0',t,l = Uztzl , where II = t2 = respective layer thicknesses = 3 mm
or 0'l = 0'2 and E18, = E282, where
El and E2 are the respective Young’s modulus, and,

8, and 8, are the respective strains

L . . . .
Now, 5,. = ”' ‘ , where subscript 1 denotes either layer 1 or layer 2.

I

 

Hence, on solving the above equations, we get:

L = 1.1.(E.+E.)
'" (E.L.+E.z.)

60

Table 3: Comparing nodal temperatures from the hand calculated FD model with
the Ansys FE model

TF0 (°C) = Nodal temperature hand calculated using the finite difference scheme

TFE (°C) = Nodal temperature obtained from the Ansys finite element model

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

TF D TFE TF1) TFE TF D Tl-‘E TF D TFE
Node (°C) (°C) Node (°C) (°C) Node (°C) (°C) Node (°C) (°C)
1 45.3 34.7 17 41.5 33.7 33 42.9 33.1 49 45.2 34.6
2 41.7 33.9 18 42.6 32.5 34 39.1 29.6 50 45.1 34.4
3 45.2 34.7 19 43.1 32.8 35 45.7 34.7 51 44.9 34.1
4 44.9 34.5 20 43.4 33.3 36 44.1 34.3 52 44.8 33.9
5 44.5 33.9 21 44.2 33.7 37 45.4 34.7 53 44.4 33.7
6 44.1 33.6 22 43.8 33.7 38 45.2 34.6 54 44.3 33.5
7 43.6 32.6 23 43.6 33.6 39 45.1 34.4 55 44.1 33.4
8 37.4 27.7 24 36.5 27.6 40 44.8 34.3 56 45.1 34.4
9 34.7 26.0 25 40.2 32.4 41 44.6 34.2 57 45.0 34.4
10 35.2 26.8 26 41.5 33.7 42 44.3 34.1 58 44.7 34.4
11 36.1 27.4 27 41.9 34.1 43 44.2 34.1 59 44.2 34.0
12 36.0 27.3 28 42.7 33.0 44 42.9 32.8 60 45.1 34.4
13 43.8 32.7 29 44.2 33.7 45 43.8 33.1 61 44.3 33.5
14 37.4 27.7 30 43.8 33.7 46 42.6 32.6
15 38.2 29.4 31 43.5 33.5 47 42.1 32.4
16 40.2 32.4 32 41.3 32.8 48 45.4 34.7

 

 

 

 

 

 

 

 

 

 

 

 

61

 

Now, say layer 2 tends to expand more than layer 1 along the X axis. So compressive

stress in layer 2 would be:

 

0'2 = —E,[L3;L"'J, where L, > L", > Ll

And tensile stress in area 2 would be:

a, = +E,[L"’ _12]
L,

 

Plugging in the value, the tensile stress in layer 1 (bottom layer) was found to be 1.87 E6,
and the compressive stress in layer 2 (top layer) was found to be — 1.87 E6. The average
tensile stress in layer 1 from the FE Ansys model was found to be 2.57 E6, and the
average compressive stress in layer 2 was - 2.87 E6. It was thus seen, that considering the
very rough estimate of the model, the elemental stresses were within reasonable

comparable limits when compared with the finite element model.

3.8 CONCLUSIONS

The idea of using 8 cooling channels was useful as it could bring down the temperature in
the operating range of the device. It was observed that the high thermal gradient created
by the cooling channels results in high induced thermal stresses in the specimen. Hence
to offset this effect the concept of using a material with functionally gradient properties

was used.

From the results of the analysis, it can be inferred that the stresses in the FGM models

were considerably less than the stresses in the alumina model. Hence it can be said that

62

 

the idea of using the FGM can be justified. The stresses in the six layered FGM model
were further found to be lower than those in the three layered FGM model. Varying the
thickness of the layers of the FGM also revealed interesting results. For both cases of a
three layered and a six layered model, it was found that the model with a thinner bottom

layer had a lower value of stress than the other combinations.

3.9 RECOMNIENDATIONS FOR FUTURE WORK

As far as future work on this topic is concerned, currently attempt is being made to
powder process a three layered FGM model. Work can be done on ﬁnding what should
be the optimum ratios of layer thicknesses. Experimental validation of the heat exchanger
should of course be the final goal. Graphite can be used as the fugitive phase for making

the internal channels and the heat exchanger can then be tested under actual conditions.

63

 

4. COMPUTATIONAL STUDY OF MIXING MICROCHANNEL
FLOWS

Ever-improving technological advances in micro and nano fabrication have encouraged
renewed interest in mixing through topological or chaotic ﬂows. Particular attention has
been paid to geometries that are sufﬁciently small because neither turbulent (owing to the
low Reynolds number), nor active (owing to the difﬁculty in micromachining moveable
or easily controlled parts) mechanisms are feasible. Examples of recent successful
strategies include a “staggered herringbone” nricromixer [37] in which two species were
mixed well beyond the diffusion limit simply through patterning the bottom of a channel

200 microns in width.

Motivated by recent experimental advances [37, 38] in microﬂuidic mixers, we study the
passive mixing and ﬂow properties of a ribbed micro-channel by means of computational
ﬂuid dynamics (CFD). Such geometries overcome the low Reynolds number boundaries
to efﬁcient mixing by creating a three dimensional ﬂow that yields chaotic trajectories. It
is hoped that such CFD studies advance both the capabilities and the understanding of
such micro-pattemed mixing ﬂows. After studying several properties of the ﬂow which
are not experimentally accessible, we studied mixing in: (i) channel with two pairs of
inclined ribs (ii) identical channel with no ribs. We also study the pressure drop along the
length of the channel, which limits the velocity and therefore the time necessary for
effective mixing. Advantages of full CFD studies of micrornixing include the ability to
visualize species densities at any point in the volume without optical challenges of

confocal-microscopy.

 

“j. 7'.-- .

4.1 CHANNEL DESIGN AND SPECIFICATIONS

To study the differences in mixing rates two geometries of microchannels, one of them as
shown in ﬁgure 10, were built. All the channels are square in cross-section with sides of
1.6 mm and are all 20 mm long. The straight ribless channel only allows one-
dimensional ﬂow of the ﬂuid. On the other hand, the ribbed channel has a 3-dimensional
geometry with two pairs of inclined ribs. Near the inlet the ribs were on the top and
bottom sides of the duct, while near the outlet the inclined ribs were placed on the two
opposite side walls. Each ribbed pair consists of ten ribs in each straight duct, ﬁve on the
leading wall and ﬁve on the trailing wall, all with the same pitch. The ribs on those two
walls are staggered relative to each other with he ribs on the leading wall offset from
those on the trailing wall by a half pitch. All ribs are inclined with respect to the ﬂow at
an angle of 45 degrees. The cross-section of the rounded ribs is made up of three circular
arcs of radius R, where R equals 0.08 mm so that the rib height is 0.16 mm and the rib-

height to hydraulic-diameter is 0.1. The pitch to height ratio is 5.

4.2 MESH GENERATION

When computing ducts with ribs, it is important for the geometry of the ribs to be
captured correctly. This is because they exert considerable inﬂuence on the ﬂow and heat
transfer. Shih, et al [30] studied the effect of a ribbed U-duct with or without rotation.
Their model consisted of a U-duct with inclined ribs on just two Opposite faced walls.
They developed a mesh capturing the effects of the side ribs to perform computational

ﬂuid dynamics simulations.

65

 

We used the code written by Shih et al [30], which was modiﬁed to represent our
geometry. The code by Shih et a1 [30], which was in NASA plot3D format, was ﬁrst
converted to ASCII dat format readable by Gambit a CFD meshing tool. Then a code in
C++ was written [Appendix A] to cut the grid and modify it according to the required
geometry. Later, the modiﬁed grid was imported into Gambit where the boundary

conditions were applied.

The domains of the smooth and ribbed ducts were replaced by H-H structured grids
(Figure 35). The number of grid points in the streamwise direction from inlet to outlet
was 255 for the smooth duct and 761 for the ribbed duct. Whether smooth or ribbed, the
number of grid points in the cross-stream plane was 65 X 65. The number of grid points
and their distribution were obtained by satisfying a set of rules such as aligning the grid
with the ﬂow direction as much as possible, keeping the grid aspect ratio near unity in

regions with recirculating ﬂow.

4.2 MODELLING DETAILS, INITIAL AND BOUNDARY CONITIONS

After modifying the mesh using our C++ code, the model was imported into Gambit V.
Using Fluent 6.0 speciﬁc boundary conditions, channel walls were deﬁned, the two inlets
were deﬁned as “velocity inlets” and the outlet was deﬁned as “outﬂow”. The volume of

the channel was deﬁned as “ﬂuid. The meshed models are shown in Figure 35.

The models were then checked and exported to Fluent 6.0 for performing ﬂow and
mixing computations. The Multiphase Solver in Fluent was used to solve for mixing of

two ﬂuids entering the microchannel. Initially two liquids with identical properties as

66

 

that of water were used to analyze mixing. Later, the mixing simulations were run for a
variety of combinations of different liquids. The properties of the liquids used are listed
in Table 4. A segregated steady-state solver with Mixture Multiphase model was used.

Standard operating conditions of pressure were speciﬁed.

 

 

 

 

 

 

 

Figure 35: Meshing of ribbed-duct

At inlet 1, a volume fraction of 1 was speciﬁed for liquid 1 and a volume fraction of 0

was speciﬁed for liquid 2. Similarly, at inlet 2, a volume fraction of 1 was speciﬁed for

67

'1'»

..l I III}.

liquid 2 and 0 for liquid 1. Equal velocity magnitudes ranging from 0.1 m/s to 1.5 m/s
were speciﬁed to investigate mixing. For the channel height and width speciﬁed, and the
hydraulic diameter calculated as 1.6 mm, the velocities correspond to Reynolds numbers
from 10 to 1500 depending on ﬂuids used. The outlet of the channel was given a zero-

boundary condition for all ﬂow variables.

Appropriate solution controls and residual monitors were set. After setting zero initial
conditions the model was iterated. Depending on the velocity magnitudes, about 600-
1500 iterations were required to achieve convergence. The end of the simulation was

veriﬁed by a quick check to see that the ﬂow ﬁeld for the entire channel was solved.

4.3 SIMULATION RESULTS

4.3.1 FLOW PATTERN

Velocity vectors and coordinate positions were plotted in Fluent to obtain vector plots.
To obtain the plot, a horizontal center plane was deﬁned for the ribbed and ribless ducts,
and the top view of the plane was plotted in Fluent. These plots were used to check for
the formation of secondary ﬂow and ﬂow separation for an indirect veriﬁcation of
chaotic advection. Figure 36, shows the velocity vectors plotted on the center plane along
the length of the ribbed duct for a water-water mixing at Re = 800. Figure 36, in this

thesis is presented in color.

68

 

veracny iagr tutfe irriisiur
FLL

er arms in." T SEEM

1 HT 13.1 130. tip. segregated rrritrure. larni

 

Figure 36: Velocity vectors of the ribbed duct for water-water mixing at Re = 800

A high speed core velocity is formed in the ribbed ducts with the highest velocities were
in the ribbed region. The chaotic high speed ﬂow in the ribbed region indicates mixing of

the two ﬂuids.

4.3.2 MIXING

For investigating mixing along the channels, contour plots of concentration or volume
fraction of liquid 2 were plotted in Fluent. The blue color at inletl indicates a volume
ﬁaction of 0 and red color at inlet 2 indicates a volume faction of 0 for liquid 2. Green
indicates a volume fraction of 0.5 which is the desired result for perfect mixing. The plots

show how the volume fraction approaches 0.5 as the two liquids mix as they ﬂow along

69

 

Lh m!“ Ihﬁ

 

the length of the channel. For the case of water-water mixing, it can be seen that the two
ﬂuids mix almost completely in the ﬁrst set of ribs (a) in the ribbed duct. In the ribless
duct (b), the ﬂuids hardly mix at all outside of their interfacial area along the center of the
channel. The water-glycerin mixture takes considerable longer time to mix and they are
nearly mixed only at the end of the second set of ribs. Figures 37 to 42 show the volume
fraction contours of one of the mixing ﬂuids. These contour plots in this thesis are
presented in color. A more detailed analysis of these simulation results and differences is

provided in the next section.

 

Figure 37: Contour plots of volume fraction of liquid 2 in water-water mixing

in the ribbed duct with inlet velocity = 0.5m/s (isometric view)

70

 

’1)‘t."ifi'?éill‘ll)l‘ rrlraspe ‘r

 

Figure 38: Contour plots of volume fraction of liquid 2 in water-water mixing
in the ribbed duct with inlet velocity = 0.5m/s (cross-sectional view)

 

_ all ILZ‘I'PI :ili. Z"‘.'>:'..l'~r,1 v-l 1..” ll

Figure 39: Contour plots of volume fraction of liquid 2 in water-water mixing
in a ribless duct with inlet velocity = 0.5m/s (cross-sectional view)

71

 

Guitars 0' "GIMME 'raclizxi rpl‘saseel

 

FLUE? IT E. ‘: rid, (ll). <etjregaretl.[ milk: 1
Figure 40: Contour plots of volume fraction of liquid 100 times more viscous than
water when mixed with water in a ribbed duct with inlet velocity = 0.5m/s
(cross-sectional view)

IE“ l«I [It];

.
I t .
PLUE" 1T {:1 13:1. (ll). segregated mixture. lair-.1

"clitoris :l “chime tr;

 

Figure 41: Contour plots of \wlume fraction of liquid 300 times more viscous than
water when mixed with water in a ribbed duct with inlet velocity = 0.5m/s
(cross-sectional view)

 

‘0'0U' re 'izec! Jl' ti; raise:

 

Figure 42: Contour plots of volume fraction of water-glycerin mixing in a ribbed
duct with inlet velocity = 0.5m/s (top view)

4.4 ANALYSIS OF RESULTS

To better characterize mixing in these channels and to provide a quantitative treatment of
analysis, “cross-sections” were deﬁned for all the models. Four cross-sections were
deﬁned — one each before and after each set of ribs, to observe how mixing occurs along
the entire depth of the duct and to see if recirculation in the plane perpendicular to the

direction of ﬂow occurs.

4.4.1 CHARACTERIZING MIXING

Mixing characteristics were analyzed in the cross-sections and simulations were run in

each of the two types of ducts for different inlet velocities and different combinations of

73

 

ﬂuid mixtures. For each run, volume fraction plots were obtained and the data values of
volume fraction of liquid2 calculated in the zone of nodes included in the cross-sections
were exported to MATLAB to obtain statistics. The summary of all runs is tabulated in
Tables 4 and 5. In Tables 4 and 5, the ﬁrst three columns tabulate the inlet velocities and
the corresponding Re values. The average values of volume fraction of liquid2 calculated
for the four cross-sections are tabulated. Clearly average values were found to be poor
indicators of mixing as can be seen for the case of the ribless duct where the values are
0.5. ‘Mixing’ here takes place only at the interface of the two liquids and the volume

fractions remain fairly constant along the duct at regions further apart from the interface.

4.4.2 MIXING RATES AND EFFICIENCIES

Since we are trying to characterize the mixing rate by the values of the volume fraction
which is a conserved quantity, as the volume fraction of a fully mixed ﬂuid is 0.5, the
standard deviation in composition was considered to be an approximate estimate for
mixing. This was actually conﬁrmed by literature by Baker et al [29], who noted that the
most widely used measure of composition uniformity is the Coefﬁcient of Variation,
(CoV), which is the ratio of the standard deviation in composition, (RMSD) and the mean

composition (AVG). These statistics are deﬁned below:
2x.-
AVG =—— ; Mean
N
- 2
DEVSQ = Z[X, — X r] ; Square of deviations

DEVSQ
N —1

RMSD = ; Root Mean Square Deviation (o)

 

74

 

I‘m-.1 r. H. .

_ RMSD

CoV ; Coefﬁcient of Variation

 

Here X ,. represents the value of volume fraction of liquid2 calculated for the i’” node in

the sample and N represents the total number of data points present in the sample of the

cross-section. CoV of 1% or less (<0.01) was given as a general rule of thumb in the

 

literature for blending of colors to uniformity. A CoV of 1%, in other words, indicates

1.11“. il..-‘

99% complete mixing.

 

‘15

Using this deﬁnition for mixing rate, values of Gov were calculated in MATLAB

[Appendix B] for each of the runs. The values are tabulated in Table 4 for a few cases.

It is seen that in the case of ribbed ducts, the CoV drops rapidly from 5% at the lowest
Reynold’s number to less than 1% at Re > 50. The CoV in fact drops to less than 0.1% at

Re > 100. This indicates an extremely efﬁcient mixing.

In the ribless ducts, the CoV is constant around 0.8. Although no mixing is expected to
take place, the CoV is not 1 because of the limited mixing at the interface of the liquids as
they move along the straight channel. Thus mixing is only about 20% complete in the

ribless duct, an extremely poor rate when compared to the ribbed ducts.

75

Table 4: Simulation results tabulated for the ribbed duct

 

 

 

 

 

 

 

 

 

 

Density of water 1000 kg/m3 Channel Height 1.6 mm

Kinematic viscosity of water 0.001kg/ms Channel Width 1.6 mm

Density of liquid A 1000 kg/m3 Hydraulic Diameter 1.6 mm

Kinematic viscosity of liquid A 0.01kg/ms Row: Reynolds No. of water

Density of liquid B 1000 kg/m3 c = Cross-Section

Kinematic viscosity of liquid B 0.03kg/ms vf2 = Volume fraction of liquid 2

Density of glycerin 1259.9 kg/m3

Kinematic viscosity of glycerin 0.799kg/ms

Average of vf2 Statistics for c 4
Mixtures 3:3 cl c2 c3 c4 DEVSQ RMSD CoV

Water-Water 400 0.604 0.5 0.5 0.5 2.4E-05 3.9E-04 7.9E-04
Water-Water 800 0.592 0.5 0.5 0.5 5.2E-06 1.8E-04 3.6E-04
Water- Liquid A 400 0.437 0.5 0.5 0.5 0.0043 0.0052 1E-02
Water- Liquid A 800 0.467 0.5 0.5 0.5 4.8E—05 5.5E-04 LIE-03
Water- Liquid B 400 0.615 0.443 0.51 0.5 0.1346 0.0291 5.8E-02
Water-Liquid B 800 0.608 0.466 0.51 0.5 0.0417 0.0162 3.2E-02
Water-Glycerin 400 0.746 0.674 0.64 0.5 0.206 0.036 7.2E-02
Water-Glycerin 800 0.738 0.65 0.61 0.5 0.0358 0.015 2.9E-02

 

 

 

 

 

 

 

 

 

 

76

 

 

Table 5: Simulation results tabulated for the ribless duct

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Average of vf2 Statistics for c 4

Mixtures 3: cl c2 c3 c4 DEVSQ RMSD CoV
Water-Water 400 0.5 0.5 0.5 0.5 26.079 0.4050 8.1E-01
Water-Water 800 0.5 0.5 0.5 0.5 26.079 0.4050 8.1E-Ol
Water—LiquidA 400 0.54 0.52 0.52 0.52 31.798 0.4472 8.6E-01
Water-LiquidA 800 0.53 0.52 0.52 0.52 31.798 0.4472 8.6E-Ol
Water-LiquidB 400 0.58 0.56 0.55 0.55 36.405 0.4785 8.7E-Ol
Water-LiquidB 800 0.57 0.56 0.55 0.55 37.2467 0.484 8.8E-01
Water-Glycerin 400 0.57 0.56 0.55 0.55 42.298 0.5170 9.4E-Ol
Water-Glycerin 800 0.55 0.55 0.54 0.54 41.4045 0.5103 94513—01

 

4.5 RESULTS AND DISCUSSION

4.5.1 Effectiveness of Microchannel System

The above analysis reveals some interesting results. The value of CoV calculated at a
cross-section indicates the amount of mixing that has taken place between the two ﬂuids
when they pass through the cross-section. In case of the ribbed duct the mixing rates were
found to be much higher than the ribless duct. Taking the ratios of CoV calculated at
cross-section 4 for the two ducts for Re of water = 800, the results indicate that the ribbd

duct produced about 100 times more mixing than the ribless duct.

77

 

 

We expect mixing in the ribless duct to decrease rapidly with increasing Reynolds
number due to corresponding decrease in residence times. Thus it has an inverse
relationship with diffusion-dominated mixing and Reynolds number [28]. However, since
diffusion is not modeled here, the mixing rate is expected to be independent of Reynolds
number. The results obtained in the simulation strongly correlates with this expectation.
The CoV values obtained for the case of water-water mixing are constant at 0.8 for all

Reynolds numbers.

4.5.2. Validation of results

Although the models have not been directly validated with any experimental or
theoretical results, their correlation of the general trends in mixing and with those of Liu
et al [28] builds conﬁdence in the simulation models. Since the criteria for design of the
micromixers were deﬁned in terms of the occurrence of chaotic advection, it was decided
to look at the conditions and indicators of chaotic mixing in simulation models as the

next step in validation.

The ﬂows around the ribs should be three-dimensional with secondary ﬂows being
generated in the channel cross-section in addition to the bulk ﬂow along the axis of the
channel. The secondary ﬂows in combination with the axial ﬂow distort and stretch
material interfaces and can produce chaotic advection which leads to rapid mixing (Liu et
al, [28]). Thus secondary ﬂows are expected to occur in the ribbed duct. The ribless duct
only allows one dimensional ﬂow 0 the ﬂuid and hence no chaotic advection is expected

to occur.

78

 

1"'r.~ r: ‘ v
. i

3‘ "ulterior 5 I: oldie-:1 Eiy ".ielcuzrtj; l-.~lr.«.r:_trntun:ie .; mi. .ture :1 rjmis. :1

FLUEI‘JT E- 1 [3d, [it], 5.91'

 

Figure 43: Velocity vectors in a cross-section in the ribbed area

Secondary ﬂows occur in the ribbed duct because the ribs stir and rotate the ﬂuid. This
was conﬁrmed by the plot of velocity vectors at a rib cross-section (Figure 43). The
ﬁgure in this thesis is presented in color. It was concluded that the rapid mixing and

strong secondary ﬂow is consistent with the occurrence of chaotic advection.

4.6 FABRICATION OF MICROCHANNELS

The same technique was used to manufacture the ribbed internal channels as described
earlier (Chapter 3). Graphite was used as the fugitive phase in the shape of desired
channels during compaction, which decomposes during sintering leaving cooling,

channels (Figures 7 and 8). Ribs were machined onto the graphite pieces using CNC

79

 

controlled machines. The idea of using a fugitive phase was really advantageous as we
could machine out the complicated rib geometry with ease using a CNC controlled
machine. Hence, we were able to successfully process internal channels of the given

geometry.

4.7 FUTURE WORK

Immediate work is proposed to validate the results of the CFD analysis by performing a
series of detailed experiments on the ribbed ducts. An arrangement could be set up such
that we have water ﬂowing into the channel through Inlet 1 and blue ink through Inlet 2.
The experiments can then be used to track particles and image along various depths of the

channel, and the results compared with the ribless ducts.

4.8 CONCLUSIONS AND EXTENSIONS

While CFD is clearly useful for illuminating the details of chaotic mixing ﬂows, a more
exciting application would be optimization. In particular, a systematic study over many
parameter values (groove depth, width, spacing, Reynolds number, etc...) could be of
great beneﬁt both to those interested in building better mixers and those interested in
understanding how and why they work. An even more appealing numerical study would
be an automated iterative parametric optimization, using Fluent as a “black box” to
calculate and return mixing length (or any other appropriate mixing metric) to an
optimization routine suited for multidimensional parameter space searching, e.g., a
Nelder-Mead simplex method. One may also try other experimentally plausible
geometries, limited only by the capabilities of current micro-machining and powder

processing techniques.

80

 

I {j _ '
m-

APPENDIX A

C++ code to modify the CFD grid that was borrowed from Dr. T. Shih

#include <fstream>
#include <cassert> '9‘“
#include <iostream>

#include <stdio.h>

 

#include <stdlib.h> t
#include <unistd.h> r
#include <string>
#include <iomanip>
#include <vector>
using namespace std;
void ﬁndminmax(vector<double> vectorarray, double &min, double &max)
{
for(int p = 0; p < vectorarray.size(); p-H-)
{
if (1) == 0)
{
min = vectorarray[p];
max = vectorarray[p];

}

81

else

if( min > vectorarray[p])
min = vectorarray[p];
if( max < vectorarray[p]) r

{

max = vectorarray[p];

r l

 

int main()
{
double i,j,k,m;
vector<double> ivalue;
vector<double> jvalue;
vector<double> kvalue;
vector<double> mvalue;
double ivaluemin,ivaluemax,jvaluemin,jvaluemax,kvaluemin,kvaluemax;
double precisioni, precisionj, precisionk;

string inputFilenameB;

82

string commentsl, commentsZ, comments3, comments4, commentsS, comments6,
comments7, comments8;

cout << "please enter the input ﬁle name: ";

 

getline(cin, inputFilenameB); // reading line

ifstream instreamB;

instreamB.open(inputFilenameB.data()); "
assert(instreamB.is_open()); ,
getline(instreamB,comments1); // reading line E

getline(instreamB,commentsZ);
getline(instreamB,comments3);
getline(instreamB,comments4);
getline(instreamB,commentsS);
getline(instreamB,comments6);
getline(instreamB,comments7);

getline(instreamB,comments8);

string Irnatch = "1:";
string Jmatch = "1:";
string Kmatch = "K=";
string comm = ".";

string::size_type posl = comments7.ﬁnd (Irnatch,0);

string::size_type posll = comments7.ﬁnd (comm,posl);

83

string::size_type pos2 = comments7.find (J match,0);
string::size_type pos2] = comments7.ﬁnd (comm,pos2);
string::size_type pos3 = comments7.ﬁnd (Kmatch,0);
string::size_type pos31 = comments7.ﬁnd (comm,pos3);
string Inumber = comments7.substr(posl+2,pos11-posl-2);
string Inumber = comments7.substr(p082+2,pos21-pos2-2);

string Knumber = comments7.substr(pos3+2,p083l-pos3-2);

int Ivalue = atoi(Inumber.c_str());
int J value = atoi(Jnumber.c__str());

int Kvalue = atoi(Knumber.c__str());

cout << "Ivalue is " << Inumber << " J value is " << I number << " Kvalue is " <<

Knumber << endl;

cout << " Total Number is " << Ivalue * I value * Kvalue << end];
for( int p = 0; p < Ivalue * J value * Kvalue ; p++)
I

instreamB >> i;
ivalue.push_back(i);
}

for( int p = 0; p < Ivalue * Jvalue * Kvalue ; p++)

84

 

I
instreamB >> j;
jvalue.push_back(j);
}

for( int p = 0; p < Ivalue * J value * Kvalue ; p++)
{ 5""
instreamB >> k; E;
kvalue.push_back(k);

l

 

for( int p = 0; p < Ivalue * I value * Kvalue ; p++)
{
instreamB >> m;

mvalue.push_back(m);

}

ﬁndminmax(ivalue,ivaluemin,ivaluemax);

ﬁndminmaxﬁvaluej valuemin,jvaluemax);
ﬁndminmax(kvalue,kvaluemin,kvaluemax);

cout << " i value Min and Max = " << setprecision(10) << ivaluemin << " " <<

ivaluemax ;

cout << " j value Min and Max = " << setprecision(10) << jvaluemin << " " <<

jvaluemax ;

85

cout << " k value Min and Max = " << setprecision(10) << kvaluemin << " " <<
kvaluemax ;
instreamB.oloseO;

// printing output ﬁle

//fﬂush(cin);

cout << "Enter the name of the Output File: ";

string outputFile;

getline(cin,outputFile);

ofstream outStream(outputFile.data());

assert(outStream.is_open());

cout << "please enter the Magniﬁcation Factor for X: ";
cin >> precisioni;
cout << "please enter the Magniﬁcation Factor for Y: ";
cin >> precisionj;
cout << "please enter the Magniﬁcation Factor for Z: ";

cin >> precisionk;

outStream << commentsl << endl; // reading line
outStream << commentsZ << endl;
outStream << comments3 << endl;
outStream << comments4 << endl;

outStream << commentsS << endl;

86

 

outStream << comments6 << endl;
outStream << comments7 << endl;

outStream << comment88 << endl;

for(int p = 0; p < Ivalue * Jvalue * Kvalue;p++)

I

outStream << scientiﬁc << setprecision(9) << (ivalue[p] - ivaluemin)/(ivaluemax -
ivaluemin)

* precisioni << " " ;

if ((p+l)%5==0)

outStream << "\n" ;

}

outStream << "\n";

for(int p = 0; p < Ivalue * J value * Kvalue;p++)

I

outStream << scientiﬁc << setprecision(9) << (ivalue[p] - jvaluemin)/(ivaluemax -
jvaluemin)

* precisionj << " " ;

if ((p+l)%5==0)

outStream << "\n" ;

}

outStream << "\n";

for(int p = O; p < Ivalue * Jvalue * Kvalue;p++)

87

 

'. Zﬂ__"l." "‘ 1' )-

I
a..-

I:

I

outStream << scientiﬁc << setprecision(9) << (kvalue[p] - kvaluemin)/(kvaluemax -
kvaluemin)

* precisionk << " " ;

if ((p+1)%5==0)

outStream << "\n" ;

}

outStream << "\n";

for(int p = 0; p < Ivalue * Jvalue * Kvalue;p++)

I

outStream << scientiﬁc << setprecision(9) << mvalue[p] << ;
if ((p+l)%5==0)
outStream << "\n" ;

}

outStreamcloseO;

88

 

APPENDIX B

Matlab code to calculate CoV

% Calculate coefﬁcient of variation from Fluent output ﬁle (node, x, y, z, vf2) for a

cross-section

clc;
clear;
[ﬁlename1, pathnamel] = uigetﬁle (‘*.dat’, ‘ Choose ﬁle’);

[node, x, y, z, vf2] = textread (fullﬁle (pathnamel, ﬁlenamel), ‘%f %f %f %f %f ’ );

sd = std (vf2)
avg = mean (vf2)
CoV = std (vf2) / mean (vf2)

DEVSQ = (length (vf2) -1) * (sd) "2

89

 

APPENDIX C

Properties of the FGM (from Mori—Tanaka scheme)

Table 1: Assuming Y-PSZ as the matrix phase and alumina as the particle phase

Y-PSZ 41.03 Thermal
Matrix Particle Young’s Poisson’s conductivity
(%) (%) Modulus (N/m’) Ratio CTE (/°C) (W/m/K)

100 0 2.00000E+1 1 0.3 1.03E-05 3

80 20 2.23300E+1 1 0.28767 9.78E-06 5.7018
60 40 2.49500E-l-1 1 0.27477 9.31E-06 8.6916
50 50 2.63684E+l 1 0.2682 9.08E-06 10.3106

40 60 2.787OOE+1 1 0.2612 8.87E-06 12.022
20 80 3.1 1693E+1 1 0.24677 8.47E-06 15.7626
0 100 3.50000E+l 1 0.23 8. 1013-06 20

 

Table II: Assuming alumina as the matrix phase and Y-PSZ as the particle phase

Y-PSZ 4120. Thermal
Particle Matrix Young’ s Poisson’s conductivity
(%) (%) Modulus (N/m’) Ratio CTE (/°C) (W/m/K)

100 0 2.00000E+1 1 0.3 1.03E-05 3
80 20 2.25000E+1 1 0.2865 9.75E-06 5.857
60 40 2.52120E+11 0.2729 9.27E-06 8.9491
50 50 2.66570E+1 1 0.2662 9.05E-06 10.5933
40 60 2.81640E+l 1 0.2593 8.84E-06 12.3091
20 80 3.13900E+1 1 0.2454 8.45E-06 15.9767
0 100 3.50000E+11 0.23 8.10E-06 20
Table 111: Final averaged values
Thermal
Y-PSZ A1203 Young’ s Poisson’s conductivity
(%) (%) Modulus (N/mz) Ratio CTE (/°C) (W/m/K)
100 0 2.00000E+1 1 0.3 1.03E-05 3.00E+00
80 20 2.24150E+l 1 0.287085 9.77E-06 5.78E+00
60 40 2.50810E+1 1 0.273835 9.29E-06 8.82E+00
50 50 2.65127E+1 1 0.2672 9.07E-06 1.05E+01
40 60 2.80170E+1 1 0.26025 8.85E-06 1.22E+01
20 80 3.12796E+1 1 0.246085 8.46E-06 1.59E+01
0 100 3.50000E+11 0.23 8.10E-06 2.00E+01

90

REFERENCES

1). Benson, R. 8., and J.W. Ponton, “Process Miniaturization— A Route to Total
Environmental Acceptability?” Trans. IChemE, (V71, Part A, 1993), pp160-168.

2). Koeneman, P. B., I. J. Busch-Vishniac, and K. L. Wood, 1997, “Feasibility of Micro
Power Supplies for MEMS,” J. MicroElectoMechanical Sys. (V6, n4,), pp355-362

3). Paul, B. K., and RB. Peterson, "Microlamination for Microtechnology-based Energy,
Chemical, and Biological Systems," ASME IMECE, Nashville, Tennessee, November
15-20, 1999, pp45-52.

4). Martin, P. M., W. D. Bennett, and J. W. Johnson, “Microchannel Heat Exchangers for
Advanced Climate Control,” Proc. SPIE (v2639, 1996), pp82-88

5). Matson, D. W., P. M. Martin, W. D. Bennett, D. C. Stewart, and J. W. Johnson,
“Laser Micromachined Microchannel Solvent Seperator,” Proc. SPIE, (v3223, 1997),
pp253-259.

6). Matson, D. W., P. M. Martin, A. Y. Tonkovich and G. L. Roberts, “Fabrication of a
Stainless Steel Microchannel Microcombustor Using a Lamination Process,” Proc. SPIE
(V3514, 1998). Pp386-392.

7). Tonkovich, A.Y.; Zilka, J.L.; LaMont, M.J.; Wang, Y. and Wegeng, R.S.
“Microchannel reactors for Fuel Processing Applications,” Chemical Engineering
Science v 54 n 13 Sep l3-Sep 16 1998 1999 p 2947-2951 0009-2509

8). Paul, B.K. and T. Terhaar, “Comparison of two passive microvalve designs for
rnicrolarnination architectures,” J Micromech. Microengr. (VIC, 2000), pp15-20.

9). Martin, P. M., D. W. Matson, W. D. Bennett, and D. J. Hammerstrom, “Fabrication of
plastic rnicroﬂuidic components,” Proc. SPIE (V3515, 1998), pp172-176

10). Wegeng, R. S., C. J. Call, and M. K. Drost, “Chemical Systems Miniaturization,"
1996 AIChE Spring Meeting, New Orleans, 1996.

ll). Kakac, S. and H. Liu, 1997, Heat Exchangers-Selection, Rating, and Thermal
Design, CRC Press, Boca Raton.

12). Kanlayasiri, K., and BK. Paul, 2004, “A Nickel Aluminide Microchannel Array
Heat Exchanger for Hi gh-Temperature Applications,” accepted by J. Mfg Proc..

13). Brooks, K. P., P. M. Martin, M. K. Drost, and C. J. Call, “Mesoscale Combustor/

Evaoprator Development,” ASME IMECE Conference, Nashville, TN, November 1999,
pp37-43.

91

 

1F
1

14). Wang, J. H., Yeh, H. Y. and Shyu, R. J ., 1991, “Thermal-Hydraulic Characteristic of
Microheat Exchangers,” American Society of Mechanical Engineers, Dynamic Systems
and Control Division (Publication) (DSC v32 Dec 1-6 1991), pp33l-339.

15). Peng, X. F., B. X. Wang, G. P. Peterson, and H. B. Ma, “ Experimental Investigation
of Heat Transfer in Flat Plates with Rectangular Microchannels,” International Journal of
Heat and Mass Transfer (v38, nl, 1995), pp127-137.

l6). Friedrich, C. R. and S. D. Kang, 1994, “Microheat Exchangers Fabricated by
Diamond Machining,” Precision Engineering, 16, 1, pp56-59.

17). Paul, B. K., R. B. Peterson, and W. Wattanutchariya, “The Effect of Shape Variation
on the Performance of High-Aspect-Ratio, Metal Microchannel Arrays,” IMRET 3:
Proceedings of the Third International Conference on Microreaction Technology,
Springer-New York, 1999, pp53-6l.

18). Jovanovic, G., J. Zaworski, T. Plant, and B. Paul, Proc. 3rd Int]. Conf. on
Microreaction Technology, (Frankfurt, Germany), 1999.

19). Little, W. A., “Microminiature Refrigerators for Joule-Thomson Cooling of
Electronic Chips and Devices,” Advances in Cryogenic Engineering, (v33, 1990), pp.
1325-1333.

20). Paul, B.K., H. Hasan, J. Thomas, R. Wilson, and D. Alman, “Limits on Aspect Ratio
in Two-ﬂuid Micro-scale Heat Exchangers,” Transactions of NAMRC XXIX, Gainesville,
FL, 2001.

21). Gray, K.J., 2000, “Diamond and Related Materials”, 9: (2), pp. 201-204

22). Jagannadham, K., 1999, J. Vac. Sci. Technol, A17: (2), pp. 373-379

23). Hui, P., Tan, HS. and Lye, Y.S., 1997, IEEE Tran. Comp. Pack. & Manuf. Tech.
Part A, 20: (4); pp. 452-457

24). Bower, MB. and Mudawar, I., 1994, J .Electron Package, 116, pp. 290-297.

25). Gillot, C., Meysenc L., Schaeffer, C. and Bricard, A., 1999, IEEE Trans. Comp.
Pack. Tech., 22: (3), pp. 384-389

26). Mizuno, Y., Kawasaki, A., and Watanabe, R., 1995, Metal]. Mater. Trans, 26 B,
p.75.

27). Zhang, L.M., Yuan, R.Z., Oomori, M. and Hirai, T., 1995, J. Mat. Sci. Let., 14, 22,
p.1620.

28). Liu et al., “Passive Mixing in a Three-Dimensional Serpentine Microchannel”,
Journal of Microelectromechanical Systems, Vol. 9, No. 2, June 2000, pp. 190-197.

92

 

7W1“ . Qm‘l.-L\
l0

29). Baker et al., “Laminar Flow in Static Mixers with Helical Elements”, The Online
CFM Book, http://bakerorg/cfm, pp. 1-11

 

30). Shih, T.I.P. et al., “A Numerical Study of Flow and Heat Transfer in a Smooth and
Ribbed U-Duct With and Without Rotation”, ASME Journal of Hear Transfer, Vol. 123,
April 2001, pp. 219-232

31).German, R.M., 1984, Powder Metallurgy Science, Metal Powder Industries
Federation, Princeton, New Jersey.

32). Yeh, T. and Sacks, M.D., 1998, J. Am. Ceram. Soc., 71, 12, pp. C484-V487

33). Mori, T., and Tanaka, K., 1973, “Average Stress in Matrix and Average Elastic

Energy of Materials with Misfitting Inclusions,” Acta Metallurgica et Materialia, Vol. 21,
pp. 571-574

34). Benveniste, Y., 1987, “A New Approach to the Application of Mori-Tanaka’s
Theory in Composite Materials,” Mechanics of Materials, Vol.6, pp. 147-157.

35). Shin, H.W., Kwon, P. and Case, E.D., ASM Fall 2000 Proceedings.

36). Roeder, B.A., Sun, C.T., “Dynamic penetration of alumina/aluminum laminates:
experiments and modeling”, International Journal of Impact Engineering 25 (2001), pp.
169-185.

37). Stroock, A.D., Dertinger, S.K.W., Ajdari, A., Mezic, 1., Stone, H.A., and Whitesides,
G.M., “Chaotic mixer for microchannels.” Science Magazine, 295, 2002.

38). Stroock, A.D., Dertinger S.K.W., Whitesides G.M., Ajdari, A. “Patterning ﬂows
using grooved surfaces.” Analytical Chemistry, 74 (20): 5306-5312, 2002

39) Paul, Brian, Subramanian, Ramkumar and Kwon, Patrick “Understanding limits on
ﬁn aspect ratios in two-ﬂuid micro-channel arrays produced by diffusion bonding”,
submitted to the ASME Journal of Manufacturing Science and Engineering, special issue
on micro-mechanical manufacturing, 2004

93

 

 

   

llililljjlliilliijlllililjil