Feasibility Study for a Novel Micropump Design That Utilizes Shape Memory Alloy
Diaphragm and Blood Flow
by
Mark Anthony Thrasher
A dissertation submitted to the Graduate Faculty of
Auburn University
in partial fulfillment of the
requirements for the Degree of
Doctor of Philosophy
Auburn, Alabama
August 9, 2010
Keywords: micropump, shape memory alloy, blood, lattice Boltzmann
Copyright 2010 by Mark Anthony Thrasher
Approved by
Daniel W. Mackowski, Co-chair, Professor of Mechanical Engineering
David F. Dyer, Co-chair, Professor of Mechanical Engineering
George T. Flowers, Professor of Mechanical Engineering
Amnon J. Meir, Professor of Mathematics and Statistics
ii
Abstract
Many different types of micropumps have been designed to pump a variety of
fluids at the micro-liter level. Some of the major short-comings of these designs are leaks
due to moving parts, lack of tolerance to bubbles, high operating voltages, structural
complexity, slow response time, and insufficient pressure head. This research presents a
novel design that utilizes a two-way shape memory alloy (SMA) diaphragm within the
micropump. Heating of the SMA diaphragm is accomplished via resistive heating, while
the cooling is accomplished with a novel method that utilizes blood flow. Blood is used
as the cooling fluid since it is readily available in most medical applications, such as drug
delivery or lab-on-a-chip (LOC). Additionally, with humans there is a constant supply of
blood at a constant temperature.
This study explains and illustrates how an SMA micropump operates. With the
primary liquid flow in this application being blood, a detailed model of the non-
Newtonian behaving blood is included. The hydrodynamics are numerically computed
using the Lattice Boltzmann Method (LBM). Two different non-Newtonian blood
models are presented. Phase changes in the shape memory alloy are computed using the
Cosine Model. All of the thermal analysis is performed via finite volumes.
Ultimately, this study shows that an SMA micropump is feasible when pumping
fluids in the micro-liter range of flow. This work is a feasibility study and addresses
iii
operation parameters in a general sense. All of the design parameters are reasonable for a
micro-liter micropump.
iv
Acknowledgments
The author would like to thank Dr. George Flowers for his committee assistance
and Dr. Daniel Mackowski and Dr. David Dyer for taking over as my chair during a
rather tumultuous time in my research. I am very thankful for that. Also, thanks to Dr.
Amnon Meir for all the knowledge I gained while taking his classes. I want to thank my
wife who went through the PhD program the same time as me. Lastly, thanks to my two
kids, Jessica and Brandon, for tolerating their parents during this stage of our lives.
v
Table of Contents
Abstract ............................................................................................................................... ii
Acknowledgments.............................................................................................................. iv
List of Tables ..................................................................................................................... ix
List of Figures ..................................................................................................................... x
Chapter 1 Introduction ........................................................................................................ 1
1.1 Research Motivation ................................................................................................. 1
1.2 SMA Micropump Operation ..................................................................................... 2
1.3 Expected Research Contributions ............................................................................. 3
1.4 Dissertation Organization ......................................................................................... 3
Chapter 2 Literature Review ............................................................................................... 5
2.1 Micropumps .............................................................................................................. 5
2.1.1 Bubble Tolerance and Self-Priming .................................................................. 8
2.1.2 Diaphragm Deflection Analysis ....................................................................... 10
2.1.3 Micropump Summary ...................................................................................... 12
2.2 Shape Memory Alloys ............................................................................................ 13
2.2.1 Basic Operation Overview of Shape Memory Alloys ..................................... 14
vi
2.2.2 Shape Memory Alloy Phase Change Temperatures ........................................ 16
2.2.3 SMA Transformation Kinetics ......................................................................... 18
2.2.4 SMA Mechanical Properties ............................................................................ 20
2.2.5 SMA Physical Properties ................................................................................. 23
2.2.6 SMA Chemical Properties ............................................................................... 24
2.2.7 SMA Biocompatibility ..................................................................................... 25
2.3 Lattice Boltzmann Method ..................................................................................... 27
2.3.1 General Overview ............................................................................................ 28
2.3.2 Boltzmann Equation ........................................................................................ 30
2.3.3 Lattice Boltzmann BGK .................................................................................. 31
2.3.4 Stability of Lattice Boltzmann Method ........................................................... 35
Chapter 3 Shape Memory Alloy Diaphragm Model ......................................................... 47
3.1 Micropump Overview ............................................................................................. 47
3.2 SMA Micropump Diaphragm ................................................................................. 48
3.2.1 Phase Transformation Model for SMA Diaphragm ........................................ 48
3.2.2 Thermal Model for SMA Diaphragm .............................................................. 49
3.2.3 SMA Diaphragm Deformation Geometry ....................................................... 50
3.2.4 SMA Diaphragm Kinematics ........................................................................... 52
3.3 Numerical Analysis of SMA Diaphragm ................................................................ 55
3.3.1 Transient Displacement Analysis .................................................................... 55
vii
3.3.2 Transient Velocity Analysis ............................................................................. 56
Chapter 4 Lattice Boltzmann Boundary Conditions ......................................................... 67
4.1 Boundaries Overview.............................................................................................. 67
4.2 Velocity and Pressure Boundaries using Lattice Boltzmann Method .................... 67
4.2.1 Velocity Boundary Conditions ........................................................................ 68
4.2.2 Pressure Boundary Conditions ......................................................................... 73
4.3 Moving Walls and Curved Boundaries for SMA Diaphragm ................................ 75
Chapter 5 Hydrodynamics of Micropump ........................................................................ 77
5.1 Hydrodynamic Overview ........................................................................................ 77
5.2 Lattice Boltzmann Benchmarking .......................................................................... 78
5.2.1 Units Example .................................................................................................. 78
5.2.2 Womersley Flow for Benchmarking ................................................................ 84
5.2.3 Convergence Analysis (Grid Independence) ................................................... 89
5.3 Blood Characterization ........................................................................................... 89
5.3.1 General Properties and Composition of Whole Blood .................................... 89
5.3.2 Non-Newtonian Nature of Blood Viscosity ..................................................... 90
5.4 Boundary Implementation ...................................................................................... 94
5.4.1 Inlet Velocity ................................................................................................... 94
5.4.2 Outlet Pressure ................................................................................................. 96
5.4.3 All Other non-Fluid Boundaries ...................................................................... 96
viii
5.5 Numerical Results ................................................................................................... 96
5.5.1 Blood Flow Velocity ........................................................................................ 97
5.5.2 Viscosity .......................................................................................................... 98
5.5.3 Hydrodynamic Boundary Layer on SMA Diaphragm ..................................... 99
Chapter 6 Thermal Results of Micropump ..................................................................... 132
6.1 Thermal Overview ................................................................................................ 132
6.2 General Approach for Thermal Solver ................................................................. 132
6.3 Introduction of Finite Volume Method for Thermal Analysis .............................. 133
6.4 Energy Analysis of SMA Diaphragm ................................................................... 136
6.5 Thermal Flux Analysis for Cooling of SMA Diaphragm ..................................... 138
6.5.1 Thermal Flux Analysis with Conduction Only (SMA to Blood) ................... 139
6.5.2 Thermal Flux Analysis with Conduction and Convection (SMA to Blood) . 140
6.6 Temperature Profile along SMA Diaphragm ........................................................ 140
6.7 Multicycle Thermal Analysis of Micropump ....................................................... 142
6.8 Micropump Analysis Based upon Thermal Modeling .......................................... 142
Chapter 7 Conclusions .................................................................................................... 153
Bibliography ................................................................................................................... 155
ix
List of Tables
Table 2.1 Classification and Operation of Micropumps ..................................................... 8
Table 2.2 Critical Pressure Evaluation.............................................................................. 10
Table 2.3 Properties of Diaphragm ................................................................................... 11
Table 2.4 Mechanical Data for Trained and Untrained Nitinol SMA .............................. 21
Table 2.5 Physical Properties of Nitinol SMA ................................................................. 25
Table 5.1 Parameters for Womersley Flow Solution Listed in Lattice Units ................... 88
Table 5.2 CY Blood Model Parameters ............................................................................ 94
Table 5.3 Curve-Fit Equations for Blood Velocity Profile ............................................... 95
Table 6.1 SMA Transformation Temperatures ............................................................... 137
Table 6.2 Energy Requirements during Cool-Down of SMA Diaphragm ..................... 138
Table 6.3 Minimum and Maximum SMA Micropump Capacity for This Design ......... 143
x
List of Figures
Figure 1.1 Schematic of Shape Memory Alloy Micropump ............................................... 4
Figure 2.1 Breakdown of the Different Types of Pumps .................................................. 37
Figure 2.2 Illustration of SMA Wire Recovering Original Shape .................................... 38
Figure 2.3 Crystallographic view of Shape Memory Effect ............................................. 39
Figure 2.4 SMA Composition as a Function of Temperature ........................................... 40
Figure 2.5 Dependence of Transition Temperatures upon Applied Stress ....................... 41
Figure 2.6 Shape Memory Alloy Stress-Strain Curve to Ultimate Strength..................... 42
Figure 2.7 Stress-Strain Curves for Trained and Untrained Nitinol SMA ....................... 43
Figure 2.8 Modulus of Elasticity as a Function of Temperature ...................................... 44
Figure 2.9 Work Extraction from SMA Actuator ............................................................. 45
Figure 2.10 Detailed Breakdown of SMA Actuator Cycle ............................................... 46
Figure 3.1 SMA Thermal Response throughout the Phase Transformation (Typical results
with an SMA mass of 0.00258 kg) ................................................................................... 57
Figure 3.2 Basic Geometry Variables for SMA Diaphragm............................................. 58
Figure 3.3 Change in SMA Geometry after One Time Step (parameters are not to scale)
........................................................................................................................................... 59
xi
Figure 3.4 Further SMA Geometric Analysis ................................................................... 60
Figure 3.5 SMA Diaphragm Displacement Analysis ....................................................... 61
Figure 3.6 Motion of Diaphragm from Maximum Deflection Point at the Center of
Diaphragm. (The larger the negative number, the greater the penetration into the cavity.)
........................................................................................................................................... 62
Figure 3.7 X Velocity along the Width of the SMA Diaphragm ...................................... 63
Figure 3.8 Maximum X Velocity of Diaphragm at each Time Step ................................. 64
Figure 3.9 Maximum Y Velocity along Width of SMA Diaphragm ................................ 65
Figure 3.10 Maximum Y Velocity of Diaphragm at each Time Step ............................... 66
Figure 4.1 Different Boundaries for LBM ........................................................................ 76
Figure 5.1 Schematic of Blood Micropump with Penetrating SMA Diaphragm ........... 100
Figure 5.2 Comparison of Womersley Flow (Diamonds-Fourier Series, Solid Line-LBM)
......................................................................................................................................... 101
Figure 5.3 Lattice Boltzmann Relative Velocity Error for Varying Square Grids ......... 102
Figure 5.4 Fourier Series Relative Velocity Error for Varying Square Grids ................ 103
Figure 5.5 Cross-sections at Center of First Chamber of Micropump ............................ 104
Figure 5.6 Cross-sections at Center of Last Chamber of Micropump ............................ 105
Figure 5.7 Rheological Behavior of Whole Blood ......................................................... 106
Figure 5.8 Strain Rate for 2-D Channel Flow ................................................................. 107
Figure 5.9 Viscosity Profile for Shear-Thickening Fluid ............................................... 108
xii
Figure 5.10 Viscosity Profile for Shear-Thinning Fluid ................................................. 109
Figure 5.11 Velocity Profile for Shear Thickening Fluid ............................................... 110
Figure 5.12 Velocity Profile for Shear-Thinning Fluid .................................................. 111
Figure 5.13 Velocity Profiles Normalized for Maximum Velocity ................................ 112
Figure 5.14 Viscosity Predictions from CY Blood Model ............................................. 113
Figure 5.15 Typical Velocity Profile for Pumped Blood ................................................ 114
Figure 5.16 Axial Velocity (lattice units) throughout Blood Micropump ...................... 115
Figure 5.17 Transverse Velocity (lattice units) throughout Blood Micropump ............. 116
Figure 5.18 Streamlines of Blood Flow Through Micropump with SMA at 2% Strain . 117
Figure 5.19 Streamlines of Blood Flow Through Micropump with SMA at 8% Strain . 118
Figure 5.20 Axial Velocity of Blood for Pump Cross-Section at Diaphragm Center .... 119
Figure 5.21 Transverse Velocity of Blood for Pump Cross-Section at Diaphragm Center
......................................................................................................................................... 120
Figure 5.22 Maximum Blood Velocity at top of Ramp and Exit of Micropump ........... 121
Figure 5.23 Velocity Profiles at end of One Second of Blood Flow .............................. 122
Figure 5.24 Strain Rate of Blood in Micropump ............................................................ 123
Figure 5.25 Viscosity of Blood Based upon CY Model ................................................. 124
Figure 5.26 Modified Relaxation Factor (Tau) throughout Micropump ........................ 125
Figure 5.27 Axial Velocity Boundary Layer at First Quarter of SMA Length ............... 126
Figure 5.29 Axial Velocity Boundary Layer at Midpoint of SMA Length .................... 128
xiii
Figure 5.30 Transverse Velocity Boundary Layer at Midpoint of SMA Length............ 129
Figure 5.31 Axial Velocity Boundary Layer at 3rd Quarter of SMA Length .................. 130
Figure 5.32 Transverse Velocity Boundary Layer at 3rd Quarter of SMA Length ......... 131
Figure 6.1 Space Discretization for Finite Volumes ....................................................... 144
Figure 6.2 Cardinal Directions for Implementation of Finite Volume Method .............. 145
Figure 6.3 Heat Flux from SMA (I,j) to Blood (I,j-1) .................................................... 146
Figure 6.4 Comparison of Energy Flux between SMA and Blood ................................. 147
Figure 6.5 Temperature along SMA Diaphragm during Cooling ................................... 148
Figure 6.6 Temperature along SMA Diaphragm during Cooling .................................. 149
Figure 6.7 Temperature throughout SMA at Varying Deformation ............................... 150
Figure 6.8 Thermal Cycles with/without Blood Flow .................................................... 151
Figure 6.9 Thermal Cycles using Complete Model of Heart Beat ................................. 152
1
Chapter 1 Introduction
1.1 Research Motivation
Currently, there is a need for high precision pumping (few ?l/stroke) at low flow
rates (0-500 ml/min). One such application for these micropumps is drug delivery in
which a very sensitive therapeutic efficacy is maintained. Some drugs have therapeutic
levels that require precise control and this may not be achieved with traditional devices
such as syringes. For drug delivery, micropumps may be placed subcutaneously;
therefore, they need to be biocompatible. Another application is labs-on-a-chip (LOC)
where ?liters of fluid are pumped through intricate medical testing stations on
miniaturized chips. As an example, one such medical application is a lab-on-a-chip
(LOC) that allows soldiers to do field testing for exposure to different chemicals.
Currently, there are micropumps based on many different methods of operation.
Some of the most prevalent methods are piezoelectric, electrostatic, thermopneumatic,
magnetic, electrohydrodynamic, and magnetohydrodynamic. In Chapter 2, the
advantages and disadvantages of these pumps are addressed. One of the newest actuators
considered for micropumps is a shape memory alloy (SMA) diaphragm pump. By
alternate heating and cooling of the diaphragm, precise micro-quantities of fluid at high
pressures are achievable. Additionally, SMA diaphragm pumps are not as prone to leaks
as many of the mechanical micropumps with sliding or rotating parts.
2
In the very recent years, there have been a few proposed designs of micropumps
utilizing SMA diaphragms. Some of the designs rely upon bias springs while others use
opposing diaphragms that are heated and cooled alternately. When SMA?s are used in
actuators, the slow cooling time is the limiting factor in the actuator pumping frequency.
However, none of the current SMA micropumps address different ways to optimize the
cooling phase of the cycle. After a literature search, there are no analytical/numerical
models of SMA micropumps.
1.2 SMA Micropump Operation
The SMA micropump operates by alternately heating and cooling the SMA
diaphragm. This heating and cooling allows the shape ?memory? of the diaphragm to
vary from 1% to 8% strain. The SMA diaphragm separates the medical fluid in the top
chamber from the blood flow in the lower chamber, as shown in Figure 1.1. The
alternate heating and cooling of the SMA diaphragm pushes a prescribed amount of
medical fluid out of the upper chamber with each cycle of operation. Simple check
valves made from a pliable rubber material ensure that the medical fluid (in the upper
chamber) progresses from left to right as the SMA diaphragm contracts and expands.
This research looks at a novel design of an SMA diaphragm pump that maximizes
the contact of a cooling fluid (blood) with the diaphragm. Any type of cooling fluid may
be used, but blood is considered here since medical applications and/or testing are of
interest. The micropump uses blood as a coolant, and consequently, enables precision
pumping of any medical fluid.
3
1.3 Expected Research Contributions
The expected contributions of this research are three-fold. First, this research
proposes a new design of SMA micropump that utilizes separate chambers for a pumped
fluid and the cooling fluid. Second, blood flow pumped by the human heart is used to
enhance the cooling of the SMA diaphragm, which has not been previously investigated.
Third, the SMA micropump has a unique design that overcomes many weaknesses of
popular pumps currently available. Some of these advantages are low voltage
requirements, simple design, biocompatibility, high pressure head, and high compression
ratio (bubble tolerance).
1.4 Dissertation Organization
A literature search is done in Chapter 2 for micropumps, SMAs, and the Lattice
Boltzmann method (LBM). This new micropump design is shown in greater detail in
Chapter 3. The detailed implementation of the Lattice Boltzmann boundaries is
presented in Chapter 4. Afterward, the numerical modeling for the fluid flow of blood is
discussed in Chapter 5. Finally, a thermal analysis for the predicted flow field and SMA
diaphragm is shown in Chapter 6.
Figure 1.1 Schematic of Shape Memory Alloy Micropump
4
5
Chapter 2 Literature Review
2.1 Micropumps
Microelectromechanical Systems (MEMS) are miniaturized devices that are
capable of actuating, sensing, and/or controlling systems (Shoji and Esashi 1994;
Abdelgawad, Hassan et al. 2005). Typically, they involve the integration of mechanical
and electrical technologies. MEMS were introduced in the 1960?s; specifically, a gold
resonating MOS gate structure (Nathanson 1967). Since the advent of MEMS, the desire for
further capabilities and miniaturization continues to increase.
Of particular interests in this ongoing research are the study of microfluidics and
the design of micropumps that can reside within MEMS or operate independently. Flow
rates for macroscale pumps are about 100 ml/min or greater (Nguyen, Huang et al. 2002).
The need for micropumps is when the volume of liquid pumped can be as low as a few
?l/min.
A particular area needing optimized micropumps is the medical field. Many
medical applications, commonly referred to as BioMEMS, (Grayson, Shawgo et al. 2004)
require micropumps that are capable of delivering small and precise quantities of
fluid/drugs (Maillefer, Gamper et al. 2001; Shawgo, Grayson et al. 2002; Voskerician,
Shive et al. 2003; Grayson, Shawgo et al. 2004; Kan, Yang et al. 2005; Yoshida 2005).
These accurate delivery quantities have a direct impact on the efficacy of the drugs
6
(Bakken and Heruth 1991). Typical medical applications using micropumps are listed as
follows:(Shawgo, Grayson et al. 2002):
? Micro-quantity drug delivery
? Milliliter-scaled fluids for lab-on-a-chip analysis(Ikuta, Takahashi et al. 2003;
Erickson 2005)
? Assisting damaged heart valves(Stepan, Levi et al. 2005)
Depending upon the application, there are several advantageous traits that may be
required for a micropump. Some of the obvious quantifiable characteristics are the
displacement volume, operating frequency, and pressure head. The combination of the
displacement volume and operating frequency determine the maximum achievable flow
rate possible. Other secondary issues that may be of importance are the ability of a
micropump to pump fluid with high precision and reliability. Also, the weight to power
ratio of the micropump may be of significance.
There are many different types of micropumps and many different modes of
operation. Modern microfluidics is traced back to a gas chromatograph at Stanford
University and an ink-jet printer (Erickson and Li 2004). Looking at micropumps from a
very general overview, there are two different types of micropumps, namely,
displacement and dynamic. Dynamic pumps displace fluid using methods such as
centrifugal, electrohydrodynamic, electroosmotic (Erickson and Li 2003), and
magnetohydrodynamic. Displacement pumps use diaphragms and pistons in combination
with valves, chambers, etc. (Laser and Santiago 2004). Some of the most prevalent
diaphragm materials are silicon, metal, and plastic. Some of the common drivers for the
diaphragm are piezoelectric, thermopneumatic, electrostatic, and pneumatic.
7
Micropumps with sliding or rotating parts are more prone to leaking than micropumps
with deformable diaphragms. A breakdown of the different types of pumps is shown in
Figure 2.1.
The type of micropump of particular interest is one that utilizes shape memory
alloy diaphragms (Benard, Kahn et al. 1997; Benard, Kahn et al. 1998; Kahn, Huff et al.
1998). These diaphragms are constructed by applying SMA to a rubber diaphragm using
a sputtering deposit method (Makino, Mitsuya et al. 2000; He, Huang et al. 2004).
Section 2.2 delves into detail regarding shape memory alloys. There have been some
experimental testing of SMA diaphragms, although quite limited (Makino, Mitsuya et al.
2000).
A summary of the most popular pumps is shown in Table 2.1. Their mode of
operation and advantages/disadvantages are briefly listed.
Pump Category Mode of Operation Advantages (+) Disadvantages (-)
Displacement Pumps
? Electrostatic Two opposite
electrostatic plates have
voltages alternately
applied.
+ generates appreciable forces
+ operates at high frequencies
- requires high voltage
- complex structure
? Piezoelectric A piezoelectric material
is deposited onto a
membrane and a high
voltage is applied. The
quartz material
mechanically deforms.
+ operates at high frequencies
+ good efficiency
- requires high voltage
- small displacement stroke
? Thermo-pneumatic A trapped cavity of air is
heated with a resistive
heater. The
+ large displacement stroke
- complex structure
- slow response time
8
expanding/contracting air
results in the diaphragm
moving.
- low efficiency
? Bimetallic Two dissimilar metals are
joined and alternately
heated.
+ simple design
+ large forces generated
- slow response time
- small displacement stroke
? Shape Memory Alloy Nitinol SMA is
alternately heated and
cooled to induce memory
effect.
+simple design
+large displacement stroke
+large forces generated
- low efficiency
Dynamic Pumps
? Magnetohydrodynamic Uses electric and
magnetic fields to drive a
conductive fluid.
+ simple design
+ flow can be bi-directional
- bubbles from ionization
- miniaturization difficult
- joule heating
? Electroosmotic Operates by applied
potential across a
microchannel
+ generate high pressures
+no mechanical actuator
- requires high voltage
- joule heating
? Chemical Generation of bubbles
from electrolysis
+ simple design
- bubbles collapse/reliability
- low flow rate
Table 2.1 Classification and Operation of Micropumps
2.1.1 Bubble Tolerance and Self-Priming
Bubble tolerance refers to the ability of a micropump to continue working when a
bubble is introduced into the pumped fluid. Since the bubbles are gas, they can be
9
compressed to the point that the pressure in the pump is not sufficient to actuate the check
valves. Hence, the volume stroke needs to be large enough to overcome the fluid
bubbles.
The compression ratio for a micropump is defined as the ratio of the stroke
volume to dead volume. This equation is (Richter)
g2013 = ?g3023g3023
g3116
(2.1)
where ?g1848 is the stroke volume and g1848g2868 is the dead volume. The larger the compression
ratio, the more tolerant the micropump is to bubbles. For the scenario of a liquid pump
with no bubbles, the minimum compression ratio for operation is (Richter, Linnemann et
al. 1998)
g2013g3039 > g2018uni007C?g1868g3030g3045g3036g3047uni007C (2.2)
Where ?g1868g3030g3045g3036g3047 is the critical pressure needed to actuate the check valves. The
compressibility (g2018) of water is 0.5 g1876 10?8 g1865^2/g1840. For an ideal gas, the minimum
compression ratio to operate a gas pump is (Richter)
g2013g3034g3028g3046 > g2869g3082 ?g3043g3278g3293g3284g3295g3043
g3116
(2.3)
where g2011 is the adiabatic coefficient and g1868g2868 is atmospheric pressure. If there is back
pressure, the critical pressure to operate the check valves is increased by this amount.
10
Experimental work has been performed for a typical micropump with passive
silicon check valves. The check valve studied is a cantilever style with the dimensions of
1700 ?m x 1000 ?m and a thickness of 15 ?m. The critical pressure required to open the
valves is measured for two scenarios, specifically, a full liquid reservoir with no bubbles
and an empty liquid reservoir. These critical pressures are listed in Table 2.2.
Additionally, the minimum compression ratios are computed from equations 2.2 and 2.3,
and they are listed in Table 2.2.
Liquid Micropump Liquid Micropump (self priming/ bubble tolerant)
Critical Pressure
(open valves)
10 hPa
75 hPa
Minimum
Compression Ratio
? > 0.000005 ? > 0.075
Table 2.2 Critical Pressure Evaluation
2.1.2 Diaphragm Deflection Analysis
The shape memory effect of Nitinol allows for strains as high as 8% during
operation. In the case of the SMA micropump, the maximum SMA deflection for the 2
cm diaphragm is approximately 0.4 cm. This can be compared to other diaphragm
materials (such as silicon, glass, or a malleable metal) by focusing on stress limits and
driver forces. Assuming a circular diameter diaphragm, the driver force, g1868g3028, can be
calculated from this equation (Young and Budynas):
11
g3043g3276g3031g3279g3120g2869g2874g3006
g3300g3047g3279g3120
= g2873.g2871g2871(g2869g2879g3092g3118)g3052g3116g3047
g3279
+ g2870.g2874(g2869g2879g3092g3118)g4672g3052g3116g3047
g3279
g4673g2871 (2.4)
where
g1868g3028 = g1873g1866g1861g1858g1867g1870g1865 g1853g1868g1868g1864g1861g1857g1856 g1856g1870g1861g1874g1857g1870 g1858g1867g1870g1855g1857 g1868g1857g1870 g1873g1866g1861g1872 g1876 ?g1871g1857g1855g1872g1861g1867g1866g1853g1864 g1853g1870g1857g1853
g1856g3031 = g1856g1861g1853g1868?g1870g1853g1859g1865 g1856g1861g1853g1865g1857g1872g1857g1870
g1831g3052 = g1877g1867g1873g1866g1859?g1871 g1865g1867g1856g1873g1864g1873g1871
g1872g3031 = g1856g1861g1853g1868?g1870g1853g1859g1865 g1872?g1861g1855g1863g1866g1857g1871g1871
g1877g2868 = g1855g1857g1866g1872g1857g1870g1864g1861g1866g1857 g1856g1861g1871g1868g1864g1853g1855g1857g1865g1857g1866g1872
? = poisson?s ratio
Assuming a diaphragm is silicon, the typical mechanical properties are shown in Table
2.3. The driver force can be readily computed.
Variables Values (for silicon)
g2186g2186 2 cm
g2161g2207 150,000 M-Pa
g2202g2186 0.004 cm
g2207g2777 0.4 cm
g2247 0.17
Table 2.3 Properties of Diaphragm
12
Solving eq. 2.4 with the values in Table 2.3, the driver force per unit area is 15000
pounds. The internal stress in the diaphragm may also be computed using eq. 2.5 (Young
and Budynas)
g3097g3031g3279g3118g2872g3006
g3300g3047g3279g3118
= g2872(g2869g2879g3092g3118)g3052g3116g3047
g3279
+ g2869.g2875g2871(g2869g2879g3092g3118)g4672g3052g3116g3047
g3279
g4673g2870 (2.5)
The computation shows that the stress in the diaphragm would be over 6,400,000 psi.
Based on these computations, diaphragms like silicon and glass will not allow stroke
volumes comparable to a shape memory diaphragm of the same size. The only type of
diaphragm that will allow high deformations like the SMA is a Mylar or silicone rubber
low-modulus material.
2.1.3 Micropump Summary
In general, electrostatic and piezoelectric micropumps have a very small
compression ratio due to their very small stroke volume. In order to increase the
compression ratio to acceptable values, the reservoirs for them have to be very small and
micromachined. Typical electrostatic and piezoelectric micropumps used with the valves
analyzed above have compression ratios of 0.002 and 0.017, respectively (Richter,
Linnemann et al. 1998). The SMA micropump has a high compression ratio in
comparison to these pumps due to the large volume stroke of the SMA diaphragm. The
compression ratio for the SMA micropump analyzed in this dissertation is 0.24 . This
provides very high bubble tolerance in comparison to other pumps with small stroke
volumes. The large compression also permits self-priming for an SMA micropump.
13
The other major limiting factor for micropumps (listed in Table 2.1) is the
requirement of high voltages. One of the most popular pumps, piezoelectric, suffers from
the high voltage requirement. This makes subcutaneous applications not possible.
2.2 Shape Memory Alloys
Shape Memory Alloys (SMAs) are adaptive materials that have the capability of
converting thermal energy directly into mechanical work via a material phase
transformation(Perkins 1975). This phase transformation is based upon the alteration of
the crystalline structure of select alloys that have been properly mixed and heat-treated.
A listing of some of the primary elements found in many SMAs that exhibit the shape
memory alloy effect includes Copper, Zinc, Aluminum, Nickel, Cadmium Gallium and
Titanium(Wayman 1980). The most prevalent SMA by far is Nitinol (Nickel Titanium
Naval Ordnance Laboratory).
Nitinol is essentially a 50-50 stoichiometric blend of Nickel and Titanium, and
small variations in the blend ratio result in vastly different operating properties and
ranges. Select heat treatments may also be performed in order to alter crystalline
structures (Schetky 1981). Perhaps one of the greatest benefits of shape memory alloys is
the relatively high power-to-weight ratio(Ikuta 1990). Or alternately stated, it has a high
work density(Shin, Mohanchandra et al. 2005). For example, an application of pumping
fluids at high pressure may only require a small, light Nitinol shape memory alloy pump
while other traditional pumps are larger and heavier.
14
Nitinol has many useful engineering properties including high fatigue limit,
corrosion resistance(Ford and White 1996)and biomedical compatibility(Duerig, Pelton
et al. 1999). Due to the superior characteristics of Nitinol, it is the only SMA considered
hereafter.
2.2.1 Basic Operation Overview of Shape Memory Alloys
There are two different modes of operation for shape memory alloys, namely,
one-way shape memory effect (OWSME) and two-way shape memory effect (TWSME)
(Adnyana 1984). The two-way shape memory effect involves thermal ?training? such
that no other external forces are needed for the SMA operation. It operates strictly by
simply adding heat and cooling. The two-way SMA mode of operation is utilized in this
new micropump design; however, the one-way is explained due to its simplicity.
The one-way shape memory effect (SME) occurs when one of the select shape
memory alloys (Nitinol) is mechanically deformed while at a temperature lower than the
martensite phase transformation temperature and then heated above the austenite phase
transformation temperature. When the alloy is heated above the critical austenite phase
transformation temperature, the original undeformed specimen is recovered. This process
is best described as a diffusionless transformation in which there is no migration of atoms
within the alloy specimen while transitioning between the two phases, martensite and
austenite. Martensite is the low temperature phase that exists during the mechanical
deformation. Austenite is the high temperature phase that recovers the original
?memory? of the specimen (Miles 1989). The unique trait that distinguishes shape
memory alloys from ordinary metals is their ability to reverse the transformation from
15
austenite to martensite. This can only be accomplished if the material is
crystallographically reversible. Crystallographic reversibility is closely related to ordered
structures. For example, most shape memory alloys have superlattice structures with
atoms in a body-centered cubic (BCC) arrangement (Funakubo 1987).
To graphically illustrate how a shape memory alloy operates, consider the SMA
wire shown in Figure 2.2. The wire is initially below the transition temperature
(martensite). An external force is then applied to ?plastically? elongate the martensitic
wire to some predetermined length (up to a maximum of 8%). Upon the subsequent
removal of the external force, the deformed wire is heated (e.g., resistive heating) until
the wire is above the transition temperature (austenite). This martensite to austenite
transformation produces a contraction force that enables the full recovery of the original
length of the wire. Note that the deforming force for the wire does not have to be
removed in order for the strain recovery. SMA?s can typically recover about three times
the martensite deforming stress (25,000 psi). The fatigue limit of the SMA wire is
directly related to the strain recovered during each deform/recover cycle (Curry 1979).
For a nominal operating strain value slightly greater than 8%, the fatigue limit varies
from ten to one hundred cycles. For a nominal operating strain that is less than 7%, the
fatigue limit can be in excess of one million cycles.
This same cycle is represented at a more crystallographic level in Figure 2.3. The
wire is initially cooled below the martensite transition temperature under zero stress
conditions. Consequently, this results in 24 variants (orientations) of martensite being
produced (Wayman 1980). In the absence of an external force, these 24 variants are very
stable. The two ways to crystallographically change the martensitic specimen are to
16
apply an external force or to heat above the austenite transition temperature. The
application of a force results in only one orientation of martensite prevailing because of
the movement of interfaces (Wayman, 1980). Upon removal of the external force, the
deformed martensite remains elongated in what appears to be a permanent plastic
deformation. However, heating above the austenite transition temperature recovers (or
reverses) the ?plastic? deformation.
Regardless of whether the martensite has one orientation or many orientations,
heating above the austenite transition temperature causes the original austenite to reform.
No motion of the SMA wire occurs if it is heated in the multi-oriented martensite state
(prior to deforming to single orientation). If mechanical work is to be extracted from the
system, the multi-orientated martensite must be subjected to a deforming force prior to
heating.
2.2.2 Shape Memory Alloy Phase Change Temperatures
The characteristic transition temperatures for Nitinol are easily varied by slightly
changing the Nickel/Titanium ratio in the SMA. A 1% shift in the nominally 50-50 blend
of Nitinol can shift transitions by several degrees (20-50 C). The heat-treating process
for Nitinol is discussed here, however, suffice it to say that the transition temperature for
one-way shape memory in Nitinol is almost entirely controlled by metal composition, not
heat-treating(Tuominen and Biermann 1988).
The temperatures at which phase transformation from martensite to austenite
occur are referred to as austenite start (Tas) and austenite finish (Taf). Likewise, the
temperatures at which the reverse phase transformation from austenite to martensite
17
occur are referred to as martensite start (Tms) and martensite finish (Tmf) (Pelton, Duerig
et al. 2005). A complete thermal cycle is depicted in Figure 2.4 where the martensite
volume fraction (?) is plotted against temperature (Thrasher, Shahin et al. 1992). The
cooling portion of the thermal cycle is considered the forward transformation while the
heating is considered the reverse transformation (Stoeckel 1989).
As evident in Figure 2.4, the phase transformation does not occur instantaneously,
but rather it occurs gradually over a 10-20 oC temperature range (Tas-Taf or Tms-
Tmf)(Shahin, Meckl et al. 1994). Another notable feature evident in Figure 2.4 is the
hysteresis in the SMA thermal cycle. A typical hysteresis value for Nitinol is
approximately 20 oC. Unlike the transition temperature that is easily shifted when
changing the stoichiometric (50-50) Nickel-Titanium alloy, the hysteresis is constant and
not strongly dependent on the exact ratio of Nickel and Titanium. The hysteresis is
specifically defined as the temperature difference between Taf and Tms or Tas and Tmf.
The thermal cycle shown in Figure 2.4 is very repeatable over an unlimited
number of thermal cycles if it occurs under a zero-stress condition. However, the
introduction of an external force (stress) results in a small ?almost linear? shift in the
transition temperatures (Tas, Taf, Tms, Tmf). The relationship between martensite and
austenite transformation temperatures is shown in Figure 2.5 as a function of externally
applied stress(Golestaneh 1984). Each of the four lines represents the shift of the
corresponding transition temperature, and the slope ?C? of each line is approximately
equal. All of the subscripts ?o? represent the zero stress condition.
The transformation lines for the transformation from martensite to austenite have
some nonlinearity associated with them at low stress conditions (Tanaka 1986). At these
18
lower stresses, the transformation lines for Tas and Taf drop off more rapidly than the
linear martensite transformation lines (Tms and Tmf). This nonlinearity is not well
understood and consequently is often neglected. For Nitinol, the slope, ?C?, of the stress-
induced transition lines is 1000 psi/oC (Tanaka 1990).
2.2.3 SMA Transformation Kinetics
Thermodynamics is very useful in determining the driving force for a phase
transformation, but it doesn?t reveal anything regarding the rapidity of the transformation.
The Arrhenius rate equation is derived from kinetic theory to estimate the probability of
an atom overcoming the activation free energy barrier (Porter and Easterling 1981). The
exponential function gives good results for the martensite fraction (?) around the
beginning of the austenite to martensite transformation (Liang 1990), but the results at the
end of the transition are not as good. The exact same thing is said for the reverse
transformation. This shortcoming is circumvented by using a simpler model that utilizes
cosine functions instead of exponential functions(Liang and Rogers 1990). This cosine
model easily compensates for the asymptotic nature at both ends of the transformation
region (Vitiello, Giorleo et al. 2005) (Figure 2.4).
The martensite fraction for the martensite to austenite transformation is
{ }1)](cos[21 +?= asTTAa?
(2.1)
And the martensite fraction for the reverse, austenite to martensite, transformation is
19
{ }1)](cos[21 +?= mfTTMa?
(2.2)
where,
)(
asTafTA
a ?= pi (2.3)
and )(
mfTmsTM
a ?= pi (2.4)
These equations apply if the SMA specimen is heated and cooled while under zero stress
conditions. Using the slope, ?C?, of the stress-shifted transition temperatures in Figure
2.5, stress-compensating forms of the previous equations yield
2])(cos[2
MAbasTTAaM ???? ++?=
(2.5)
for the martensite to austenite transformation and
2
1])(cos[
2
1 A
MbmfTTMaA
???? +++??=
(2.6)
20
for the austenite to martensite transformation, where bA =-aA/C and bM=-aM/C. ?A and ?M
are the martensite fractions at the beginning of the martensite to austenite transformation
and the austenite to martensite transformation, respectively. These martensite fractions
are always ?0? or ?1? if the transformations are complete for each leg of the thermal
cycle.
2.2.4 SMA Mechanical Properties
2.2.4.1 Stress-Strain
Some of the most fundamental and useful mechanical properties of interest are
derived from the stress-strain plot. A stress-strain plot provides a quick visual insight
such that the modulus of elasticity, yield stress and ultimate strength are easily
ascertained. Figure 2.6 shows the stress-strain curve for a shape memory alloy (Nitinol).
Obviously, the stress-strain curve for Nitinol is not typical and requires further insight for
explanation.
Four distinct regions of the Nitinol stress-strain curve are indicated with letters in
Figure 2.6. The strain-strain data is obtained using an Instron stress machine.
Additionally, the assumption that the Nitinol specimen is cooled below the martensite
transformation temperature is enforced.
In the (a-b) region of the graph, multi-oriented martensite is elastically deformed
up to the first yield point (b). At point (b), the multi-oriented martensite begins
detwinning and is reoriented until there is a single orientation of martensite (c). For
typical materials, this region (b-c) is plastic deformation that is unrecoverable. However,
21
this is the region that makes shape memory alloys unique to other materials because this
strain (around 8%) is fully recoverable once it is heated above the austenite
transformation temperature. Region (c-d) is elastic loading of the deformed, single
oriented martensite. At point (d), permanent unrecoverable slippage occurs. This
slippage continues until the ultimate fracture point is reached at point (e). The region (c-
d-e) is the part of the stress-strain curve that is analogous for other materials.
Data for the (a-b-c) region of Nitinol SMA in martensitic phase is shown in Figure 2.7
(Thrasher 1991). There is an obvious difference between the performance of the
?trained? and ?untrained? SMA specimens. An SMA specimen that is only annealed and
not cycled is considered as being ?untrained?. The training process of an SMA is
accomplished by elongating the wire and then heating it while under load until full
recovery is complete. This process is repeated for several cycles until the SMA
performance is reproducible for each cycle. The graph illustrates the importance of
?training? an SMA before designing and implementing the SMA in an actuator. Table
2.4 contains the mechanical data for both trained and untrained Nitinol SMA specimens.
Martensitic SMA
(Nitinol)
Modulus of Elasticity
(psi)
Yield Strength
(psi)
Untrained
20e6
17,500
Trained
20e6
23,500
Table 2.4 Mechanical Data for Trained and Untrained Nitinol SMA
22
The modulus of elasticity of Nitinol as a function of temperature is shown in
Figure 2.8. The modulus of elasticity changes by a factor of approximately three between
the martensite and austenite crystal structure. The dynamic stiffness of the SMA
materials is very useful in structural vibration applications (Rogers, Liang et al. 1991).
2.2.4.2 Recovery Force
A recovery force (or stress) is generated when SMA specimens are being heated
to convert the crystalline structure from martensite to austenite. A fundamental
knowledge of the generated force is important since it provides the upper limit of net
mechanical work that can be extracted from an SMA actuator in one cycle of operation.
A simple scenario to depict a typical cycle of SMA operation is shown in Figure 2.9. The
recovery stress recovered is represented by W2 while W1 is the required weight that is
just heavy enough to deform the SMA.
At state 1 in Figure 2.9, a shape memory alloy (SMA) wire of length ?L? is
attached to a fixed structure. A weight (W1) just large enough to deform the wire is
applied at state 2 to deform the SMA by an amount ??L?. It is assumed that the elastic
portion of the deformation is reflected in the ?L but it is never recovered as long as W1 is
applied. Furthermore, a stop will be placed at the bottom of the stroke in order to ensure
that the deformation ?L does not exceed the maximum recovery strain limit (6-8%).
After W1 reaches the stop, an additional weight, W2 , is applied. This represents the
recovery force/stress that the SMA generates during the return stroke. The SMA wire is
heated during the process 3-4 to produce a mechanical work output. After complete
contraction of the SMA wire (minus the elastic strain of austenite), weight W2 is removed
23
(step 5). Subsequently, the SMA is cooled until the stress generated by W1 is sufficiently
large enough to deform the wire back to the stop (step 2). This cyclical process may be
repeated until the endurance limit of the SMA wire is reached. As mentioned earlier, if
the strain is not excessive (<7 %), it should have an infinite endurance limit.
To facilitate understanding, another representation of this cyclical process is
shown graphically in Figure 2.10 in the stress-strain diagram. Step 1 to 2 is merely the
stress-strain curve for martensitic SMA. The amount of strain (up to 8% maximum) is
controlled by placing a stop under the weight. At step 3, an additional weight is added,
and this weight does not increase the stress in the SMA because the stop absorbs the extra
loading. After applying the extra weight (W2), the wire is heated and results in an
isostrain region while going from state 3 to 4. Motion of the two weights does not occur
until the recovery stress in the SMA exceeds the stress caused by the combined weights.
Once contraction commences, it will continue from state 4a to 5 under constant stress
condition. The contraction continues until the only strain left is that of the elastic
austenite SMA. The limiting strain energy is for austenite rather than martensite since
the SMA will have austenite crystalline structure at this high temperature. Removing the
extra weight at step 4 causes the strain to progress along the elastic line of the austenite.
Once at step 5, the SMA is cooled and the cycle repeats.
2.2.5 SMA Physical Properties
Table 2.5 lists some of the most commonly used physical properties for Nitinol
SMA (Ikuta 1990). The values are considered applicable for all blends of Nitinol close to
50-50 (USNitinol 1989). Both SI and English units are listed.
24
2.2.6 SMA Chemical Properties
For many applications, chemical properties such as resistance to corrosion or
oxidation may be of significant importance. There are some applications in which the
SMA material is used under water. These issues become more relevant under these
circumstances.
Resistance to corrosion has been tested by subjecting a Nitinol sample to a 96-hour spray
of salt water at 35 oC (Jackson, Wagner et al. 1972). No corrosion is evident. Also, no
corrosion is evident when a Nitinol sample is exposed to 60 days of high-velocity sea
water spray at 26 ft/sec (Jackson, Wagner et al. 1972). This demonstrates the capability
of using Nitinol submerged in water.
Physical Properties for Nitinol SI
(English)
Specific Heat 250 J/kg-K
(0.0597 BTU/lbm-F)
Latent Heat of Transformation 30,000 J/kg
(12.899 BTU/lbm)
Density 6.45 g/cm3
(402.66 lbm/ft3)
Melting Point 1200-1300 oC
(2192-2372 oF)
Electrical Resistivity 76 microhm-cm (Mart.)
(193 microhm-in)
82 microhm-cm (Aus.)
(208 microhm-in)
Thermal Conductivity 0.17 watt/cm-C
(9.823 BTU/h-ft-F)
25
Coefficient of Thermal Expansion 6.6 x 10-6/C (Mart.)
(3.7 x 10-6/F)
11.0 x 10-6/C (Aus.)
(6.11 x 10-6/F)
Table 2.5 Physical Properties of Nitinol SMA
The rate of oxidation for Nitinol increases very rapidly when the temperature is
increased above 600 oC. A sample subjected to this temperature experiences a weight
gain of only 0.002 grams when exposed for 10 hours. At 1000 oC, the sample has a total
weight gain of 0.1 grams for the same 10-hour exposure time.
2.2.7 SMA Biocompatibility
Shape memory alloys, most notably Nitinol, were first used in medical
applications in the early 1970?s. The first application was an intrauterine contraceptive
device (Hodgson 2005). Some common applications developed since then are orthopedic
implants, bone staples, spinal spacers (Petrini, Migliavacca et al. 2005), stents, inner ear
implants, tendon suture material (Kujala, Pajala et al. 2004) and artificial anal sphincters
(Yun and et al. 2003). Another growing field in medical applications is the regeneration
mechanism for tissues and organs. Porous biocompatible implants are created by laser
sintering of Nitinol. This was initially done with pure titanium (Shishkovsky 2008).
With a growing demand for Nitinol applications in the human body, the biocompatibility
is of utmost importance and must be evaluated.
As previously mentioned, the most prevalent and popular shape memory alloy is
Nitinol. The two components of Nitinol are nickel (Ni) and titanium (Ti). Titanium is
26
extremely corrosion resistant and has no known biocompatibility problems (Gil and
Planell 1998). Pure titanium and Nitinol have no adverse effects on the human body and
do not alter the physiological or behavioral functions of cells. Conversely, nickel in large
quantities is detrimental to cells and tissue. Nickel can exist in the human body in
elemental form in trace quantities only. These trace quantities are essential and can enter
the human body via lungs, skin, or orally. From food alone, the human body ingests
between 80 and 400 ?g/day/person (Es-Souni, Es-Souni et al. 2005).
Unlike pure titanium and Nitinol, nickel in pure form is not very corrosion
resistant or biocompatible. Additionally, as mentioned above, nickel (g1840g1861g2871) ions in
appreciable quantities may be very harmful to the human body. Some of the major
concerns for pure nickel are cytotoxicity, carcinogenicity, and tissue necrosis.
Electrochemical and cell culture tests reveal that Ni in its pure form is not biocompatible.
The most common method of controlling the release of Ni from Nitinol is surface
oxidation (Firstov, Vitchev et al. 2002). Many oxidation methods may be utilized,
however, simple oxidation with air is sufficient. The oxidation temperature varies in the
range of 300-800 C for various lengths of time. The oxidation layer (TiO2) may be
scanned with scanning electron microscopy (SEM). The titanium is a more reacting
element than nickel, and it segregates to the surface. This forces the Ni enriched layer to
form beneath the outer oxidation layer. Oxidation at temperatures above 500 C or higher
produces a Ni-free outer layer and renders the Nitinol as biocompatible (Wang 1998).
With proper oxidation of the surface, the Nitinol implant releases less Ni than the best
stainless steels (316L) (Ryh?nen, Kallioinen et al. 1999; Mantovani 2000; Hodgson 2005;
Michiardi, Aparicio et al. 2007). An alternative method for surface oxidation is
27
electropolishing in methanol-sulphuric acid electrolyte solution (Barison, Cattarin et al.
2004). A solution of g1834g2870g1841g2870 may be used to improve the bond between the g1846g1861g1841g2870and the
underlying Nitinol matrix (Hu, Chu et al. 2007).
A separate focus on biocompatibility is the bioactivity between the Nitinol and
human tissue. In orthopedic implants, it is important for bone tissue and Nitinol to bond.
This bonding may be accelerated if a layer of hydroxyapatite (HA) is applied to the
Nitinol bone implant. This HA layer is created by treating the implant with HNO3 and
NaOH aqueous solution (Chen, Yang et al. 2004). Bioactivity may also be improved by
using a technique of plasma immersion ion implementation and deposition (PIII&D)
(Chu 2006) (Yeung, Chan et al. 2007).
The last concern for biocompatibility is the effect of the elevated SMA
temperature on the human body. Research has shown that as blood gets warmer,
coagulation (clotting) is reduced. Once the temperature of the blood is in excess of 55 C,
there is no coagulation (Simpson and Rasmussen). The last issue is the effect of high
blood temperature on surrounding cells and tissue in the body. The body tissue can
handle temperatures as high as 42 C before damage occurs (Yun and et al. 2003). It is
shown in the thermal section that the average blood temperature exiting the micropump is
below 42 C.
2.3 Lattice Boltzmann Method
28
2.3.1 General Overview
The study of fluid dynamics is a rather complex problem because the flow may
vary from simple laminar flow in a pipe to complex vortex shedding over an airfoil.
Using conservation equations on small differential volumes, the continuity and
momentum equations (Navier-Stokes) can be derived. However, there are no analytic
solutions based upon these equations for a general fluid flow. Traditionally, the fluid
flow is analyzed numerically by discretizing the nonlinear partial differential equations
using finite elements, finite differences, etc.
An alternate approach to solving the Navier-Stokes equations is using particle-
based methods (Reider and Sterling 1995). Some of the more popular particle methods
are as follows:
? Direct Simulation Monte Carlo (DSMC)
? Molecular Dynamics (MD)
? Lattice Gas Method (Bernsdorf 2001)
There are many advantages and disadvantages for each of these methods. Suffice it to
say that the focus of this research is the implementation of Lattice Boltzmann method
(LBM).
The immediate predecessor to the Lattice Boltzmann method is the Lattice Gas
Automata. Unlike the Lattice Gas method which looks at the Boolean state before and
after the particle collisions (Luo 2003; Kang 2005), the Lattice Boltzmann Method uses a
distribution function that is a real number (Love and Boghosian 2004). This eliminates
the noise problem due to binary encounters with the Lattice Gas (Qian, D'Humieres et al.
1992). Some of the primary advantages of LBM are
29
? Highly Parallelizable (Satofuka and Nishioka 1999; Desplat,
Pagonabarraga et al. 2001; Pohl, Deserno et al. 2004)
? Easily adapted to complex geometries (Martys 2001; Li, LeBoeuf et al.
2004)
? Equally suited for rarefied and/or continuum flows
? Great for high Reynolds number flows
? Easy computer implementation
? Multicomponent Fluids (Nie, Qian et al. 1998)
Likewise, several of the most significant disadvantages are
? Velocity of flow limited to < Mach 0.3
? ?Weakly? compressible (Sun and Hsu 2003; Kataoka and Tsutahara 2004)
? Very small time steps due to explicit nature of equations
? Conditional stability (unless using Entropic LBM) (Sterling and Chen
1996; Lallemand and Luo 2000; Hinton, Rosenbluth et al. 2001; Sun and
Hsu 2004)
For the study of the SMA micropump, most of these disadvantages are not relevant.
Some general applications for LBM are amphiphilic fluids (Chen, Boghosian et
al. 2003), Womersley flow (Rohlf and Tenti 2001; Artoli, Hoekstra et al. 2006), Kelvin-
Helmholtz instabilities (Zhang, He et al. 2001), interface and surface tension (Zhang, He
et al. 2000), pulsatory channel flows (Majdalani and Chibli 2002), two-phase flows
(Swift, Osborn et al. 1995; Do-Quang, Aurell et al. 2000; Watanabe and Ebihara 2002),
binary fluid systems (Swift, Orlandini et al. 1996; Dupin, Spencer et al. 2004), obstructed
flows (Takada and Tsutahara 1998; Shi, Lin et al. 2003; Li, Shock et al. 2004),
30
cylindrical Couette (Halliday, Hammond et al. 2001; Aoki, Yoshida et al. 2003),
electrorheology (Melchionna and Succi 2004), static mixer reactor (Kandhai, Vidal et al.
1999), blood flow (Fang, Wang et al. 2002; Hirabayashi, Ohta et al. 2004), bifurcation
(Sone and Doi 2000), electron transport (Date and Shimozuma 2001), microchannels
(Arkilic, Schmidt et al. 1997; Rector, Cuta et al. 1998; Rasmussen and Zaghloul 1999;
Quadir, Mydin et al. 2001; Lim, Shu et al. 2002; Nie, Doolen et al. 2002; Sinton,
Erickson et al. 2002; Erickson and Li 2003; Li and Kwok 2003; Li and Kwok 2003; Sert
and Beskok 2003; Li and Kwok 2004; Raabe 2004; Lee and Lin 2005; Verberg, Yeomans
et al. 2005; Xuan and Yao 2005; Zhang, Qin et al. 2005), porous structures (Inamuro,
Yoshino et al. 1999; Guo and Zhao 2002; Yoshino and Inamuro 2003; Pan, Prins et al.
2004; Tang, Tao et al. 2005), bubble flow (Sankaranarayanan, Shan et al. 1999; Holdych,
Georgiadis et al. 2001), Elder problem (Thorne and Sukop 2004), adsorption
breakthrough (Agarwal, Verma et al. 2005), buoyancy-driven flows (De La Fuente,
Causon et al. 2003; Zhou, Zhang et al. 2004), turbulence (Luo, Qi et al. 2002),
viscoelastic media (Qian and Deng 1997) and cavity flow (Miller 1994; Hou, Zou et al.
1995; Onishi, Chen et al. 2001; Wu and Shao 2004).
2.3.2 Boltzmann Equation
The Boltzmann equation provides a statistical description of a fluid by using a
density distribution function (Yu, Mei et al. 2003). This distribution function, f, has a
dependence upon space, velocity and time (Nekovee, Coveney et al. 2000). Neglecting
external forces acting upon a volume of fluid, the Boltzmann equation (Doolen 1987)
implementing the distribution function and written in differential form is
31
?=?+?? ftf (f) (2.7)
where
f = density distribution function
? = collision function
Stated in simpler terms, the distribution function, f(r,?,t)drd?, reflects how many
molecules at a time, t, positioned between r and r + dr which have velocities in the range
of ? and ? + d?.
The collision function from kinetic theory is a rather complicated function. For
local equilibrium, the gains and losses due to collisions must be equal. There needs to be
an exact balance such that the collision term drives toward zero.
2.3.3 Lattice Boltzmann BGK
The most popular numerical implementation of the Boltzmann equation is the
BGK model. This model uses a set of lattice velocities and a simplified collision function
to drive towards local equilibrium. The relaxation time of the collision function is
numerically linked to the viscosity in order for Navier-Stokes to be recovered.
32
By using a single relaxation time, ?, for the collision term (He and Luo 1997; Martys,
Shan et al. 1998; Fang, Wan et al. 2002) in the Boltzmann equation, the Boltzmann
equation can be approximated as
)(1 eqffftf ?=??+?? ?? (2.8)
where the Maxwellian equilibrium distribution is given by
??
?
?
???
? ??=
RT
u
RTf
eq
2
)(exp
2
2?
pi
? (2.9)
? is the molecular velocity while u is the macroscopic velocity. Since the relaxation is a
simple parameter (and everything relaxes at the same rate), the conservation of the
hydrodynamics must be ensured by the ?g1858g3032g3044? part of the collision operator. The link
between the relaxation time and viscosity is based upon the continuum assumption.
These are discussed in greater detail next.
The macroscopic variables, ? and u, are defined as (Abe 1997; Buick and Greated 2000;
Shu, Niu et al. 2005)
?= ?? fd (2.10)
?= ??? fdu (2.11)
33
From performing a formal integration of Equation 2.8 (He and Luo 1997; Kandhai,
Koponen et al. 2001) over time, ?t,
)],,(),,([1),,(),,( txftxftxftxf eqtt ???????? ??=?++ (2.12)
where ? = ?/?t. The feq term is calculated using a Taylor series expansion. After a Taylor
series and keeping terms up to order O(u2),
??
?
??
? ??+?+
???
?
???
? ?=
RT
u
RT
u
RT
u
RTRTf
eq
2)(2
)()(1
2exp2
2
2
22 ???
pi
? (2.13)
To recover the Navier-Stokes equations, a Chapman-Enskog multiscale expansion is
used. The details of that derivation are omitted.
For the purpose of this research, the LBM element of choice is the 2-D element
that has 9 velocities. These are shown pictorially as
0
8 7
6 5
4
3
2
1
34
where the 9 velocities are listed as
?
?
?
?
??
?
?
?
?
?
?
?
??
?
?
?
+?==
?==
=
=
4)5(2,8,7,6,5,)sin,(cos2
)1(2,4,3,2,1,)sin,(cos
0),0,0(
pi?pi????
?pi????
?
?
c
ce (2.14)
Speed ?c? is the lattice speed (He and Luo 1997), ?x/?t. With much algebraic
manipulation and Chapman-Enskog expansion (Qian and Zhou 2000; Guo, Zheng et al.
2002), the distribution equilibrium function is listed as
??
?
??
? ??+?+=
2
2
4
2
2 2
3
2
)(9)(31
c
u
c
ue
c
uewf eq ??
?? ? (2.15)
where
??
???
??
???
=
=
=
=
8,7,6,5,36/1
4,3,2,1,9/1(
0,9/4
?
?
?
?w (2.16)
2.3.3.1 Algorithm for Lattice Boltzmann Method Solver
The numerical steps for a solution are fairly straightforward. The procedure is
1. Initialization of all parameters
2. Propagate all distributions to adjacent nodes
35
3. Calculate new velocities
4. Calculate new equilibrium
5. Evaluate collision function
6. Return to step 2
2.3.3.2 Boundaries for Lattice Boltzmann Method
There are three primary boundary types for fluid analysis: solid, velocity, and
pressure. Like any other numerical scheme, the boundaries are implemented in LBM as
either first order or second order. Some good detailed descriptions of second order
boundaries are available in literature (Chen, Martinez et al. 1996; Ginzburg and
D.d'Humieres 1996; Maier, Bernard et al. 1996; Zou and He 1997; Guo, Zheng et al.
2002; Verberg and Ladd 2002; Ginzburg and D'Humieres 2003; Rohde, Kandhai et al.
2003; Tang, Tao et al. 2005). Of the boundary treatments, some account for moving
boundaries (He and Q. 1995; Noble, Chen et al. 1995; Fang, Lin et al. 1998; Bouzidi,
Firdaouss et al. 2001; Li, Lu et al. 2004), curved boundaries (Mei, Luo et al. 2000), non-
slip (Inamuro, Yoshino et al. 1995), and velocity/pressure (Lin, Fang et al. 1996; Guo,
Zheng et al. 2002; Lim, Shu et al. 2002). This is discussed in greater detail when analysis
of the SMA moving diaphragm is considered.
2.3.4 Stability of Lattice Boltzmann Method
One of the biggest disadvantages of the LBM is the numerical stability. The
numerical instability poses problems for some simulations such as very high Reynolds
numbers or compressible problems (Ansumali and Karlin 2000; Ansumali 2004). The
36
instabilities result from the distribution function going negative in regions of phase space.
The method needs an analog to the Boltzmann H-theorem in order to ensure the entropy
always increases (Boghosian, Love et al. 2003).
One method for eliminating the numerical stability is the entropic lattice
Boltzmann method (ELBM) (Karlin, Ansumali et al. ; Ansumali and Karlin 2002;
Ansumali and Karlin 2002; Ansumali, Karlin et al. 2003; Boghosian, Love et al. 2004;
Ansumali and Karlin 2005; Gorban and Karlin 2005). The ELBM begins by positing an
H-function depending on the single particle distribution function (Boghosian, Yepez et al.
2001). The relaxation time is then dynamically adjusted such that H-theorem is
guaranteed (Boghosian and Boon 2005). The smallest relaxation time that doesn?t
increase the entropy is the one selected (Boghosian, Love et al. 2004).
The ELBM method maintains all of the advantages of the LBM method, including
the solution of microflows (Ansumali, Karlin et al. 2006). ELBM is implemented for
both 2-D and 3-D flows (Boghosian 2003). The most recent research for the ELBM
method is that for thermal flows (Chen and Teixeira 2000; Chikatamarla, Prasianakis et
al. 2005). The instabilities mentioned above have not adversely affected any of the
ongoing research of micropumps. Therefore, the LBM is the only method considered in
this dissertation.
Figure 2.1 Breakdown
Displacement
Reciprocating
Diaphragm Piston
37
of the Different Types of Pumps
Micropumps
Rotary
Dynamic
Centrifugal Electrohydrodynamic
Electroosmotic
38
Figure 2.2 Illustration of SMA Wire Recovering Original Shape
39
Figure 2.3 Crystallographic view of Shape Memory Effect
40
Figure 2.4 SMA Composition as a Function of Temperature
41
Figure 2.5 Dependence of Transition Temperatures upon Applied Stress
42
Figure 2.6 Shape Memory Alloy Stress-Strain Curve to Ultimate Strength
43
Figure 2.7 Stress-Strain Curves for Trained and Untrained Nitinol SMA
44
Figure 2.8 Modulus of Elasticity as a Function of Temperature
45
Figure 2.9 Work Extraction from SMA Actuator
46
Figure 2.10 Detailed Breakdown of SMA Actuator Cycle
47
Chapter 3 Shape Memory Alloy Diaphragm Model
3.1 Micropump Overview
Figure 1.1 in Chapter 1 shows a basic schematic of the proposed SMA
micropump design. Sputtering Nitinol onto a silicon substrate creates the SMA
diaphragm shown in the micropump. The shape memory alloy (SMA) diaphragm is
resistively heated in order to cause contraction (martensite to austenite crystalline
conversion). During the contraction, the check valve (CV) on the left blocks flow while
the CV on the right allows outward flow of the medical fluid. The volume of the flow
with each cycle is prescribed by the shape of the SMA membrane.
After the SMA is heated and fully contracted, the resistive heating is cycled off
such that cooling begins. During the cooling, the check valves behave the opposite as
mentioned above. A continuous oscillatory flow of blood is used to enhance the cooling
of the SMA membrane. As the SMA membrane is cooled and returned to its initial state,
the CV on the left allows the medical fluid to be pulled into the upper chamber while the
CV on the right blocks flow. This process is repeated as needed for the required flow rate
of the medical fluid.
The functionality of the micropump depends almost entirely upon the dynamics of
the SMA diaphragm and the check valves location. The SMA diaphragm response is
directly linked to the material phase transformation, while the rapidity of the phase
transformation is controlled by either the resistive heating power density (Watts/kg) or
48
the cooling from fluid contact. This chapter focuses upon the heating; hence, the
convection/conduction with the fluid is discussed later. In order to model the cooling
phase of the SMA diaphragm, knowledge of fluid temperature and velocity is essential.
This is discussed in a later chapter.
3.2 SMA Micropump Diaphragm
The SMA diaphragm used in this micropump is ?trained? as a two-way memory.
This means that no bias forces are required and the entire motion of the SMA diaphragm
is solely controlled by the input or extraction of heat for the phase transformation
kinetics. The rest of this chapter focuses on these transformations.
3.2.1 Phase Transformation Model for SMA Diaphragm
The kinetics considered in this discussion are for the heating of the SMA
diaphragm. The complete cycle for the SMA entails being cooled and then heated in
order for the strain in the deformed SMA to be recovered. The method of heating is
resistive heating. Cooling is accomplished via conduction/convection from the
diaphragm into the pumped fluid. Details of the cooling portion of the cycle are
discussed in Chapter 6.
Figure 2.3 shows the Martensite fraction as a function of temperature. Due to the
asymptotic nature of the curve, the best model for martensite fraction, ?, as a function of
temperature is the Cosine Model (CM). CM is described with the equations 2.1 through
2.6.
49
3.2.2 Thermal Model for SMA Diaphragm
The thermal model of the micropump depends upon resistive heating and heat loss
via convection/conduction to a cooling fluid. With the assumption that the resistive
heating is much quicker than the cooling to the fluid (blood), the heating part of the cycle
is easily computed.
The internal energy generation in the SMA diaphragm due to resistive heating
(watts) is
RIE 2= (3.7)
where ?I? is the current in amps and ?R? is the electrical resistance of the Nitinol SMA.
For a specified value of ?E?, the temperature, hence the strain, of the SMA diaphragm is
computed. From the basic internal energy equation (while temporarily neglecting the
effect of cooling), the temperature change per unit of time of heating is
pmC
E
dt
dT =
(3.8)
Once the temperature reaches Tas, the specific heat is modified to account for the latent
heat of transformation, ?h?. The latent heat is assumed to vary linearly with the
martensite fraction over the temperature range from Tas to Taf.
asaf
pp TT
hCC
?+= (3.9)
50
This modified specific heat effectively slows the heating/cooling of the SMA during the
phase transformation. This is illustrated graphically in Figure 3.1 for a power input of
0.02 watts. (This wattage is typical of what is modeled later in order to achieve a heating
time about 0.05 sec during each thermal cycle). The total mass of the SMA in this
scenario is 0.00258 kg. Note that this curve temporarily neglects conduction/convection
to the cooling fluid. It is only to illustrate the slope change of the heating curve during
the phase transformation. The slope of the heating curve is drastically reduced while
going through the phase transformation from martensite to austenite.
3.2.3 SMA Diaphragm Deformation Geometry
Figure 3.2 provides a labeled geometric schematic of the SMA diaphragm while
being heated. ?S? is the actual length of the deformed SMA diaphragm while ?C? is the
length of the undeformed SMA diaphragm that is recovered while heating the SMA
above the transition temperature. The strain of the SMA is defined as ???. Hence, the
length of ?S? is
CS )1( ?+= (3.10)
where ? << 1. The value of ?C? depends upon the two-way training of the SMA and is
fixed. It is also implicitly assumed that the SMA diaphragm preserves the chord of a
perfect circle. The strain for Nitinol is usually limited to 8%, or less, if many cycles of
operation are needed. Other useful relationships are
51
?RS = (3.11)
and
??????= 2sin2 ?RC
(3.12)
Combining the previous three equations yields the relation
??????+= 2sin)1(2 ???
(3.13)
Equation 3.13 is iteratively solved for ??? and then ?R? is calculated from equation 3.11.
The distance?d? is calculated as
??????= 2cos ?Rd
(3.14)
And finally, the height, ?h?, at which the diaphragm penetrates the chamber is prescribed
by equation 3.15 as
dRh ?= (3.15)
52
The equations listed above apply to the SMA diaphragm when it is being heated
or cooled. By knowing the energy density input into the SMA, the strain ??? is easily
calculated since the strain is directly related to the martensite fraction during the heating
or cooling cycle. By knowing the kinetics of the phase transformation within the SMA,
the details of the oscillating wall movement in the cavity are known.
3.2.4 SMA Diaphragm Kinematics
Figure 3.3 conceptually shows the position of an SMA diaphragm at two different
consecutive time steps. The SMA at position 1 is initially a deformed and cooled SMA
diaphragm. As electrical current is conducted through the diaphragm, the internal heat
generation in the SMA drives the material phase transformation towards austenite, and
the SMA contracts to position 2. Using the assumption that the SMA diaphragm
preserves the chord of a perfect circle during contraction, the geometry and velocity of
the diaphragm are calculated as a function of SMA strain (which is predictable by using
phase transformation kinetics). The SMA strain, ???, depends only upon temperature;
therefore, it is controlled by the rate of resistive heating/fluid cooling. Hence, the strain
is directly connected to time and rate of heating/cooling. After each time step, the new
strain of the contracting SMA diaphragm is calculated using transformation kinetics (Sec.
3.2.1). From knowledge of the strain, the above equations for the angle, ???, and the
radius, ?R?, predict the geometry of the contracting diaphragm. Note that during the
contraction, the angle, ???, decreases while the radius, ?R?, increases. As the diaphragm
flattens near the end of the contraction cycle, the angle approaches zero while the radius
approaches infinity.
53
The next aspect of interest in the contracting diaphragm is the velocity of the
membrane. Before solving for velocity, the complete geometry of the contraction must
be known. As shown in Figure 3.3, the radius, R, and the angle, ?, are shown for two
consecutive positions. The angle, ?, is solved by implicitly using equation 3.13.
Subsequently, R is solved using equation 3.11. The next unknown to be solved for is the
distance between the circle centers, ?y. This is shown pictorially in Figure 3.4. The only
way to solve for this parameter is to look at the position in which the end points for g1844g2869
and g1844g2870 coincide. As shown in Figure 3.4, the law of cosines are utilized to solve for the
only unknown in the triangle, ?y. From implementing the law of cosines,
)cos(2 2122212 aRRRRy ?+=? (3.16)
where,
2
2?=b (3.17)
2180
1??=c (3.18)
cba ??= 180 (3.19)
Figure 3.5 graphically shows the displacements, ?x and ?y, for an arbitrary location
between two time steps. The angles, g2038g2869 and g2038g2870, represent the angle between circle center
54
point and the ?x? location of interest. For small time steps, which are typical for the
Lattice Boltzmann Method, the angles g2038g2869 and g2038g2870 are assumed as approximately equal.
This is a needed assumption when solving for ?x. To solve for velocity of the diaphragm,
?x and ?y is calculated at every ?x? location. Assuming that g1876g2869 is a location of interest for
the diaphragm at an arbitrary time step, then
)cos( 111 ?Rx = (3.20)
)cos( 221 ?Rx = (3.21)
Using trigonometric analysis,
yRRy ?+?= )sin()sin( 2211 ??? (3.22)
With knowledge of ?y and the average ? for a particular location, the final calculations
for the diaphragm motion are
?????? += 2 21 ??? (3.23)
and
55
)tan(?
?? yx = (3.24)
The velocities along the diaphragm are now simply stated as
t
xV
x ?
?= (3.25)
t
yV
y ?
?= (3.26)
3.3 Numerical Analysis of SMA Diaphragm
Prior to performing a thorough analysis of the fluid in the micropump chamber, it
is essential to know the velocity of the SMA diaphragm (moving boundary problem).
Using all the information developed in earlier sections in this chapter, the diaphragm
motion is well understood. The next two sections address the velocity and displacement
of the transient motion of the SMA diaphragm.
3.3.1 Transient Displacement Analysis
The displacement of the SMA diaphragm during contraction/expansion is directly
controlled by the heat transfer rate, which is accomplished via resistive heating and fluid
cooling. As the SMA diaphragm is heated, the phase transformation from martensite to
austenite progresses. The exact opposite results when cooling. Using the geometric
analysis above the displacement of the SMA diaphragm (at the center point of the
56
diaphragm) in a micropump is shown in Figure 3.6. The fully contracted length of the
SMA diaphragm, when displacement is zero, is approximately 2 cm for this particular
pump design.
3.3.2 Transient Velocity Analysis
As the heated diaphragm contracts, the X and Y velocities are calculated by using
the equations above to calculate the deflection of each point for each time step. Stating
the obvious, the velocities at the end points of the diaphragm are zero. Additionally, the
X velocity at the center of the diaphragm is zero due to symmetry. This is shown
graphically in Figure 3.7, which shows the X velocity along the SMA for a typical time
step. The maximum X velocity at which the diaphragm moves as a function of time is
shown in Figure 3.8. This time of maximum velocity is controlled by the rate of heating.
Similar to the X velocity, the maximum Y velocity for complete heating transient is
shown in Figure 3.9. Figure 3.10 graphically shows the Y Velocity for two different
locations along the SMA; the center point and the quarter point. All of these velocities
(both X and Y) behave as expected.
57
Figure 3.1 SMA Thermal Response throughout the Phase Transformation (Typical results
with an SMA mass of 0.00258 kg)
0
20
40
60
80
100
120
0 0.02 0.04 0.06 0.08 0.1 0.12
Te
mp
er
atu
re
(C
)
Time (s)
Heating Response of SMA throughout Phase Transformation
(Heating Power Input = .02 watts)
Begin Phase
Transformation
End Phase
Transformation
58
Figure 3.2 Basic Geometry Variables for SMA Diaphragm
Cooled, Deformed
Heated, Undeformed
s
R d
?
h Heat
C
59
Figure 3.3 Change in SMA Geometry after One Time Step (parameters are not to scale)
SMA (pos 1)
SMA (pos 2)
?
?
1
2
R1 R1
R2 R2
60
Figure 3.4 Further SMA Geometric Analysis
SMA (pos 1)
SMA (pos 2)
?2
R1
R2
?y
a
b
c
a=180-b-c
b= 2
?1c=
2180-
61
Figure 3.5 SMA Diaphragm Displacement Analysis
?
?
R1R2
1
2
?
?
x
y
SMA (pos 1)
SMA (pos 2)
? y
62
Figure 3.6 Motion of Diaphragm from Maximum Deflection Point at the Center of
Diaphragm. (The larger the negative number, the greater the penetration into the cavity.)
-0.0016
-0.0014
-0.0012
-0.001
-0.0008
-0.0006
-0.0004
-0.0002
0
0.0002
0 0.02 0.04 0.06 0.08 0.1
Di
sp
lac
em
en
t (
m)
Time (s)
Transformation Begins
Transformation
Ends
63
Figure 3.7 X Velocity along the Width of the SMA Diaphragm
-0.008
-0.006
-0.004
-0.002
0
0.002
0.004
0.006
0.008
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
X
Ve
loc
ity
(m
/s)
Position Along SMA
64
Figure 3.8 Maximum X Velocity of Diaphragm at each Time Step
0
0.001
0.002
0.003
0.004
0.005
0.006
0.007
0 0.02 0.04 0.06 0.08 0.1 0.12
X
Ve
loc
ity
(m
/s)
Time (s)
Transformation Begins Transformation Ends
65
Figure 3.9 Maximum Y Velocity along Width of SMA Diaphragm
-0.045
-0.04
-0.035
-0.03
-0.025
-0.02
-0.015
-0.01
-0.005
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Y
Ve
loc
ity
(m
/s)
Position along Width of SMA Diaphragm
66
Figure 3.10 Maximum Y Velocity of Diaphragm at each Time Step
-0.045
-0.04
-0.035
-0.03
-0.025
-0.02
-0.015
-0.01
-0.005
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
Y
Ve
loc
ity
(m
/s)
Time (s)
Y Velocity at Center of
Diaphragm
Y Velocity at 1/4th Distance of
Diaphragm
67
Chapter 4 Lattice Boltzmann Boundary Conditions
4.1 Boundaries Overview
In general for fluid analysis, boundaries around the domain of interest are solid
walls or a fluid interface. For the solid wall boundary, a velocity is prescribed. At a fluid
interface, a velocity or pressure boundary condition is imposed. Figure 4.1 shows a
channel of flow that is aligned with the horizontal (x) and vertical (y) axes. The walls
located at the top and bottom are impermeable, and they have velocity type boundary
conditions. The fluid inlet and outlet (left and right) have either a velocity or pressure
boundary condition.
The method of dealing with boundaries in this chapter is based upon the bounce-
back method (Zou and He 1997). The non-equilibrium component (from LBM
equilibrium equation) normal to the wall is bounced back. The details of implementing
this procedure are illustrated below.
4.2 Velocity and Pressure Boundaries using Lattice Boltzmann Method
Figure 4.1 shows the three basic regions that are present when doing a fluid solve
with the Lattice Boltzmann method (LBM). Region 1 is used exclusively for a solid wall
while region 2 is used for a fluid boundary. Lastly, region 3 is used for a corner between
walls or a wall and fluid.
68
The conservation of mass and momentum are the three equations used to solve for
the unknowns in each region. Using basic equations presented in Chapter 2, the density
is represented as a function of the density distributions as follows
?
=
=++++++++=
8
0
876543210
i
iffffffffff? (4.1)
The macroscopic flow velocities are also represented as a function of the density
distributions as ?=
?
???? fu where g1873g3364 is the macroscopic velocity vector. Therefore, the
velocity components are
658731 ffffffu ?++??=? (4.2)
874652 ffffffv ???++=? (4.3)
where ?u? is the horizontal velocity and ?v? is the vertical velocity. The corresponding
density distribution functions in the above equations are numbered the same as those
presented in Chapter 2.
4.2.1 Velocity Boundary Conditions
The analysis for the velocity boundary condition in this chapter is restricted to the
horizontally aligned channel shown in Figure 4.1. Each of the four walls, or fraction of
wall, is either solid or fluid. The only difference between region 1 and 2 is that region 2
69
represents the portion of the channel that is fluid. The three different regions shown in
Figure 4.1 are dealt with below.
4.2.1.1 Wall Velocity Boundary (Region 1)
In Figure 4.1, the node on the bottom wall (region 1) is depicted without the
standard nomenclature for the density distributions. The nomenclature for the
distribution functions are changed such that the node shown on the bottom wall is easily
rotated, i.e., rotate 180 degrees and use on the top wall. The new nomenclature for
unknown density distributions is uf and the nomenclature for knowns is kf . The
velocity at the wall is denoted as un or ut for the normal and tangential components,
respectively. For the particular orientation of the node on the bottom wall in region 1, the
newly assigned unknown distributions are 61 ffu = , 22 ffu = , and 53 ffu = . The known
distributions are 31 ffk = , 72 ffk = , 43 ffk = , 84 ffk = , and 15 ffk = . The standardized
normal and tangential velocities for the bottom node are yn uu = and xt uu = ,
respectively. These are prescribed velocities.
In addition to the unknown distributions listed above, ? is also unknown. Using
the conservation of mass and momentum equations yield
421331250 kkuukkuk fffffffff ++++++++=? (4.4)
421315 kkuukkt ffffffu +??+?=? (4.5)
70
421332 kkuukun ffffffu ??++?=? (4.6)
For the bounce-back of the non-equilibrium component,
eqkkequu ffff 3322 ?=? (4.7)
The equilibrium component, g1858g3032g3044, is the same as equation 2.13.
The difference in the equilibrium components is from the one term in the
equilibrium equation in which the velocity is not squared. Hence,
nnneqkequ uuuff ?? 32))3(3(932 =??=? (4.8)
therefore,
nku uff ?3232 += (4.9)
Substitution and algebraic manipulation of the above equations yield the following
explicit equations for the unknowns.
)](2[1 1 423150 kkkkk
n
ffffffu +++++?=? (4.10)
71
nku uff ?3
2
32 += (4.11)
ntkkku uuffff ?? 6
1
2
1)(
2
1
1523 ++??= (4.12)
ntkkku uuffff ?? 6
1
2
1)(
2
1
1541 +??+= (4.13)
4.2.1.2 Fluid Velocity Boundary (Region 2)
In region 2, the unknowns after streaming in LBM are ? , 1uf , 2uf , and 3uf . As
similar to the bottom wall, the three conservation equations and bounce-back are used to
solve for the unknowns. The new assignment of the distribution functions and velocities
for nodes in region 2 are now 51 ffu = , 12 ffu = , 83 ffu = , 21 ffk = , 62 ffk = , 33 ffk = ,
74 ffk = , and 45 ffk = . Similarly, the normal and tangential velocity components for the
fluid inlet are xn uu = and yt uu = .
After assigning the knowns and unknowns, the equations for the unknowns are
identical to the velocity solution for the bottom wall. This greatly simplifies the
computer code needed to handle the different sides (orientations) of the flow channel.
4.2.1.3 Corner Velocity Boundary (Region 3)
After the streaming step in the LBM, the obvious unknowns at the corner node are
1uf , 2uf , 3uf , 4uf , and 5uf . For the lower left corner node specifically, 61 ffu = ,
72
22 ffu = , 53 ffu = , 14 ffu = , and 85 ffu = . Note that if the pressure is specified at the
inlet, then the density is known. The other knowns are 0f , 1kf , 2kf , and 3kf , where
31 ffk = , 72 ffk = , and 43 ffk = . Applying the bounce-back from the distributions that
are normal to the inlet and bottom wall,
11414 )( keqkequku fffff =?+= (4.14)
33232 )( keqkequku fffff =?+= (4.15)
Note that the velocities (ux, uy) at the corner are zero. Using the momentum conservation
equations for both ?x? and ?y? directions,
23 ku ff = (4.16)
51 uu ff = (4.17)
where
)(21 23312405 kukkuuu ffffffff ???????= ? (4.18)
73
The numerical computation for the other corner nodes is identical once the initial uf and
kf assignments are made. Each corner node is rotated relative to the lower left corner
node and then the unknowns properly assigned. As previously stated, this simplifies the
computer code needed for different orientations.
4.2.2 Pressure Boundary Conditions
The analysis for the pressure boundary condition in this chapter is restricted to the
horizontally aligned channel shown in Figure 4.1. Each of the four walls, or fraction of
wall, is solid or fluid. The only difference between region 1 and 2 is that region 2
represents the boundary of the channel that is fluid. The algorithm developed below
allows any of the walls to be fluid with a pressure applied.
4.2.2.1 Fluid Pressure Boundary (Region 2)
For region 2 once again, the assignments for the generic distribution functions and
velocities are the same as that above. In the case of a prescribed pressure on a boundary,
the tangential velocity is zero, i.e., ut = 0. This is the asumption made in the modeling of
the blood pump in this research. All of the blood is assumed to leave ?normal? to the exit
opening. The unknowns now are 1uf , 2uf , 3uf , and nu . The conservation equations are
543213210 kkkkkuuu fffffffff ++++++++=? (4.19)
432321 kkknuuu fffufff +++=++ ? (4.20)
74
425131 kkkkuu ffffff +?+?=? (4.21)
By combining the mass equation and the momentum in the normal direction, the normal
velocity is
?
)](2[1 423510 kkkkk
n
ffffffu +++++?= (4.22)
The unknown distribution, 2uf , is obtained by using bounce-back, namely,
eqkkequu ffff 3322 ?=? (4.23)
Solving for 2uf yields
nku uff ?3
2
32 += (4.24)
Lastly, the other two unknowns are obtained using the momentum equations. They are
nkkku uffff ?61)(21 5141 +??= (4.25)
75
and nkkku uffff ?61)(21 5123 +?+= (4.26)
4.3 Moving Walls and Curved Boundaries for SMA Diaphragm
In the previous section, the moving wall is always located exactly on a grid point.
The bounce-back routine is implemented to solve for both pressure and velocity boundary
conditions. In a general situation, the boundaries of the LBM domain don?t always align
perfectly on structured grids. Additional numerical techniques are needed to resolve
these special cases. Since this research is only a feasibility study, it assumes all
boundaries fall exactly on given nodes.
As shown in Chapter 3, the motion of the SMA diaphragm is well understood
when implementing the phase transformation kinetics. The added complication with the
SMA diaphragm is the variable velocity, or acceleration, during the SMA deformation.
This is accounted for at each time step during the solve.
76
Figure 4.1 Different Boundaries for LBM
Region
1
Fluid
Flow
Region
2
Region
3
Inlet Outlet
1uf 2uf 3uf
1kf 5kf
2kf 4kf3kf
1uf 2uf 3uf
4uf
5uf
1kf
2kf 3kf
1uf
2uf
3uf
1kf
2kf
3kf
4kf 5kf
77
Chapter 5 Hydrodynamics of Micropump
5.1 Hydrodynamic Overview
There are two basic fluid flows to consider during the operation of the SMA
micropump, namely, the medical fluid in the top chamber (as shown in Figure 1.1) and
blood flow in the lower chambers. The modeling of the medical fluid in the upper
chamber is a basic conservation of mass for fluid entering or exiting the check valves.
The motion of the SMA diaphragm dictates this flow. Other than conservation of mass
calculations, the behavior of the medical fluid in the top of the micropump is not
considered.
The blood flow in the lower chambers is fully modeled using the Lattice
Boltzmann Method (LBM). The blood entering the micropump is driven by the human
heart, and is transient in nature. The other transient part of the analysis is the motion of
the SMA diaphragm. The position of the SMA diaphragm is directly dependent upon the
addition or removal of heat, as presented in Chapter 2. Figure 5.1 shows the relative
motion of the diaphragm from the minimum position (1%) to the maximum position
(8%). Chapter 6 addresses all of the thermal aspects of the model and how they impact
the SMA diaphragm. The remainder of this chapter is focused primarily upon the fluid
analysis for the blood flow in the micropump.
78
5.2 Lattice Boltzmann Benchmarking
In order to validate the Lattice Boltzmann Method (LBM), which is used for all of
the fluid analysis, a thorough benchmarking is included. A Womersley flow is chosen as
the benchmark since it is a good test for the transient ability of the LBM method. After
the Womersley benchmark, grid independence is tested for the micropump. Lastly, the
ability to handle non-Newtonian fluid flow is investigated.
5.2.1 Units Example
Before proceeding with the LBM benchmark, a units example for the micropump
is included. Proper selection of units is needed in order to ensure stability of the
numerical algorithm. This units example computes the Reynolds number of the blood
flow in order to provide insight into what type of flow is expected (creep, laminar, or
turbulent).
5.2.1.1 Micropump Physical Dimensions
A schematic of the proposed micropump design is shown in Figure 5.1. The
physical dimensions are 4 cm in length and 1 cm in height. The height of the internal
ramp is 0.7 cm high and 1 cm wide at the base. The blood inlet on the center third of the
left side of the pump is 0.3333 cm. Likewise, the blood exit on the top third on the right
side is 0.3333 cm. The last internal component of the micro pump is the SMA diaphragm
that protrudes into the last half of the pump. The maximum protrusion, located at the
center of the diaphragm, varies from 0.17 cm to 0.35 cm. The minimum distance (0.17
cm) is based upon an SMA strain of 1% while the maximum distance (0.35 cm) is based
79
upon an SMA strain of 8%. These are the operation parameters of the diaphragm
depending upon the heating/cooling stage of the micropump. These details are addressed
in Chapter 6.
5.2.1.2 Blood Properties
The viscosity of the pumped blood varies, depending upon the shear rate. This is
addressed later in the chapter. Suffice it to say that the blood behaves non-Newtonian
and that the maximum viscosity (at zero shear rate) is
g2021g3040g3028g3051 = 1.6 g3030g3040g3118g3046 (5.1)
while the minimum viscosity (at infinite shear rate) is
g2021g3040g3036g3041 = 0.035g3030g3040g3118g3046 (5.2)
Blood is predominantly water and is treated as an incompressible fluid with a density of
water. The velocity of blood through a typical artery in the human body is transient in
nature, with a maximum velocity of 100 cm/s. The average velocity of the blood is 30
cm/s.
Before doing a solve with the lattice Boltzmann method, all parameters for length,
time, and viscosity are converted to lattice units. The first step is to define all relevant
80
reference properties for the micropump. The reference length is chosen as the height of
the micro pump, which is
g1864g3045g3032g3033 = 1 g1855g1865 (5.3)
and the reference velocity is chosen as the average blood velocity, which is
g1873g3045g3032g3033 = 30 g1855g1865/g1871 (5.4)
Knowing the reference length and velocity, the reference time is obviously
g1872g3045g3032g3033 = g3039g3293g3280g3281g3048
g3293g3280g3281
= 0.0333333 g1871 (5.5)
The viscosities for blood are already stated above; thus the minimum and
maximum Reynolds number (Re) are computed for the blood flow. The minimum Re is
when the blood is at maximum viscosity. This is computed as
g1844g1857g3040g3036g3041 = g3048g3293g3280g3281g3039g3293g3280g3281g3092
g3288g3276g3299
= g4672g2871g2868
g3278g3288
g3294 g4673(g2869 g3030g3040)
g2869.g2874 g3030g3040g3118/g3046 = 18.75 (5.6)
And lastly, the maximum Re number is
g1844g1857g3040g3028g3051 = g3048g3293g3280g3281g3039g3293g3280g3281g3092
g3288g3284g3289
= g4672g2871g2868
g3278g3288
g3294 g4673(g2869 g3030g3040)
g2868.g2868g2871g2873 g3030g3040g3118/g3046 = 857 (5.7)
81
so,
18.75 < g1844g1857 < 857 (5.8)
Based upon the Re number, the blood flow is laminar and well behaved since the
maximum value is well below 2000 (typical number for turbulence in internal flow).
5.2.1.3 Dimensionless Variables and Discretization of Space and Time
The time, space, and velocity in physical units are converted to dimensionless
numbers using the equations
g1872g3031 = g3047g3291g3283g3300g3294g3047
g3293g3280g3281
(5.9)
g1876g3031 = g3051g3291g3283g3300g3294g3039
g3293g3280g3281
(5.10)
g1873g3031 = g1873g3043g3035g3052g3046 g3047g3293g3280g3281g3039
g3293g3280g3281
(5.11)
g2021g3031 ? g2869g3019g3032 (5.12)
A typical grid used for this micro pump is 27 grid points for the pump height and 102
grid points for the pump length. Therefore, the space interval is
82
g2012g3051 = g3039g3293g3280g3281g3052 g3034g3045g3036g3031 = g2869(g2870g2875g2879g2869) = 0.03846 (5.13)
For stability (Abdelgawad, 2005), the time interval is typically chosen such that
g2012g3047 ? g2012g3051g2870 (5.14)
For the analysis presented here, the time is chosen as approximately the discrete space
squared, which is
g2012g3047 = 0.002 (5.15)
This implies that the number of iterations required to reach one reference time (for this
example) is
g1840 = g2869g3083
g3295
= 500 (5.16)
5.2.1.4 Convert to Lattice Boltzmann Units
The final manipulation of units is the conversion of everything to lattice units
(lbu). These are listed as
g1873g3039g3029g3048 = g1873g3031 g3083g3299g3083
g3295
(5.17)
83
g2021g3039g3029g3048 = g2869g3019g3032 g3083g3295g3083
g3299g3118
(5.18)
For a final numerical computation example, suppose a velocity is imposed on a boundary
as
g1873g3043g3035g3052g3046 = 15g3030g3040g3046 = 0.5 (dimensig145ng142ess) (5.19)
then
g1873g3039g3029g3048 = 0.5g1499 g3083g3295g3083
g3299
= 0.0260 (5.20)
g2021g3039g3029g3048 = g3083g3295g3083
g3299g3118
g2869
g3019g3032 = 0.0721 (5.21)
The final calculation is for the relaxation time. As stated from Chapter 4,
g2028 = g3092g3287g3277g3296g3030
g3294g3118
+0.5 = 3g2021g3039g3029g3048 +0.5 (5.22)
g2028 = 0.7163 (5.23)
84
5.2.2 Womersley Flow for Benchmarking
In order to provide a robust benchmark for the LBM method, a Womersley flow
is analyzed. A Womersley flow is chosen since it will test the transient accuracy of the
LBM method. Additionally, a Womersley flow is somewhat analogous to the transient
velocity of blood being pumped through the blood micropump. Hence, the transient
solve of the LBM is quite important for the micropump. The LBM model is compared to
a Fourier Series solution of the same flow.
5.2.2.1 Womersley Flow Theory
The Womersley flow solved is a periodic flow (oscillatory pressure) in a 2-D
channel. A sinusoidal pressure boundary condition is applied at the entrance of the
channel. Some of the basic assumptions for the flow are
? Impervious walls
? Incompressible flow
? Fully developed laminar flow
? Isothermal Newtonian flow with constant fluid properties
The general form of the continuity equation is
g3105g3096g3105g3047 + g3105g3105g3051(g2025g1873)+ g3105g3105g3052(g2025g1874)+ g3105g3105g3053(g2025g1875) = 0 (5.24)
With impervious walls and incompressible flow, the continuity equation above reduces to
85
g3105g3105g3051(g2025g1873) = 0 (5.25)
Since g2025 is constant, this requires the axial velocity component (u) to be constant in x
direction. Hence, the velocity varies only with y and time.
g1873 = g1873(g1877g481g1872) (5.26)
The Navier-Stokes equations for 2-D channel flow, with the above assumptions applied,
yields
g3105g3048g3105g3047 = ?g2869g3096 g3105g3043g3105g3051 +g1874 g3105g3118g3048g3105g3052g3118 (5.27)
0 = ?g2869g3096 g3105g3043g3105g3052 (5.28)
The second N-S equation implies that pressure is only a function of x. Therefore,
g1868 = g1868(g1876g481g1872) (5.29)
5.2.2.2 Fourier Series Solution
A Fourier series for the unknowns, u and p, is written as
g1873 = g1873g2868 +uni2211 g1873g3041g1857(g3036g3104g3041g3047)?g3041g2880g2869 (5.30)
86
?g2869g3096 g3105g3043g3105g3051 = g1868g2868 +uni2211 g1868g3041g1857(g3036g3104g3041g3047)?g3041g2880g2869 (5.31)
Substituting these series into the N-S equations above yields
uni2211 g1857(g3036g3104g3041g3047)g3428g3036g3104g3041g3048g3289g2879g3049
g3279g3118g3296g3289
g3279g3300g3118 g2879g3043g3289g3432g2879g3043g3116g2879g3049
g3279g3118g3296g3116
g3279g3300g3118 g2880g2868?g3041g2880g2869 (5.32)
To ensure that this equation is always zero, two equations must be true, namely,
?g1868g2868 ?g1874g3031g3118g3048g3116g3031g3052g3118 = 0 (5.33)
and
g1861g2033g1866g1873g3041 ?g1874g3031g3118g3048g3289g3031g3052g3118 ?g1868g3041 = 0 (5.34)
These ODE?s are integrated and solved.
Since this is a Womersley flow and the pressure is a sinusoidal function, the
pressure is written as
g3105g3043g3105g3051 = g1844g1857g3427g1827g1857g3036g3104g3047g3431 (5.35)
where A is the amplitude and ? is the frequency. An analytical solution for the velocity
is obtained as (He and Luo 1997)
87
g1873(g1877g481g1872) = g1844g1857g3429g1861 g3002g3104g34371?
g3030g3042g3046g3428g3090(g3118g3300g3261g3300g2879g2869)g3432
g3030g3042g3046g3090 g3441g1857
g3036g3104g3047g3433 (5.36)
where
g2019g2870 = ?g1861g2009g2870 (5.37)
and
g2009g2870 = g3013g3300
g3118g3104
g2872g3092 (5.38)
g2021 is the viscosity and g2009 is Womersley number. g1838g3052 is the height of the channel.
5.2.2.3 Lattice Boltzmann Solution
The LBM method developed in the previous chapters is used in order to
benchmark the solution with the Fourier Series solution. The convergence for different
grids is analyzed as well. Implementation of a typical 2-D channel problem is described
below in lattice units for simplicity.
The pressure is applied to the 2-D channel entrance based upon the plane
poiseuille flow equation. The Poiseuille equation is
88
g1873g3040g3028g3051 = (g3017g3117g2879g3017g3118)g3013g3300
g3118
g2876g3092g3013g3299 (5.39)
Where g1838g3051 and g1838g3052 are lattice grids. All of the key input parameters for the model are
listed in lattice units in Table 5.1
Simulation Parameters Lattice Units
Grid Size 21 x 21 and scaled up to 161 x 161
Viscosity 0.03925
g2028 (relaxation factor) 0.6178
Pressure Gradient (dP/dx) 0.001
g1873g3040g3028g3051 0.0668
Period (T) 1000
Womersley Number (?) 4.2
Table 5.1 Parameters for Womersley Flow Solution Listed in Lattice Units
Figure 5.2 shows the results of the LBM and Fourier Series solution for a grid of
41 x 41. The axial velocity is plotted for three different points in the sinusoidal pressure
cycle, namely, T/4, T/2, and 3T/4. The results of the two different methods agree quite
well. The relative accuracy for both methods is shown in Figures 5.3 and 5.4.
89
5.2.3 Convergence Analysis (Grid Independence)
From the Womersley benchmark problem, the LBM method agrees very well with
the Fourier Series solution. Additional grid refinement results are performed for the
micropump itself. To test the relative accuracy of varying grids, the axial velocity from
two cross-sections of the micropump is plotted as a function of grid size. The cross-
sections analyzed are the center of the first chamber and second chamber. These results
are shown in Figures 5.5 and 5.6. As expected, the solutions are well behaved and
converge with increasing grid size.
5.3 Blood Characterization
With the newly proposed shape memory alloy micropump, blood is being utilized
as a cooling fluid. In order to determine the feasibility of this design, the flow
characteristics of blood have to be determined. Specifically, the physical and rheological
properties of blood need to be understood.
5.3.1 General Properties and Composition of Whole Blood
The primary constituents of whole blood are plasma and red blood cells (RBC).
The white blood cells and platelets constitute a significantly less volume fraction than
RBC?s. The volume fraction that the RBC?s occupy is referred to as the hematocrit. For
typical men, the hematocrit is in the range of 40-50% (Abdelgawad, 2005). The plasma
is primarily water (90%) and protein (7%) (Abdelgawad, 2005). The plasma itself
behaves in a Newtonian manner and has a typical viscosity that is approximately 1.2 x 10-
3 Pa-s. The other blood constituents are non-Newtonian and are addressed below.
90
5.3.2 Non-Newtonian Nature of Blood Viscosity
The flow of whole blood does not behave in a Newtonian manner. The blood
flow regime is divided into three primary regions that are dependent upon strain rate. At
very low strain rates, blood has a constant viscosity. On the other extreme, high strain
rate, the blood has a constant viscosity, but it is significantly less than the viscosity in the
low strain rate region. At the intermediate strain rate region, the viscosity varies as a
function of the strain rate. Hence, blood behaves non-Newtonian. When a fluid
viscosity increases as a result of increasing strain rate, it is considered as ?shear-
thickening?. Conversely, if the viscosity decreases as a result of increasing strain rate, it
is considered as ?shear-thinning?. The behavior of blood falls into the latter category. In
Figure 5.7, shown for qualitative purpose only, the three regions mentioned above are
illustrated. The non-Newtonian region results from the RBC?s layering under shear
stress.
5.3.2.1 Power-Law Approximation
The simplest model to represent non-Newtonian behavior is the Power Law
(Abdelgawad, 2005). The resulting equation to relate viscosity (g2020) to strain rate (g2011g4662) is
g586 = g586g2868uni007Cg576g4662uni007C(g2924g2879g2869) (5.40)
where strain rate in tensor notation for two dimensions is
91
g576g4662 = g2869g2870g4672g2986g2931g3209g2986g2934g3210 + g2986g2931g3210g2986g2934g3209g4673 (5.41)
If the exponent n is set equal to unity, Newtonian behavior is recovered and g2020 = g2020g2868. g2020g2868
in the above equation is referred to as the consistency index. If n > 1, the flow is shear-
thickening. Conversely, if n < 1, the flow is shear-thinning. As previously mentioned,
blood is Newtonian at very low and very high strain rates while shear-thinning at
intermediate strain rates.
To account for the non-Newtonian fluid, the LBM adjusts the viscosity at each
lattice site based upon the strain rate at that node. At each time step in the LBM scheme,
the strain rates are calculated and subsequently the viscosities. The strain rate is
computed using equation 5.41 where (Abdelgawad, 2005)
g3105g3048g3328
g3105g3051g3329 =
g2869
g2874uni2211 g1857g3036g3081g1873g3080(g1876 +g1857g3036g543g1876)
g2877
g3036g2880g2869 (5.42)
Once the viscosity is calculated based upon the strain rate and power-law, the LBM
relaxation time is calculated using this equation from Chapter 2
g2028 = 3g2021 +.5 (5.43)
This time marching and viscosity adjusting scheme is continued throughout the solution.
The benchmarking of the LBM code is done for a 2-D channel with a height of H.
The entire 2-D grid is composed of 32 x 32 lattice nodes. Periodic boundary (inlet
conditions match exit conditions) is applied at the channel input and output while non-
92
slip boundaries are applied at the upper and lower channel walls. All of the results for
velocity are compared to the exact channel analytical solution (Bird) given as
g1873(g1877) = g1873g3040g3028g3051 g46821?g4672g2870uni007Cg3052uni007Cg3009 g4673
g3289g3126g3117
g3289 g4683 (5.44)
where
g1873g3040g3028g3051 = g4672 g3041g3041g2878g2869g4673g4672g3007g3091
g3116
g4673
g3117
g3289 g4672g3009
g2870g4673
g3289g3126g3117
g3289 (5.45)
The consistency index, g2020g2868, is chosen to be the Newtonian viscosity that corresponds to a
relaxation factor of g2028 = 0.75. In lattice units, this Newtonian viscosity is 0.083333. ?F?
is the pressure gradient needed to generate a final velocity, g1873g3033 = 0.01, for Newtonian
flow in a 2-D channel with a viscosity of g2020g2868. In equation form,
g1832 = g3031g3017g3031g3051 = g2876g3091g3116g3048g3281g3009g3118 (5.46)
Figure 5.8 shows a typical strain rate for flow in a 2-D channel. As expected, the
maximum shear occurs near the wall and decreases to a minimum at the center of the
channel. The corresponding viscosities for shear-thickening and shear-thinning fluids are
shown in Figures 5.9 and 5.10. The steady-state velocity profiles for non-Newtonian
fluids are compared to the exact analytical solutions in Figures 5.11 and 5.12. In Figure
5.13, these are normalized against the maximum velocity at the center of the channel.
93
The Power-Law method is simple and useful for benchmarking the LBM code for
non-Newtonian flow. The Power-Law yields a linear relationship between viscosity and
strain rate on a log-log scale. The straight line has a slope equal to ?n?. However, this
model is too simplistic for the complex nature of blood flow. Other blood models exist
that improve upon the Power-Law method.
5.3.2.2 Carreau-Yasuda Blood Model
Perhaps the most highly regarded model for blood flow is the Carreau-Yasuda
(CY) model (Abdelgawad, 2005). The governing equation for the viscosity is
g2020(g2011g4662) = g2020g2998 +(g2020g2868 ?g2020g2998)(1+(g2019g2011g4662)g3028)g3289g3127g3117g3276 (5.47)
The parameters for the CY blood model are shown in Table 5.2. Using these equation
parameters, the dynamic viscosity output (?) units are ?Pa-s?. As mentioned in the
previous section, the shear rate is numerically calculated from
g2011g4662 = g2869g2870g4672g3105g3048g3328g3105g3051g3329 + g3105g3048g3329g3105g3051g3328g4673 (5.48)
The viscosity as a function of shear rate is shown graphically in Figure 5.14. It is evident
that the CY blood model predicts shear-thinning behavior for the blood.
94
CY Blood Model
Parameters
SI units
g2020g2998 0.0035 (Pa-s)
g2020g2868 0.16 (Pa-s)
g2019 8.2
a 0.64
n 0.2128
Table 5.2 CY Blood Model Parameters
5.4 Boundary Implementation
All of the boundaries are numerically implemented using the equations developed
in Chapter 4. All walls have no slip while the fluid boundaries are either velocity or
pressure. These fluid boundaries are specified below.
5.4.1 Inlet Velocity
The micropump schematic is shown in Chapter 1. The blood flow inlet on the left
side of the micropump is oscillatory in nature due to the repetitious pumping of the
human heart. The simplest approach to modeling the blood flow is to assume a
sinusoidal nature of the flow. The maximum pressure corresponds to the systolic blood
pressure while the minimum pressure corresponds to the diastolic pressure. However,
this assumption of a simple sinusoidal flow does not sufficiently capture the dynamic
interaction of the heart and arteries.
95
Experimental data is available for the flow of blood in typical size arteries in the
human body. Figure 5.15 shows the velocity profile of blood being pumped through an
artery of size 2.5 mm. (Abdelgawad, 2005). By knowing the transient nature of the blood
flow, the velocity boundary condition at the entrance of the micropump is implemented.
For each time step in the LBM solve, the lattice time and physical time are computed so
that the velocity at the pump entrance is adjusted accordingly. The blood entrance
channel on the left side of the micropump is the same size as typical arteries (2.5 mm).
In order to extract the velocity data from Figure 5.15 numerically, the region
spanning one second of time is segmented into four sections (a-b, b-c, c-d, d-e). Each
section is then numerically modeled with a polynomial of varying order. The equations
to represent each of the four sections are shown below in Table 5.3. The variable ?x? is
the time, which varies from zero up to one second in physical units. The inlet velocity
profile is assumed to be that of plug flow in which the velocity across the opening is
constant.
Section of Blood
Velocity Profile
Polynomial (x is seconds)
a-b (0 - .12 s) 87.5g1876g2870 ?3.1g1876 +0.1g8916
b-c (.12 - .28 s) 14g89174g1876g2872 ?12448g1876g2871 +3826g1876g2870 ?518g1876 +26.56
c-d (.28 - .36 s) ?87.5g1876g2870 +56.5g1876 ?8.72
d-e (.36 ? 1 s) 0.g89155g1876g2871 ?1.76g891g1876g2870 +0.856g1876 +0.157
Table 5.3 Curve-Fit Equations for Blood Velocity Profile
96
5.4.2 Outlet Pressure
On the exit boundary (right side of micropump), the tangential and normal
velocities are unknown. Using all of the conservation equations developed in Chapter 4,
a pressure boundary condition is applied. At the exit, the pressure is arbitrarily assigned
as unity in the lattice units.
5.4.3 All Other non-Fluid Boundaries
For all of the solid walls (no fluid flow), the 2nd order boundary conditions
developed in Chapter 4 are applied. All of the non-fluid nodes located within the walls of
the micropump are handled using the rather simple bounce-back method. This includes
the SMA diaphragm and the ramp in the center of the micropump.
5.5 Numerical Results
At each time step of the transient analysis for the fluid, both viscosity and velocity
are calculated. They are mutually dependent upon each other since viscosity is dependent
upon strain rate, and consequently, the velocity gradient. Additionally, the fluid
modeling in this chapter accounts for the protrusion of the SMA diaphragm into the
second chamber of the micropump. The strain state of the SMA diaphragm is dependent
upon the thermal model. Since the results here are only hydrodynamics, strain state of
the diaphragm is arbitrarily specified when appropriate. Lastly, the location of the SMA
diaphragm is approximated by the existing x/y grid. Since this is a feasibility study, the
SMA outline will be treated in a stair-step fashion. In future work, interpolating the SMA
position between grids may be appropriate.
97
5.5.1 Blood Flow Velocity
With all of the boundary conditions applied to the pump schematic shown in
Figure 1.1, the hydrodynamics of the blood micropump are numerically solved. Figure
5.16 is a contour plot for the x-component velocity of the blood flow. The lighter regions
represent higher velocities while the black regions represent zero velocity. A similar
contour plot is shown in Figure 5.17 for the y-component of blood velocity. For
simplicity of generating data, the inlet velocity is kept constant at 0.035 lattice units. The
results behave exactly as expected for a uniform velocity across the inlet. The flow
profile at the exit is parabolic and the maximum flow velocity is about 1.5 times greater
than the velocity input.
Figures 5.18 and 5.19 show the streamlines for blood flow through the
micropump in two different configurations. In Figure 5.18, the strain of the SMA is at its
minimum operating point, 1%. In Figure 5.19, the strain is at the maximum operating
point, 8%.
Looking at the blood flow in more detail, results are shown at the pump cross-
section at the center of the SMA diaphragm. This represents the location for the greatest
penetration of the SMA diaphragm into the blood stream. Velocity profiles are shown for
three diaphragm states; strain at 2%, 5%, and 8%. In Figure 5.20, the axial-component of
the flow shows a flow pattern that is typical for internal flow. Finally, the transverse-
component of the flow is shown for the three different strain states in Figure 5.21.
Using the regression curves from Table 5.3, the velocity at the inlet of blood flow
is computed based upon the physical time, varying between 0 and 1 seconds. Figure 5.22
shows the peak velocity as a function of time for two locations along the pump,
98
specifically, the top of the ramp and the outlet. The transient velocity very closely
replicates the shape of the velocity input profile. The maximum input of the blood flow
at the inlet is barely over 1 m/s, and this is constant over the entire inlet opening.
Therefore, the maximum velocity at the top of the ramp and outlet are about 1.5 m/s (1.5
times inlet velocity).
To validate conservation of mass, a constant input velocity of 0.18 m/s is
uniformly applied to the inlet. Figure 5.23 shows the axial velocity profile for the inlet,
top of ramp, and exit. With the assumption of incompressible flow and steady-state, the
area under the curves for the different locations is equal due to conservation of mass.
The area under these velocity curves is computed using trapezoidal integration and mass
is conserved as expected. Note that the opening at the top of ramp is slightly larger
than the exit, hence, a lower peak velocity.
5.5.2 Viscosity
The strain rate needed for the CY blood model is shown in Figure 5.24. The
largest strain rates are in locations as expected. In the high velocity areas near walls and
corners, the strain rate is the maximum. The viscosity is shown in the contour plot in
Figure 5.25. The ?lighter? shades of gray indicate the areas of maximum viscosity.
Since the viscosity is directly dependent upon the strain rate, these results are as
expected. The modified relaxation factor (Tau) is computed for each local viscosity and
is shown in Figure 5.26 (Wang and Bernsdorf 2009).
99
5.5.3 Hydrodynamic Boundary Layer on SMA Diaphragm
Since the next chapter addresses the thermal modeling of the SMA diaphragm, the
hydrodynamic boundary layer near the SMA is important. In order to get an idea of what
to expect for the boundary layer, three different areas of the diaphragm are analyzed for
four different diaphragm strains. The three locations of interest are the midpoint of the
SMA length, 1st quarter of the SMA length, and 3rd quarter of the SMA length. In
Figures 5.27 through 5.32, both the axial and transverse components are graphed
according to location and SMA strain. The boundary layer does vary, but is quite similar
in several areas. The transverse velocity at the midpoint seems a little sporadic for the
different strains, but this is due to the velocity being almost zero at that location. Hence,
they should be negligible in the thermal analysis.
Figure 5.1 Schematic of Blood Micropump with
Blood Flow
Blood Flow
Blood Flow
100
Penetrating SMA Diaphragm
SMA 2%
SMA 8%
SMA 5%
Figure 5.2 Comparison of Womersley Flow (Diamonds
101
-Fourier Series, Solid Line
-LBM)
102
Figure 5.3 Lattice Boltzmann Relative Velocity Error for Varying Square Grids
0.00E+00
2.00E-04
4.00E-04
6.00E-04
8.00E-04
1.00E-03
1.20E-03
1.40E-03
1.60E-03
0 20 40 60 80 100
Re
lat
ive
E
rror
Grid Size
103
Figure 5.4 Fourier Series Relative Velocity Error for Varying Square Grids
0.00E+00
5.00E-05
1.00E-04
1.50E-04
2.00E-04
2.50E-04
3.00E-04
3.50E-04
4.00E-04
0 50 100 150 200
Re
lat
ive
E
rror
Grid Size
104
Figure 5.5 Cross-sections at Center of First Chamber of Micropump
0
0.1
0.2
0.3
0.4
0.5
0.6
0 0.2 0.4 0.6 0.8 1 1.2
Axi
al
Ve
loc
ity
(lat
tic
e u
nit
s)
Distance Across Micropump
20x80
25x100
40x160
50x200
105
Figure 5.6 Cross-sections at Center of Last Chamber of Micropump
0
0.1
0.2
0.3
0.4
0.5
0.6
0 0.2 0.4 0.6 0.8 1 1.2
Axi
al
Ve
loc
ity
(lat
tic
e u
nit
s)
Distance Across Micropump
20x80
25x100
40x160
50x200
Figure 5.
106
7 Rheological Behavior of Whole Blood
107
Figure 5.8 Strain Rate for 2-D Channel Flow
0
0.0005
0.001
0.0015
0.002
0.0025
0.003
0.0035
0.004
0 5 10 15 20 25 30 35
Typical Strain Rate, n = 1.2
108
Figure 5.9 Viscosity Profile for Shear-Thickening Fluid
0.01
0.012
0.014
0.016
0.018
0.02
0.022
0.024
0.026
0.028
0 5 10 15 20 25 30 35
Viscosity Profile, n = 1.2
109
Figure 5.10 Viscosity Profile for Shear-Thinning Fluid
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0 5 10 15 20 25 30 35
Viscosity Profile, n = 0.8
110
Figure 5.11 Velocity Profile for Shear Thickening Fluid
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0 10 20 30 40
Velocity Profile, n = 1.2
u numerical
u exact
111
Figure 5.12 Velocity Profile for Shear-Thinning Fluid
0
0.0002
0.0004
0.0006
0.0008
0.001
0.0012
0.0014
0.0016
0.0018
0 10 20 30 40
Velocity Profile, n = 0.8
u numercal
u exact
112
Figure 5.13 Velocity Profiles Normalized for Maximum Velocity
0
0.2
0.4
0.6
0.8
1
1.2
0 10 20 30 40
Normalized Velocity Profiles
n = 1.2
n = 0.8
113
Figure 5.14 Viscosity Predictions from CY Blood Model
0.01
0.1
1
10
0.01 0.1 1 10 100
Kin
em
ati
c V
isc
os
ity
(c
m^
2/
s)
Shear Rate (1/s)
Blood Viscosity from CY Blood Model
114
Figure 5.15 Typical Velocity Profile for Pumped Blood
0
0.2
0.4
0.6
0.8
1
1.2
0 0.2 0.4 0.6 0.8 1
Ve
loc
ity
(m
/s)
Time (s)
Blood Velocity in 2.5 mm Artery
a
b
c
d
e
115
Fig
ure
5.
16
A
xia
l V
elo
cit
y (
lat
tic
e u
nit
s)
thr
ou
gh
ou
t B
loo
d M
icr
op
um
p
116
Fig
ure
5.
17
Tr
an
sve
rse
V
elo
cit
y (
lat
tic
e u
nit
s)
thr
ou
gh
ou
t B
loo
d M
icr
op
um
p
Figure 5.18 Streamlines of Blood Flow Through Micropump with SMA at 2% Strain
117
Figure 5.19 Streamlines of Blood Flow Through Micropump with SMA at 8% Strain
118
119
Figure 5.20 Axial Velocity of Blood for Pump Cross-Section at Diaphragm Center
1, 00
0.05
0.1
0.15
0.2
0.25
0.3
0 0.2 0.4 0.6 0.8 1
Bl
ood
Ve
loc
ity
(m
/s)
Transverse Location Across Pump
strain=2%
strain=5%
strain=8%
120
Figure 5.21 Transverse Velocity of Blood for Pump Cross-Section at Diaphragm Center
-0.005
-0.004
-0.003
-0.002
-0.001
0
0.001
0.002
0 0.2 0.4 0.6 0.8 1
Bl
ood
Ve
loc
ity
(m
/s)
Transverse Location Across Pump
strain=2%
strain=5%
strain=8%
121
Figure 5.22 Maximum Blood Velocity at top of Ramp and Exit of Micropump
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
0 0.2 0.4 0.6 0.8 1 1.2
Ve
loc
ity (
m/
s)
Time (s)
Max Exit Vel (m/s)
Max Ramp Vel (m/s)
122
Figure 5.23 Velocity Profiles at end of One Second of Blood Flow
-0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
0 2 4 6 8 10 12
Inlet Velocity (m/s)
Top Ramp Velocity (m/s)
Exit Velocity (m/s)
123
Figure 5.24 Strain Rate of Blood in Micropump
124
Figure 5.25 Viscosity of Blood Based upon CY Model
125
Fig
ure
5.
26
M
od
ifie
d R
ela
xa
tio
n F
act
or
(T
au
) th
rou
gh
ou
t M
icr
op
um
p
126
Figure 5.27 Axial Velocity Boundary Layer at First Quarter of SMA Length
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0 2 4 6 8 10 12
Axi
al
Ve
loc
ity
(lat
tic
e u
nit
s)
Boundary Layer from Edge of SMA
SMA 0%
SMA 3%
SMA 5%
SMA 8%
SMAWall
127
Figure 5.28 Transverse Velocity Boundary Layer at First Quarter of SMA Length
-0.09
-0.08
-0.07
-0.06
-0.05
-0.04
-0.03
-0.02
-0.01
0
0 2 4 6 8 10 12
Tr
an
sve
rse
Ve
loc
ity
(lat
tic
e u
nit
s)
Boundary Layer from Edge of SMA
SMA 0%
SMA 3%
SMA 5%
SMA 8%
SMAWall
128
Figure 5.29 Axial Velocity Boundary Layer at Midpoint of SMA Length
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0 2 4 6 8 10 12
Axi
al
Ve
loc
ity
(lat
tic
e u
nit
s)
Boundary Layer from Edge of SMA
SMA 0%
SMA 3%
SMA 5%
SMA 8%
SMA Wall
129
Figure 5.30 Transverse Velocity Boundary Layer at Midpoint of SMA Length
-0.008
-0.007
-0.006
-0.005
-0.004
-0.003
-0.002
-0.001
0
0 2 4 6 8 10 12
Tr
an
sve
rse
Ve
loc
ity
(lat
tic
e u
nit
s)
Boundary Layer from Edge of SMA
SMA 0%
SMA 3%
SMA 5%
SMA 8%
SMA Wall
130
Figure 5.31 Axial Velocity Boundary Layer at 3rd Quarter of SMA Length
0
0.02
0.04
0.06
0.08
0.1
0.12
0 2 4 6 8 10 12
Axi
al
Ve
loc
ity
(lat
tic
e u
nit
s)
Boundary Layer from Edge of SMA
SMA 0%
SMA 3%
SMA 5%
SMA 8%
SMA Wall
131
Figure 5.32 Transverse Velocity Boundary Layer at 3rd Quarter of SMA Length
0
0.005
0.01
0.015
0.02
0.025
0.03
0 2 4 6 8 10 12
Tr
an
sve
rse
Ve
loc
ity
(lat
tic
e u
nit
s)
Boundary Layer from Edge of SMA
SMA 0%
SMA 3%
SMA 5%
SMA 8%
SMA Wall
132
Chapter 6 Thermal Results of Micropump
6.1 Thermal Overview
In the previous chapter, the hydrodynamics for the blood flow are presented.
Knowing the blood flow velocities throughout the SMA micropump, it is possible to
solve the energy equation for both conduction and convection from the SMA diaphragm
to the blood. The thermal results are solved using Finite Volumes. The rest of this
chapter presents those results.
6.2 General Approach for Thermal Solver
The SMA diaphragm in the micropump is a moving boundary that protrudes into
the cavity where the blood flows. The strain of the diaphragm when completely heated is
1%, which implies that the protrusion of the diaphragm into the blood cavity is minimal.
During the cooling phase of the thermal cycle, the strain increases until 8% is achieved.
At this point, the crystalline structure of the SMA is completely reversed.
The thermal model will only include the right chamber of the micropump
(chamber with SMA at top). The entrance for the thermal model is the opening at the top
of the ramp. The exit for the thermal model is the right side of the micropump. The
blood flow velocity in the SMA chamber is imported from the fluid solver. The blood
velocity is high enough such that the entrance boundary at the top of the ramp is treated
as having a constant temperature of 37 degrees C. The temperature at the micropump
133
exit is numerically calculated using the conduction/convection from the finite volume
solver.
The micropump domain for the thermal solver has the same discretization as the
fluid solver. The time step for the Lattice Boltzmann Method (LBM) fluid solver is
determined by the grid size and other parameters. Once the fluid and thermal time steps
are reconciled, the fluid/thermal solution is obtained. Based on preliminary analysis, the
cooling cycle takes about one second for the average blood flow velocity. Knowing this,
it is assumed that the thermal modeling time step will set as 0.05 seconds. Since 15000
steps correspond to one second for the fluid solver, the thermal solver is performed every
750 lattice time steps (0.05 s).
The SMA diaphragm is divided into discrete elements that match the grids of the
micropump domain. During the cooling, nodes for the diaphragm are shifted down for
each time step as dictated by the derivations in chapter 3. The adjacent fluid nodes are
shifted down with the diaphragm. The SMA diaphragm is marked as a solid, and the
bounce-back method is used to model the obstruction of the diaphragm. The temperatures
of the blood in close proximity to the SMA diaphragm are shifted down such that a
complete thermal solve is permitted for both convection and conduction. The substrate
that the Nitinol SMA diaphragm material is attached to is adiabatic by choice of material.
Using an adiabatic material for the substrate eliminates the need to do a thermal analysis
for the upper chamber of the micropump.
6.3 Introduction of Finite Volume Method for Thermal Analysis
The thermal modeling of the heat exchange between the SMA diaphragm and the
blood flow in the micropump is accomplished using the finite volume method (Patankar).
134
This method is rather straightforward and uses energy balances on a small volume of
domain of interest. The solver is transient and incorporates both conduction and
convection. Radiation is negligible and ignored.
The transient energy equation with internal heat generation is
g2025g1829g3043 g3105g3021g3105g3047 + g3105g3105g3051
g3285
g3435g2025g1829g3043g1873g3037g1846g3439 = g3105g3105g3051
g3285
g3436g1863 g3105g3021g3105g3051
g3285
g3440+g1843 (6.1)
where ?Q? is internal energy generation. Using a finite volume (control volume)
approach to solving the energy equation, the discretization for the finite volume depicted
around point (i,j) is graphically illustrated below in Figure 6.1. The bold arrows indicate
the heat flow in and out of the control volume. Doing an energy balance on each of the
four sides of the control volume yields the equation
g2025g3036g481g3037g1829g3043g3284g481g3285 g3436g3021g3284g481g3285
g3287g3126g3117g2879g3021
g3284g481g3285g3287
?g3047 g3440?g1876?g1877+g1863g3036g2878g3117g3118g481g3037 g3436
g3021g3284g481g3285g3287g3126g3117g2879g3021g3284g3126g3117g481g3285g3287g3126g3117
?g3051 g3440?g1877+g1863g3036g481g3037g2878g3117g3118 g3436
g3021g3284g481g3285g3287g3126g3117g2879g3021g3284g481g3285g3126g3117g3287g3126g3117
?g3052 g3440?g1876 +
g2025g1873g3036g2878g3117
g3118
g1829g3043g1846g3036g2878g3117
g3118
g3039g2878g2869?g1877+g2025g1874
g3037g2878g3117g3118g1829g3043g1846g3037g2878g3117g3118
g3039g2878g2869 = g1863
g3036g2879g3117g3118g481g3037 g3436
g3021g3284g3127g3117g481g3285g3287g3126g3117 g2879g3021g3284g481g3285g3287g3126g3117
?g3051 g3440?g1877+g1863g3036g481g3037g2879g3117g3118 g3436
g3021g3284g481g3285g3127g3117g3287g3126g3117 g2879g3021g3284g481g3285g3287g3126g3117
?g3052 g3440?g1876 +
g2025g1873g3036g2879g3117
g3118
g1829g3043g1846g3036g2879g3117
g3118
g3039g2878g2869?g1877+g2025g1874
g3037g2879g3117g3118g1829g3043g1846g3037g2879g3117g3118
g3039g2878g2869 +g1843 (6.2)
In order to simplify the above equation, a cardinal direction system as shown in Figure
6.2 is utilized. The conduction term coefficients after dividing by the specific heat (g1829g3043)
are
135
g1830g3032 = g3038g3280g3004
g3291
?g3052
?g3051, g1830g3050 =
g3038g3298
g3004g3291
?g3052
?g3051, g1830g3046 =
g3038g3294
g3004g3291
?g3051
?g3052, g1830g3041 =
g3038g3289
g3004g3291
?g3051
?g3052 (6.3)
while the convection term coefficients are
g1832g3032 = g2025g1873g3032?g1877, g1832g3050 = g2025g1873g3050?g1877, g1832g3046 = g2025g1874g3046?g1877, g1832g3041 = g2025g1874g3041?g1877 (6.4)
Algebraic manipulation of the above equations in conjunction with the upwind
differencing scheme (UDS) yields
g1853g3017g1846g3017 = g1853g3006g1846g3006 +g1853g3024g1846g3024 +g1853g3015g1846g3015 +g1853g3020g1846g3020 +g1854 (6.5)
where
g1853g3006 = g1830g3032 + g1764?g1832g3032g4810g1765
g1853g3024 = g1830g3050 + g1764g1832g3050g4810g1765
g1853g3015 = g1830g3041 + g1764?g1832g3041g4810g1765 (6.6)
g1853g3020 = g1830g3046 + g1764g1832g3046g4810g1765
g1853g3043g2868 = g3096?g3051?g3052?g3047
g1854 = g1853g3043g2868g1846g3043g2868 +g1843
g1853g3017 = g1853g3006 +g1853g3024 +g1853g3015 +g1853g3020 +g1853g3043g2868
The g1764g481g1765 indicates to take the maximum value. These equations are iterated upon until
converged to the numerical accuracy specified.
136
6.4 Energy Analysis of SMA Diaphragm
For each cycle of operation of the SMA micropump, heat has to be added and
subsequently removed from the SMA diaphragm. Addition of heat is easily
accomplished via resistive heating. For this analysis, it is assumed that there is ample
electrical power to heat the SMA in a timely manner suitable for the pump operation.
The limiting factor for micropump operation frequency is anticipated to be the cooling of
the SMA diaphragm. Hence, the modeling focuses on the cooling required for the
austenite-martensite phase transformation.
The nominal temperature for blood in a human body is 37 C. Since the blood
serves as a cooling agent, the SMA transformation temperatures need to be sufficiently
higher than the blood temperature. The operating temperature for the SMA diaphragm is
controlled by altering the Nickel-Titanium ratio in the Nitinol. For this analysis, the
transformation temperatures of the Nitinol SMA are shown in Table 6.1. The
nomenclature for the transformation temperatures are explained in Chapter 2. The
Cosine Model is used for the transformation computations, and it is also presented in
Chapter 2.
There are two different types of energy removal during the cooling of the Nitinol
SMA, namely, the typical internal energy change of the austenite and then the latent heat
for the phase change from austenite to martensite.
137
Transformation
Parameters
Temperature (C)
g1846g3040g3033 50
g1846g3040g3046 60
g1846g3028g3046 70
g1846g3028g3033 80
Table 6.1 SMA Transformation Temperatures
The reverse is also true for heating. The internal energy change for the high-
temperature austenite is governed by the equation
g1847 = g1865g1829g3043?g1846 (6.7)
Once the phase transformation begins, the latent heat energy needed is
g1831 = ?g1865 (6.8)
where ?h? is the latent heat and ?m? is the mass of the diaphragm. With the mass of the
SMA diaphragm being known, the internal energy required while cooling from austenite-
finish temperature (g1846g3028g3033) to martensite-start (g1846g3040g3046) is readily calculated using equation 6.1.
Likewise, the energy due to latent heat is directly dependent upon the mass of the
diaphragm. As discussed in previous chapters, the dimensions of the micropump are 4
cm in length and 1 cm in height. The SMA diaphragm section of the micropump design
138
is 2 cm in length. Prior to doing an energy analysis, the width and thickness of the SMA
diaphragm is needed. It is assumed that the width of the SMA is equal to its length (2
cm) while the thickness is 0.025 mm. Now the mass is calculated since the volume and
density are specified. Knowing the mass, both ?U? and ?E? are straightforwardly
computed. The numerical values of these quantities are shown in Table 6.2.
SMA Model
Parameters
Values
SMA Mass (2 cm x 2 cm x .025 mm)
0.00006450 kg
Internal Energy (g1846g3028g3033 ?g1846g3040g3046)
0.3225 J
Latent Heat Energy (g1846g3040g3046 ?g1846g3040g3033)
1.935 J
Total Energy (g1846g3028g3033 ?g1846g3040g3033)
2.2575 J
Table 6.2 Energy Requirements during Cool-Down of SMA Diaphragm
6.5 Thermal Flux Analysis for Cooling of SMA Diaphragm
A brief analysis is performed to compare the performance of a micropump in
which the blood is not flowing (conduction only) or the blood is flowing (conduction and
convection). Since this is a preliminary analysis, the blood flow is assumed to have a
constant flow of 0.5 m/s. Later, numerical analysis models the true dynamics of the heart
139
and blood flow. A simple energy balance on the diaphragm is used such that the flux can
be integrated over time. This accounts for temperature change of the SMA diaphragm.
The energy balance accounts for the latent heat during the phase change.
Figure 6.3 shows a typical node placed on the upper SMA wall of the micropump.
The flux at node (i,j) is solved using the energy equation in conjunction with an
exponential differencing scheme. The details of the derivation are not shown (Patankar).
The equation for the flux at the SMA wall is
g1869g3020g3014g3002g4593g4593 = g3038?g3051g4674g1846g3036g481g3037 g4676?g3052g2870?g3051g4672g3017g3254g2915g2934g2926(g3017g3254)g2915g2934g2926(g3017
g3254)g2879g2869
+ g3017g3272g2915g2934g2926(g3017
g3272)g2879g2869
g4673+ ?g3051?g3052 g3017g3268g2915g2934g2926(g3017g3268)g2915g2934g2926(g3017
g3268)g2879g2869
g4677g4675 (6.9)
? g1863?g1876g3428g1846g3036g2878g2869g481g3037 g3420 ?g18772?g1876 g1842g3006eg154g146(g1842
g3006)?1
g3424+g1846g3036g2879g2869g481g3037 g3420 ?g18772?g1876 g1842g3024eg154g146 (g1842g3024)eg154g146(g1842
g3024)?1
g3424+g1846g3036g481g3037g2878g2869g3420?g1876?g1877 g1842g3020eg154g146(g1842
g3020)?1
g3424g3432
where the Peclet number (P) is computed using the ratio of the convection and
conduction terms.
6.5.1 Thermal Flux Analysis with Conduction Only (SMA to Blood)
Initially, a flux analysis is performed for a micropump that contains blood with
zero velocity. The thermal analysis begins for the cool-down period of the SMA after it
is resistively heated to the austenite-complete temperature (g1846g3028g3033). All of the
transformation temperatures are shown in Table 6.1.
During the initial cooling stage, the large temperature difference between the
SMA and blood creates a considerable heat flux. As the temperature of the blood
increases, the heat flux drastically drops during the first second of cooling. After
140
approximately one second, the latent cooling begins but the heat flux is relatively low.
Hence, the conduction model alone is not sufficient to cool the SMA diaphragm.
6.5.2 Thermal Flux Analysis with Conduction and Convection (SMA to Blood)
The next step of flux analysis is the inclusion of convection with the conduction.
The blood flow through the micropump provides a steady stream of fluid at a constant
temperature, 37 C, that is significantly lower than the heated SMA diaphragm. The
energy transferred from the SMA diaphragm to the blood is again calculated using
equation 6.9. The total energy being transported during the cool-down phase is obtained
by integrating the heat flux over time. As anticipated, the convection component of the
energy equation greatly enhances the heat transfer from the SMA diaphragm to the blood.
Figure 6.4 shows the comparison of the thermal flux when the convection is included, or
not included.
6.6 Temperature Profile along SMA Diaphragm
Since the strain of the SMA diaphragm determines the amount it protrudes into
the blood flow, a method is needed for computing the strain. If the temperature profile of
the entire SMA diaphragm is nearly constant, the strain is easily computed based upon a
single temperature. However, the situation is more complex if the temperature profile
along the diaphragm significantly varies.
Figure 6.5 shows the temperature variation along the diaphragm after 0.15
seconds of cooling by conduction only. The temperature variation is small since the
141
thermal conductivity of the Nitinol SMA is considerably larger than the conductivity of
the blood. Hence, a lumped capacitance model is acceptable.
Figure 6.6 shows the same type of temperature variation, but now both conduction
and convection to the blood are included. The inclusion of the convection to the blood
creates a much larger temperature variation along the diaphragm. This is as expected
since the heat flux due to convection is quite large.
Given that the strain of the SMA diaphragm is dependent upon the temperature, a
model is needed to determine the total strain of the diaphragm. In this research, the total
strain of the SMA diaphragm is handled by summing the individual strains of all the
diaphragm elements. The strain of each element is computed based upon the phase
transformation of that element only. Since this is a feasibility study, it is assumed that
this ?average? strain of the entire diaphragm still preserves a circular shape and the
geometric analysis of Chapter 3 applies.
Because the location of the diaphragm nodes move as the strain varies, a
sensitivity analysis of the thermal model and diaphragm strain is needed. Figure 6.7
shows a plot of the temperature of the center node of the diaphragm as a function of time
for three different diaphragm strains, namely, 1%, 3%, and 8%. As in the previous
analysis, the blood flow is again assumed to be a constant 0.5 m/s. Based upon Figure
6.7, the center temperature of the diaphragm is very similar for larger strains and
increasing time. Based upon my assumptions that the SMA moves from grid to grid
during the heating and cooling, the lack of sensitivity for different strains is good.
142
6.7 Multicycle Thermal Analysis of Micropump
Figure 6.8 shows the heating and cooling of the SMA diaphragm for three
seconds. The heating for each cycle is arbitrarily selected to be 0.05 seconds so that the
cooling phase is the obvious limiter in performance. As for previous analyses, it is
assumed that the blood flow is a constant 0.5 m/s at the pump entrance.
From Figure 6.8, it is evident that the total cycle time for both heating and cooling
is approximately 0.75 seconds. In order to show the significance of the blood flow on the
cooling of the diaphragm, the case of only conduction is shown. With conduction only,
the first cycle is about 3 seconds and drops off precipitously after the first cycle. Hence,
the flow of blood is a necessity for the pump operation.
Lastly, Figure 6.9 shows the thermal cycles when the dynamics of the human
heart are incorporated. Since the average flow of the blood from the heart is less than 0.5
m/s, the cooling time is increased. The thermal cycle for both heating and cooling is now
approximately 1.2 seconds.
6.8 Micropump Analysis Based upon Thermal Modeling
Using the geometric analysis presented for the SMA diaphragm in Chapter 3, the
volume of pumped medical fluid with each cycle of operation is readily computed. For
the pump designed in this dissertation, the length of the SMA diaphragm is 2 cm in both
length and depth. Using the typical strain limit of 8% for Nitinol with infinite fatigue
life, the volume pumped for each cycle is
g1848 = g891.5732 g1876 10g2879g2875 (g1865g2871) (6.10)
143
Assuming a fluid with density similar to water, one cubic meter is approximately 1000
liters. Each thermal cycle of the SMA micropump takes approximately 1.2 seconds;
therefore, the pump rate for the pump modeled in this research is listed in Table 6.3.
Pump Range Flow Rate
Minimum Flow
(1 cycle)
0.000g89157g1864g1861g1872g1857g1870g1871g1865g1861g1866
Maximum Flow
(cycle = 1.2 s)
0.04785g1864g1861g1872g1857g1870g1871g1865g1861g1866
Table 6.3 Minimum and Maximum SMA Micropump Capacity for This Design
144
Figure 6.1 Space Discretization for Finite Volumes
i,
j
i+1,j
i,j+1
i-1,j
i,j-1
y?
x?
ix?1?? ix
jy?
1?? jy
145
Figure 6.2 Cardinal Directions for Implementation of Finite Volume Method
N
W E
S
n
s
ew
?x
?y
P
Figure 6.3 Heat Flux from SMA (I,j) to Blood (I,j
146
-1)
147
Figure 6.4 Comparison of Energy Flux between SMA and Blood
0
20000
40000
60000
80000
100000
120000
0 0.5 1 1.5 2 2.5 3 3.5 4
Flu
x (
W
/m
^2
)
Time (s)
Conduction &
Convection
Conduction
only
148
Figure 6.5 Temperature along SMA Diaphragm during Cooling
56.9
56.95
57
57.05
57.1
57.15
57.2
57.25
57.3
57.35
0 0.2 0.4 0.6 0.8 1
Te
mp
erat
ur
e (
C)
Distance along SMA
Conduction only at time = 0.15 sec
149
Figure 6.6 Temperature along SMA Diaphragm during Cooling
40
42
44
46
48
50
52
54
56
58
60
0 0.2 0.4 0.6 0.8 1
Te
mp
erat
ur
e (
C)
Distance along SMA
Conduction & Convection at time = 0.15 sec
150
Figure 6.7 Temperature throughout SMA at Varying Deformation
30
40
50
60
70
80
90
0 0.2 0.4 0.6 0.8
Te
mp
era
tur
e (
C)
Time (s)
Strain (.01)
Strain (.03)
Strain (.08)
151
Figure 6.8 Thermal Cycles with/without Blood Flow
0.00
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0 1 2 3 4
SM
A
Di
ap
hr
agm
St
rai
n
Time (sec)
No blood flow
Blood flow (0.5 m/s)
152
Figure 6.9 Thermal Cycles using Complete Model of Heart Beat
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0 1 2 3 4 5 6 7
SM
A
Di
ap
hr
agm
St
rai
n
Time (sec)
153
Chapter 7 Conclusions
The results of this research show that a micropump operating with a SMA
diaphragm is feasible. Essentially, this micropump is powered by resistive heating and
uses natural flowing blood to cool. The analysis of the micropump shows that human
blood can be used in conjunction with resistive heating to operate the micropump. The
best cycle time with the current design is approximately 1.2 seconds.
Regarding performance, the micropump modeled in this work is not optimized for
performance or efficiency. Both the performance and efficiency may be easily improved
by reducing the mass of the SMA diaphragm. The limit of the SMA mass is governed
by the manufacturing techniques of the SMA diaphragm. A sputtering deposition
technique has been used in the past to create very small SMA membranes. The
membrane thickness modeled here is 0.1 cm, and this may be manufactured to be
considerably thinner. Other options to improve performance are relocating the internal
ramp or using a completely different method of routing the blood.
There are several things that may be done to enhance the numerical modeling in
this research. The transient velocity for a heart beat is probably not constant for different
people. Future study is needed in this area so that the micropump may be designed for
the worst-case scenario. All of the internal obstacles in this micropump (ramp and
diaphragm) are implemented using first order bounce-back technique. All of the exterior
boundaries are second order boundaries; therefore, the internal bounce-back can be
improved. The SMA diaphragm in this work assumes the diaphragm always lies on a
154
grid point. This may be improved by interpolating between grids and using an advanced
boundary technique. Lastly, a 3-D model should be implemented to account for a
circular entrance (artery) into a rectangular cavity. This allows for sophisticated
entrances and exits in the pump chamber.
This research sufficiently shows the feasibility of this micropump, and it is
reasonable to begin experimental work to validate the model. The experimental work
should generate the scalability limits of the micropump. Once a micropump is
constructed, it can be experimentally calibrated for precision and capacity.
155
Bibliography
Abdelgawad, M., I. Hassan, et al. (2005). "Numerical investigation of multistage viscous
micropump configurations." Journal of Fluids Engineering, Transactions of the
ASME 127(4): 734-742.
Abe, T. (1997). "Derivation of the Lattice Boltzmann Method by Means of the Discrete
Ordinate Method for the Boltzmann Equation." Journal of Computational Physics
131(1): 241.
Adnyana, D. N. (1984). "A Copper based alloy for shape memory wires and springs."
Wire Journal International April: 52-61.
Agarwal, S., N. Verma, et al. (2005). "A lattice Boltzmann model for adsorption
breakthrough." Heat and Mass Transfer/Waerme- und Stoffuebertragung 41(9):
843-854.
Ansumali, S. (2004). "Entropic Lattice Boltzmann Simulation of the Flow Past Square
Cylinder." International journal of Modern Physics C 15(3): 435-445.
Ansumali, S. and I. V. Karlin (2000). "Stabilization of the lattice Boltzmann method by
the H theorem: A numerical test." Physical Review E. Statistical Physics,
Plasmas, Fluids, and Related Interdisciplinary Topics 62(6 A): 7999-8003.
Ansumali, S. and I. V. Karlin (2002). "Kinetic boundary conditions in the lattice
Boltzmann method." Physical Review E - Statistical Physics, Plasmas, Fluids, and
Related Interdisciplinary Topics 66(2 2): 026311-1.
Ansumali, S. and I. V. Karlin (2002). "Single relaxation time model for entropic lattice
Boltzmann methods." Physical Review E - Statistical, Nonlinear, and Soft Matter
Physics 65(5): 056312.
Ansumali, S. and I. V. Karlin (2005). "Consistent lattice Boltzmann method." Physical
Review Letters 95(26): 260605.
Ansumali, S., I. V. Karlin, et al. (2006). "Entropic lattice Boltzmann method for
microflows." Physica A: Statistical Mechanics and its Applications 359(1-4): 289-
305.
Ansumali, S., I. V. Karlin, et al. (2003). "Minimal Entropic Kinetic Models for
Hydrodynamics." Europhysics Letters 63(6): 798-804.
Aoki, K., H. Yoshida, et al. (2003). "Inverted velocity profile in the cylindrical Couette
flow of a rarefied gas." Physical Review E 68(1 2): 16302-1.
Arkilic, E. B., M. A. Schmidt, et al. (1997). "Gaseous slip flow in long microchannels."
Journal of Microelectromechanical Systems 6(2): 167-178.
Artoli, A. M., A. G. Hoekstra, et al. (2006). "Optimizing lattice Boltzmann simulations
for unsteady flows." Computers and Fluids 35(2): 227-240.
Bakken, E. E. and K. Heruth (1991). "Temporal control of drugs: an engineering
perspective." Ann NY Academy of Science 618: 422-427.
156
Barison, S., S. Cattarin, et al. (2004). "Characterisation of surface oxidation of nickel-
titanium alloy by ion-beam and electrochemical techniques." Electrochimica Acta
50(1): 11-18.
Benard, W. L., H. Kahn, et al. (1997). Titanium-nickel shape-memory alloy actuated
micropump. Proceedings of the 1997 International Conference on Solid-State
Sensors and Actuators. Part 1 (of 2), Jun 16-19 1997, Chicago, IL, USA, IEEE,
Piscataway, NJ, USA.
Benard, W. L., H. Kahn, et al. (1998). "Thin-film shape-memory alloy actuated
micropumps." Journal of Microelectromechanical Systems 7(2): 245-251.
Bernsdorf, J. (2001). "Introduction to Lattice Gas and Lattice Boltzmann Methods."
Workshop LBM.
Boghosian, B. (2003). "Entropic lattice Boltzmann models for computational fluid
dynamics." Presentation.
Boghosian, B. and J. Boon (2005). "Lattice Boltzmann and nonextensive diffusion."
Europhysics News(November/December): 192-194.
Boghosian, B., P. J. Love, et al. (2004). "Entropic lattice Boltzmann model for Burger's
equation." Phil. Trans. R. Soc. Lond. 362: 1691-1701.
Boghosian, B., J. Yepez, et al. (2001). "Entropic Lattice Boltzmann Methods." Work in
Progress.
Boghosian, B. M., P. J. Love, et al. (2003). "Galilean-invariant lattice-Boltzmann models
with H theorem." Physical Review E 68(2 2): 025103-1.
Boghosian, B. M., P. J. Love, et al. (2004). "Galilean-invariant multi-speed entropic
lattice Boltzmann models." Physica D: Nonlinear Phenomena 193(1-4): 169-181.
Bouzidi, M., M. Firdaouss, et al. (2001). "Momentum transfer of a Boltzmann-lattice
fluid with boundaries." Physics of Fluids 13(11): 3452-3459.
Buick, J. and C. Greated (2000). "Gravity in a lattice boltzmann model." Physical Review
E 61(5): 5307-5320.
Chen, H., B. Boghosian, et al. (2003). "A Ternary Lattice Boltzmann Model for
Amphiphilic Fluids."
Chen, H. and C. Teixeira (2000). "H-theorem and origins of instability in thermal lattice
Boltzmann models." Computer Physics Communications 129: 21-31.
Chen, M. F., X. J. Yang, et al. (2004). "Bioactive NiTi shape memory alloy used as bone
bonding implants." Materials Science and Engineering: C 24(4): 497-502.
Chen, S., D. Martinez, et al. (1996). "On boundary conditions in lattice Boltzmann
methods." Physics of Fluids L2 - http://dx.doi.org/10.1063/1.869035 8(9): 2527.
Chikatamarla, S. S., N. I. Prasianakis, et al. (2005). "Entropic Lattice Boltzmann Method
for Simulation of Thermal Flows." Work in Progress.
Chu, P. K. (2006). "Bioactivity of plasma implanted biomaterials." Nuclear Instruments
and Methods in Physics Research Section B: Beam Interactions with Materials
and Atoms 242(1-2): 1-7.
Curry, D. T. (1979). "New uses for metals that remember." Machine Design(October):
113-117.
Date, H. and M. Shimozuma (2001). "Boltzmann equation description of electron
transport in an electric field with cylindrical or spherical symmetry." Physical
Review E. Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary
Topics 64(6 II): 066410-1.
157
De La Fuente, L. F., D. M. Causon, et al. (2003). "Towards simulating urban canyon
circulations with a 2D lattice Boltzmann model." Environmental Modelling and
Software 18(1): 71-79.
Desplat, J.-C., I. Pagonabarraga, et al. (2001). "Ludwig: A parallel Lattice-Boltzmann
code for complex fluids." Computer Physics Communications 134(3): 273-290.
Do-Quang, M., E. Aurell, et al. (2000). "An inventory of lattice Boltzmann models of
multiphase flows." Project 24082-61206 Parallel Scientific Computing Institute:
1-51.
Doolen, G. (1987). Lattice Gas Methods for Partial Differential Equations, Addison-
Wesley.
Duerig, T., A. Pelton, et al. (1999). "Overview of nitinol medical applications."
Proceedings of the 1998 International Conference on Martensitic Transformations
(ICOMAT 98), Dec 7-Dec 11 1998
Materials Science and Engineering A: Structural Materials: Properties, Microstructure
and Processing A273-275: 149-160.
Dupin, M. M., T. J. Spencer, et al. (2004). "A Many-Component Lattice Boltzmann
Equation Simulation for Transport of Deformable Particles." Phil. Trans. R. Soc.
Lond.(362): 1-31.
Erickson, D. (2005). "Towards numerical prototyping of labs-on-chip: Modeling for
integrated microfluidic devices." Microfluidics and Nanofluidics 1(4): 301-318.
Erickson, D. and D. Li (2003). "Analysis of alternating current electroosmotic flows in a
rectangular microchannel." Langmuir 19(13): 5421-5430.
Erickson, D. and D. Li (2004). "Integrated microfluidic devices." Analytica Chimica Acta
507(1): 11-26.
Es-Souni, M., M. Es-Souni, et al. (2005). "Assessing the biocompatibility of NiTi shape
memory alloys used for medical applications." Analytical and Bioanalytical
Chemistry 381(3): 557-567.
Fang, H.-P., R.-Z. Wan, et al. (2002). "Lattice Boltzmann model with nearly constant
density." Physical Review E - Statistical Physics, Plasmas, Fluids, and Related
Interdisciplinary Topics 66(3 2B): 036314-1.
Fang, H., Z. Lin, et al. (1998). "Lattice Boltzmann simulation of viscous fluid systems
with elastic boundaries." Physical Review E. Statistical Physics, Plasmas, Fluids,
and Related Interdisciplinary Topics 57(1): 25.
Fang, H., Z. Wang, et al. (2002). "Lattice Boltzmann method for simulating the viscous
flow in large distensible blood vessels." Physical Review E - Statistical Physics,
Plasmas, Fluids, and Related Interdisciplinary Topics 65(5 1): 51925-1.
Firstov, G. S., R. G. Vitchev, et al. (2002). "Surface oxidation of NiTi shape memory
alloy." Biomaterials 23(24): 4863-4871.
Ford, D. S. and S. R. White (1996). "Thermomechanical behavior of 55Ni45Ti nitinol."
Acta Materialia 44(6): 2295-2307.
Funakubo, H. (1987). Shape Memory Alloys, Gordon and Breach Science Publishers.
Gil, F. J. and J. A. Planell (1998). "Shape Memory Alloys for Medical Applications."
Proceedings of the Institution of Mechanical Engineers -- Part H -- Journal of
Engineering in Medicine 212(6): 473-488.
Ginzburg, I. and D. D'Humieres (2003). "Multireflection boundary conditions for lattice
Boltzmann models." Physical Review E 68(6 2): 66614-1.
158
Ginzburg, I. and D. D.d'Humieres (1996). "Local second-order boundary methods for
lattice Boltzmann models." Journal of Statistical Physics 84(5): 927-971.
Golestaneh, A. A. (1984). "Shape memory phenomena." Physics Today(April): 62-70.
Gorban, A. and I. V. Karlin (2005). "Invariance correction to Grad's equations: where to
go beyond approximations." Continuum Mech. Thermodyn. 17(4): 311-335.
Grayson, A. C. R., R. S. Shawgo, et al. (2004). "A BioMEMS review: MEMS technology
for physiologically integrated devices." Biomedical Applications for Mems and
Microfluidics
Proceedings of the IEEE 92(1): 6-21.
Grayson, A. C. R., R. S. Shawgo, et al. (2004). "Electronic MEMS for triggered
delivery." Advanced Drug Delivery Reviews 56(2): 173-184.
Guo, Z.-L., C.-G. Zheng, et al. (2002). "Non-equilibrium extrapolation method for
velocity and pressure boundary conditions in the lattice Boltzmann method."
Chinese Physics 11(4): 366-374.
Guo, Z. and T. S. Zhao (2002). "Lattice Boltzmann model for incompressible flows
through porous media." Physical Review E - Statistical Physics, Plasmas, Fluids,
and Related Interdisciplinary Topics 66(3 2B): 036304-1.
Guo, Z., C. Zheng, et al. (2002). "Discrete lattice effects on the forcing term in the lattice
Boltzmann method." Physical Review E. Statistical Physics, Plasmas, Fluids, and
Related Interdisciplinary Topics 65(4 2B): 46308-1.
Guo, Z., C. Zheng, et al. (2002). "An extrapolation method for boundary conditions in
lattice Boltzmann method." Physics of Fluids 14(6): 2007-2010.
Halliday, I., L. A. Hammond, et al. (2001). "Lattice Boltzmann equation
hydrodynamics." Physical Review E. Statistical Physics, Plasmas, Fluids, and
Related Interdisciplinary Topics 64(1 I): 011208-1.
He, Q., W. M. Huang, et al. (2004). "Characterization of sputtering deposited NiTi shape
memory thin films using a temperature controllable atomic force microscope."
Smart Materials and Structures 13: 977-982.
He, X. and L.-S. Luo (1997). "Lattice Boltzmann Model for the Incompressible Navier-
Stokes Equation." Journal of Statistical Physics 88(3-4): 927.
He, X. and L.-S. Luo (1997). "Priori derivation of the lattice Boltzmann equation."
Physical Review E. Statistical Physics, Plasmas, Fluids, and Related
Interdisciplinary Topics 55(6-A): 6333.
He, X. and L.-S. Luo (1997). "Theory of the lattice Boltzmann method: from the
Boltzmann equation to the lattice Boltzmann equation." Physical Review E.
Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 56(6):
6811.
He, X. and Z. Q. (1995). "Analysis and boundary condition of the lattice Boltzmann BGK
model with two velocity components." Rough draft.
Hinton, F. L., M. N. Rosenbluth, et al. (2001). "Modified Lattice Boltzmann method for
compressible fluid simulations." Physical Review E. Statistical Physics, Plasmas,
Fluids, and Related Interdisciplinary Topics 63(6 I): 061212-1.
Hirabayashi, M., M. Ohta, et al. (2004). "A lattice Boltzmann study of blood flow in
stented aneurism." Future Generation Computer Systems 20(6 SPEC ISS): 925-
934.
159
Hodgson, D. (2005). "Nitinol in Medical Devices." Advanced Materials & Processes
163(1): 63-63.
Holdych, D. J., J. G. Georgiadis, et al. (2001). "Migration of a van der Waals bubble:
Lattice Boltzmann formulation." Physics of Fluids 13(4): 817-825.
Hou, S., Q. Zou, et al. (1995). "Simulation of Cavity Flow by the Lattice Boltzmann
Method." Journal of Computational Physics 118(2): 329.
Hu, T., C.-l. Chu, et al. (2007). "In vitro biocompatibility of titanium-nickel alloy with
titanium oxide film by H2O2 oxidation." Transactions of Nonferrous Metals
Society of China 17(3): 553-557.
Ikuta, K. (1990). "Micro/miniature shape memory alloy actuator." Proceedings of the
1990 IEEE International Conference on Robotics and Automation: 2156-2161.
Ikuta, K., A. Takahashi, et al. (2003). Fully integrated micro biochemical laboratory
using biochemical IC chips - Cell-free protein synthesis by using a built-in
micropump chip. IEEE Sixteenth Annual International Conference on Micro
Electro Mechanical Systems, Jan 19-23 2003, Kyoto, Japan, Institute of Electrical
and Electronics Engineers Inc.
Inamuro, T., M. Yoshino, et al. (1995). "Non-slip boundary condition for lattice
Boltzmann simulations." Physics of Fluids 7(12): 2928.
Inamuro, T., M. Yoshino, et al. (1999). "Lattice Boltzmann simulation of flows in a
three-dimensional porous structure." International Journal for Numerical Methods
in Fluids
Proceedings of the 1997 10th International Congress on Nmerical Methods in Laminar
and Turbulent Flow 29(7): 737-748.
Jackson, C. M., H. J. Wagner, et al. (1972). "Nitinol-the alloy with a memory: Its
physical metallurgy, properties, and applications." NASA-SP 5110.
Kahn, H., M. A. Huff, et al. (1998). "TiNi shape-memory alloy and its applications for
MEMS." Journal of Micromechanics and Microengineering 8(3): 213-221.
Kan, J., Z. Yang, et al. (2005). "Design and test of a high-performance piezoelectric
micropump for drug delivery." Sensors and Actuators, A: Physical 121(1): 156-
161.
Kandhai, D., A. Koponen, et al. (2001). "Iterative momentum relaxation for fast lattice-
Boltzmann simulations." Future Generation Computer Systems 18(1): 89-96.
Kandhai, D., D. J.-E. Vidal, et al. (1999). "Lattice-Boltzmann and finite element
simulations of fluid flow in a SMRX static mixer reactor." International Journal
for Numerical Methods in Fluids 31(6): 1019-1033.
Kang, J. (2005). "The lattice gas model and lattice boltzmann model on hexagonal grids."
Thesis Auburn University.
Karlin, I. V., S. Ansumali, et al. "Entropic lattice Boltzmann method for large scale
turbulence simulation." Not published yet!
Kataoka, T. and M. Tsutahara (2004). "Lattice Boltzmann model for the compressible
Navier-Stokes equations with flexible specific-heat ratio." Physical Review E
69(3 2): 035701-1.
Kujala, S., A. Pajala, et al. (2004). "Biocompatibility and strength properties of nitinol
shape memory alloy suture in rabbit tendon." Biomaterials 25(2): 353-358.
Lallemand, P. and L.-S. Luo (2000). "Theory of the lattice Boltzmann method:
dispersion, dissipation, isotropy, Galilean invariance, and stability." Physical
160
Review E. Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary
Topics 61(6 A): 6546-6562.
Laser, D. J. and J. G. Santiago (2004). "A review of micropumps." Journal of
Micromechanics and Microengineering 14(6): 35-64.
Lee, T. and C.-L. Lin (2005). "Rarefaction and compressibility effects of the lattice-
Boltzmann-equation method in a gas microchannel." Physical Review E -
Statistical, Nonlinear, and Soft Matter Physics 71(4): 046706.
Li, B. and D. Y. Kwok (2003). "Discrete Boltzmann equation for microfluids." Physical
Review Letters 90(12): 124502-1.
Li, B. and D. Y. Kwok (2003). "Lattice Boltzmann model of microfluidics in the
presence of external forces." Journal of Colloid and Interface Science 263(1):
144-151.
Li, B. and D. Y. Kwok (2004). "A Lattice Boltzmann model with high Reynolds number
in the presence of external forces to describe microfluidics." Heat and Mass
Transfer/Waerme- und Stoffuebertragung 40(11): 843-851.
Li, H., X. Lu, et al. (2004). "Force evaluations in lattice Boltzmann simulations with
moving boundaries in two dimensions." Physical Review E 70(2 2): 026701-1.
Li, Y., J. LeBoeuf, et al. (2004). "Least-squares finite-element lattice Boltzmann
method." The American Physical Society 69.
Li, Y., R. Shock, et al. (2004). "Numerical study of flow past an impulsively started
cylinder by the lattice-Boltzmann method." Journal of Fluid Mechanics 519: 273-
300.
Liang, C. (1990). The Constitutive Modeling of Shape Memory Alloys, PhD
Disseratation, Virginia Polytechnic Institute and State University.
Liang, C. and C. A. Rogers (1990). "One-dimensional thermomechanical constitutive
relations for shape memory materials." Journal of Intelligent Material Systems
and Structures 1(April): 207-234.
Lim, C. Y., C. Shu, et al. (2002). "Application of lattice Boltzmann method to simulate
microchannel flows." Physics of Fluids 14(7): 2299-2308.
Lin, Z., H. Fang, et al. (1996). "Improved lattice Boltzmann model for incompressible
two-dimensional steady flows." Physical Review E. Statistical Physics, Plasmas,
Fluids, and Related Interdisciplinary Topics 54(6): 6323 L2 -
http://dx.doi.org/10.1103/PhysRevE.54.6323.
Love, P. J. and B. M. Boghosian (2004). "On the dependence of the Navier-Stokes
equations on the distribution of molecular velocities." Physica D: Nonlinear
Phenomena 193(1-4): 182-194.
Luo, L.-S. (2003). "Lattice Boltzmann Methods for Computational Fluid Dynamics."
Lecture Notes.
Luo, L.-S., D. Qi, et al. (2002). "Applications of the Lattice Boltzmann Method to
Complex and Turbulent Flows." Report.
Maier, R. S., R. S. Bernard, et al. (1996). "Boundary conditions for the lattice Boltzmann
method." Physics of Fluids 8(7): 1788 L2 - http://dx.doi.org/10.1063/1.868961.
Maillefer, D., S. Gamper, et al. (2001). A high-performance silicon micropump for
disposable drug delivery systems. 14th IEEE International Conference on Micro
Electro Mechanical Systems (MEMS 2001), Jan 21-25 2001, Interlaken, Institute
of Electrical and Electronics Engineers Inc.
161
Majdalani, J. and H. A. Chibli (2002). "Pulsatory channel flows with arbitrary pressure
gradients." 3rd AIAA Theoretical Fluid Mechanics Meeting June: 1-13.
Makino, E., T. Mitsuya, et al. (2000). "Dynamic actuation properties of TiNi shape
memory diaphragm." Sensors and Actuators, A: Physical 79(2): 128-135.
Makino, E., T. Mitsuya, et al. (2000). "Micromachining of TiNi shape memory thin film
for fabrication of micropump." Sensors and Actuators, A: Physical 79(3): 251-
259.
Mantovani, D. (2000). "Shape memory alloys: Properties and biomedical applications."
JOM Journal of the Minerals, Metals and Materials Society 52(10): 36-44.
Martys, N. S. (2001). "Improved approximation of the Brinkman equation using a lattice
Boltzmann method." Physics of Fluids 13(6): 1807-1810.
Martys, N. S., X. Shan, et al. (1998). "Evaluation of the external force term in the discrete
Boltzmann equation." Physical Review E. Statistical Physics, Plasmas, Fluids,
and Related Interdisciplinary Topics 58(5-B): 6855.
Mei, R., L.-S. Luo, et al. (2000). "An accurate curved boundary treatment in the lattice
Boltzmann method." NASA/CR 2000-209854.
Melchionna, S. and S. Succi (2004). "Electrorheology in nanopores via lattice Boltzmann
simulation." Journal of Chemical Physics 120(9): 4492-4497.
Michiardi, A., C. Aparicio, et al. (2007). "Electrochemical behaviour of oxidized NiTi
shape memory alloys for biomedical applications." Surface and Coatings
Technology 201(14): 6484-6488.
Miles, H. (1989). "Shape memory materials." Engineering July/August: 29-30.
Miller, W. (1994). "Flow in the Driven Cavity Calculated by the Lattice Boltzmann
Method." Phys. Rev. E.
Nathanson, H. C. (1967). "The Resonant Gate Transistor." IEEE Transaction Electron
Devices 14(3): 117-133.
Nekovee, M., P. V. Coveney, et al. (2000). "Lattice-Boltzmann model for interacting
amphiphilic fluids." Physical Review E. Statistical Physics, Plasmas, Fluids, and
Related Interdisciplinary Topics 62(6 B): 8282-8294.
Nguyen, N.-T., X. Huang, et al. (2002). "MEMS-micropumps: A review." Journal of
Fluids Engineering, Transactions of the ASME 124(2): 384-392.
Nie, X., G. Doolen, et al. (2002). "Lattice-Boltzmann simulations of fluid flow in
MEMS." Journal of Statistical Physics 107(112): 279-289.
Nie, X., Y.-H. Qian, et al. (1998). "Lattice Boltzmann simulation of the two-dimensional
Rayleigh-Taylor instability." Physical Review E. Statistical Physics, Plasmas,
Fluids, and Related Interdisciplinary Topics 58(5-B): 6861.
Noble, D. R., S. Chen, et al. (1995). "Consistent hydrodynamic boundary condition for
the lattice Boltzmann method." Physics of Fluids 7(1): 203.
Onishi, J., Y. Chen, et al. (2001). "Lattice Boltzmann simulation of natural convection in
a square cavity." JSME International Journal, Series B: Fluids and Thermal
Engineering 44(1): 53-62.
Pan, C., J. F. Prins, et al. (2004). "A high-performance lattice Boltzmann implementation
to model flow in porous media." Computer Physics Communications 158(2): 89-
105.
Pelton, A. R., T. W. Duerig, et al. (2005). "Nitinol medical devices." Advanced Materials
and Processes 163(10): 63-65.
162
Perkins, J. (1975). Shape Memory Effects in Alloys. New York, Plenum.
Petrini, L., F. Migliavacca, et al. (2005). "Computational Studies of Shape Memory Alloy
Behavior in Biomedical Applications." Journal of Biomechanical Engineering
127(4): 716-725.
Pohl, T., F. Deserno, et al. (2004). Performance evaluation of parallel large-scale lattice
boltzmann applications on three supercomputing architectures. IEEE/ACM
SC2004 Conference - Bridging Communities, Nov 6-12 2004, Pittsburgh, PA,
United States, Institute of Electrical and Electronics Engineers Inc., New York,
NY 10016-5997, United States.
Porter, D. A. and K. E. Easterling (1981). Phase Transformation in Metals and Alloys,
Chapman and Hall Publishing.
Qian, Y.-H. and Y.-F. Deng (1997). "Lattice BGK model for viscoelastic media."
Physical Review Letters 79(14): 2742.
Qian, Y.-H. and Y. Zhou (2000). "Higher-order dynamics in lattice-based models using
the Chapman-Enskog method." Physical Review E. Statistical Physics, Plasmas,
Fluids, and Related Interdisciplinary Topics 61(2): 2103.
Qian, Y., D. D'Humieres, et al. (1992). "Lattice BGK Models for Navier-Stokes
Equation." Europhysics Letters 17(6): 479-484.
Quadir, G. A., A. Mydin, et al. (2001). "Analysis of microchannel heat exchangers using
FEM." International Journal of Numerical Methods for Heat and Fluid Flow
11(1): 59-75.
Raabe, D. (2004). "Overview of the lattice Boltzmann method for nano- And microscale
fluid dynamics in materials science and engineering." Modelling and Simulation
in Materials Science and Engineering 12(6): 13-46.
Rasmussen, A. and M. E. Zaghloul (1999). Pumping techniques available for use in
biomedical analysis systems. 1999 IEEE 42nd Midwest Symposium on Circuits
and Sistems, Aug 08-Aug 11 1999, Las Cruces, NM, USA, Institute of Electrical
and Electronics Engineers Inc., Piscataway, NJ, USA.
Rector, D., J. Cuta, et al. (1998). "Lattice-Boltzmann simulation code development for
micro-fluidic systems."
Reider, M. B. and J. D. Sterling (1995). "Accuracy of discrete-velocity BGK models for
the simulation of the incompressible Navier-Stokes equations." Computers &
Fluids 24(4): 459-467.
Richter, M., R. Linnemann, et al. (1998). "Robust design of gas and liquid micropumps."
Sensors and Actuators A: Physical 68(1-3): 480-486.
Rogers, C. A., C. Liang, et al. (1991). "Modeling of shape memory alloy hybrid
composites for structural acoustic control." Journal of Acoustical Society of
America 89(January): 210-220.
Rohde, M., D. Kandhai, et al. (2003). "Improved bounce-back methods for no-slip walls
in lattice-Boltzmann schemes: Theory and simulations." Physical Review E 67(6
2): 066703-1.
Rohlf, K. and G. Tenti (2001). "Role of the Womersley number in pulsatile blood flow -
a theoretical study of the Casson model." Journal of Biomechanics 34(1): 141-
148.
163
Ryh?nen, J., M. Kallioinen, et al. (1999). "Bone modeling and cell-material interface
responses induced by nickel-titanium shape memory alloy after periosteal
implantation." Biomaterials 20(14): 1309-1317.
Sankaranarayanan, K., X. Shan, et al. (1999). "Bubble flow simulations with the lattice
Boltzmann method." Chemical Engineering Science
Proceedings of the 1999 4th International Conference on Gas-Liquid and Gas-Liquid-
Solid Reactor Engineering, Aug 23-Aug 25 1999 54(21): 4817-4823.
Satofuka, N. and T. Nishioka (1999). "Parallelization of lattice Boltzmann method for
incompressible flow computations." Computational Mechanics 23(2): 164-171.
Schetky, L. M. (1981). "Shape memory alloys." Scientific American 241(5): 74-82.
Sert, C. and A. Beskok (2003). "Numerical simulation of reciprocating flow forced
convection in two-dimensional channels." Journal of Heat Transfer 125(3): 403-
412.
Shahin, A. R., P. Meckl, et al. (1994). "Enhanced cooling of shape memory alloy wires
using semiconductor heat pump modules." Journal of Intelligent Material Systems
and Structures 5(1): 95-104.
Shawgo, R. S., A. C. R. Grayson, et al. (2002). "BioMEMS for drug delivery." Current
Opinion in Solid State and Materials Science 6(4): 329-334.
Shi, X., J. Lin, et al. (2003). "Discontinuous Galerkin spectral element lattice Boltzmann
method on triangular element." International Journal for Numerical Methods in
Fluids 42(11): 1249-1261.
Shin, D. D., K. P. Mohanchandra, et al. (2005). "Development of hydraulic linear
actuator using thin film SMA." Sensors and Actuators, A: Physical 119(1): 151-
156.
Shishkovsky, I. V. (2008). "Porous Biocompatible Implants and Tissue Scaffolds
Synthesized by Selective Laser Sintering from Ti and NiTi." Journal of Materials
Chemistry 18: 1309-1317.
Shoji, S. and M. Esashi (1994). "Microflow devices and systems." Journal of
Micromechanics and Microengineering 4(4): 157-171.
Shu, C., X. D. Niu, et al. (2005). "Taylor series expansion- and least square-based Lattice
Boltzmann method: An efficient approach for simulation of incompressible
viscous flows." Progress in Computational Fluid Dynamics 5(1-2): 27-36.
Sinton, D., D. Erickson, et al. (2002). "Photo-injection based sample design and
electroosmotic transport in microchannels." Journal of Micromechanics and
Microengineering 12(6): 898-904.
Sone, Y. and T. Doi (2000). "Analytical study of bifurcation of a flow of a gas between
coaxial circular cylinders with evaporation and condensation." Physics of Fluids
12(10): 2639-2660.
Stepan, L. L., D. S. Levi, et al. (2005). "A thin film nitinol heart valve." Journal of
Biomechanical Engineering 127(6): 915-918.
Sterling, J. D. and S. Chen (1996). "Stability Analysis of Lattice Boltzmann Methods."
Journal of Computational Physics: 196.
Stoeckel, D. (1989). "Fabrication and properties of Nickel-Titanium shape memory alloy
wires." Wire Journal International April: 30-40.
Sun, C. and A. Hsu (2004). "Multi-level lattice Boltzmann model on square lattice for
compressible flows." Computers and Fluids 33(10): 1363-1385.
164
Sun, C. and A. T. Hsu (2003). "Three-dimensional lattice Boltzmann model for
compressible flows." Physical Review E 68(1 2): 16303-1.
Swift, M. R., E. Orlandini, et al. (1996). "Lattice Boltzmann simulations of liquid-gas and
binary fluid systems." Physical Review E. Statistical Physics, Plasmas, Fluids,
and Related Interdisciplinary Topics L2 -
http://dx.doi.org/10.1103/PhysRevE.54.5041 54(5): 5041.
Swift, M. R., W. R. Osborn, et al. (1995). "Lattice Boltzmann simulation of nonideal
fluids." Physical Review Letters 75(5): 830.
Takada, N. and M. Tsutahara (1998). "Evolution of viscous flow around a suddenly
rotating circular cylinder in the lattice Boltzmann method." Computers & Fluids
27(7): 807-828.
Tanaka, K. (1986). "A thermomechanical sketch of shape memory effect: One-
dimensional tensile behavior." Res. Mechanica: 251-263.
Tanaka, K. (1990). "A phenomenological description on thermomechanical behavior of
shape memory alloys." Journal of Pressure Vessel Technology 112(May): 158-
163.
Tang, G. H., W. Q. Tao, et al. (2005). "Lattice Boltzmann method for gaseous
microflows using kinetic theory boundary conditions." Physics of Fluids 17(5):
058101.
Tang, G. H., W. Q. Tao, et al. (2005). "Three-dimensional lattice Boltzmann model for
gaseous flow in rectangular microducts and microscale porous media." Journal of
Applied Physics 97(10): 104918.
Thorne, D. and M. Sukop (2004). "Lattice Boltzmann Method for the Elder Problem."
CMWR 2004, FIU.
Thrasher, M. A. (1991). A Fundamental Characterization of Shape Memory Alloy
Actuators for Use in Hand Protheses. Mechanical Engineering. W. Lafayette,
Purdue University: 170.
Thrasher, M. A., A. R. Shahin, et al. (1992). "Thermal cycling of shape memory alloy
wires using semiconductor heat pump modules." Proceedings of the First
European Conference on Smart Structures & Materials.
Tuominen, S. M. and R. J. Biermann (1988). "Shape memory wires." Journal of Metals
February: 32-35.
USNitinol. (1989). "Summary of physical and mechanical properties of nominal 55-
Nitinol."
Verberg, R. and A. J. C. Ladd (2002). "Accuracy and stability of a lattice-Boltzmann
model with subgrid scale boundary conditions." Physical Review E. Statistical
Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 65(1 II): 016701-1.
Verberg, R., J. M. Yeomans, et al. (2005). "Modeling the flow of fluid/particle mixtures
in microchannels: Encapsulating nanoparticles within monodisperse droplets."
Journal of Chemical Physics 123(22): 224706.
Vitiello, A., G. Giorleo, et al. (2005). "Analysis of thermomechanical behaviour of
Nitinol wires with high strain rates." Smart Materials and Structures 14(1): 215-
221.
Voskerician, G., M. S. Shive, et al. (2003). "Biocompatibility and biofouling of MEMS
drug delivery devices." Biomaterials 24(11): 1959-1967.
165
Wang, D. and J. Bernsdorf (2009). "Lattice Boltzmann simulation of steady non-
Newtonian blood flow in a 3D generic stenosis case." Computers & Mathematics
with Applications 58(5): 1030-1034.
Wang, X.-X. (1998). "Low temperature oxidation of TiNi shape memory alloy and its
effect on Ni ion release in saline solution." Journal of Materials Science Letters
17(5): 375-376.
Watanabe, T. and K. Ebihara (2002). "Simulation of two-phase flows and numerical
evaluation of interfacial area." JSME International Journal, Series B: Fluids and
Thermal Engineering 45(3): 600-606.
Wayman, C. M. (1980). "Some applications of shape memory alloys." Journal of Metals
June: 129-137.
Wu, J.-S. and Y.-L. Shao (2004). "Simulation of lid-driven cavity flows by parallel lattice
Boltzmann method using multi-relaxation-time scheme." International Journal for
Numerical Methods in Fluids 46(9): 921-937.
Xuan, Y. and Z. Yao (2005). "Lattice boltzmann model for nanofluids." Heat and Mass
Transfer/Waerme- und Stoffuebertragung 41(3): 199-205.
Yeung, K. W. K., R. Y. L. Chan, et al. (2007). "In vitro and in vivo characterization of
novel plasma treated nickel titanium shape memory alloy for orthopedic
implantation." Surface and Coatings Technology 202(4-7): 1247-1251.
Yoshida, H. (2005). "The wide variety of possible applications of micro-thermofluid
control." Microfluidics and Nanofluidics 1(4): 289-300.
Yoshino, M. and T. Inamuro (2003). "Lattice Boltzmann simulations for flow and
heat/mass transfer problems in a three-dimensional porous structure."
International Journal for Numerical Methods in Fluids 43(2): 183-198.
Yu, D., R. Mei, et al. (2003). "Viscous flow computations with the method of lattice
Boltzmann equation." Progress in Aerospace Sciences 39(5): 329-367.
Yun, L. and et al. (2003). "Thermal responses of shape memory alloy artificial anal
sphincters." Smart Materials and Structures 12(4): 533.
Zhang, R., X. He, et al. (2000). "Interface and surface tension in incompressible lattice
Boltzmann multiphase model." Computer Physics Communications 129(1): 121-
130.
Zhang, R., X. He, et al. (2001). "Surface tension effects on two-dimensional two-phase
Kelvin-Helmholtz instabilities." Advances in Water Resources 24(3-4): 461-478.
Zhang, Y., R. Qin, et al. (2005). "Lattice Boltzmann simulation of rarefied gas flows in
microchannels." Physical Review E - Statistical, Nonlinear, and Soft Matter
Physics 71(4): 047702.
Zhou, Y., R. Zhang, et al. (2004). "Numerical simulation of laminar and turbulent
buoyancy-driven flows using a lattice Boltzmann based algorithm." International
Journal of Heat and Mass Transfer 47(22): 4869-4879.
Zou, Q. and X. He (1997). "On pressure and velocity boundary conditions for the lattice
Boltzmann BGK model." Physics of Fluids 9(6): 1591.