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.