An Accessible Platform for Numerical and Experimental Investigations of Taylor Dispersion by Abigail W. Taylor Sc.M., Brown University Thesis submitted in partial fulfillment of the requirements for the Degree of Master of Science in the School of Engineering at Brown University PROVIDENCE, RHODE ISLAND May 2019 Signature Page This thesis by Abigail W. Taylor is accepted in its present form by the School of Engineering as satisfying the thesis requirements for the degree of Master of Science. Date: Daniel M. Harris, Ph.D., Advisor Approved by the Graduate Council Date: Andrew G. Campbell, Dean of the Graduate School ii Abstract Within the field of microfluidics, scientists have discovered methods to manipulate flow at the micron scale, downsizing standard laboratory procedures. These techniques boast faster results while consuming less resources. Microfluidic devices have been used for diagnostic testing, drug delivery, chemical prepa- ration, and more. Of interest for many of these applications is the behavior of solute spreading in fully developed laminar flow, first examined by Sir Geoffrey Taylor.1 An initially localized solute under these flow conditions will exhibit enhanced spreading from the combined effects of fluid transport and molecu- lar diffusion versus purely molecular diffusion. Literature in this field documents the influence of channel geometry and aspect ratio on the magnitude and time-dependent progression of the solute’s effective disper- sion.2, 3 Historically, numerical and experimental explorations into this so called “Taylor Dispersion” problem have been exclusive to highly trained researchers. This thesis serves to improve accessibility to these studies by providing instructions for optimized numerical simulations using a user-friendly, commercially available multiphysics software, and a refined microchannel fabrication technique utilizing an inexpensive desktop craft cutter. These methods have the potential for significant immediate impact on technology. To validate the accuracy of each study, experiments are modeled after classical Taylor Dispersion studies and bench- marked with known analytical results. Two-dimensional, axisymmetric, and three-dimensional simulations are executed within COMSOL Mul- tiphysics to cover a range of channel geometries and aspect ratios; each agreeing with theoretical predictions for the solute’s distribution. The enclosed simulation results support the capabilities of COMSOL as an ac- curate and streamlined tool for numerically solving solute dispersion problems within fully laminar flows, even for users with no prior simulation experience. iii Acknowledgments First and foremost, I would like to express enormous gratitude to my thesis advisor, Dan Harris, whose guidance, mentorship, immense knowledge, and unwavering support has undoubtedly been the pillar of my success and happiness at Brown University. His encouragement allowed me to reach outside my comfort zone and achieve goals I never thought possible. Concurrently, his moral support enabled me to persevere through the occasionally daunting world of graduate school. The time and energy that Dan dedicates to every advisee, student, and project is truly above and beyond. I could not have asked for a better experience at Brown, and I owe that to Dan. As the first graduate student member of the Harris Lab, I’ve had the phenomenal experience of watching it evolve into a community of smart, collaborative, genuine, and passionate individuals. I am grateful to each of them for their continuous support, cheer, and camaraderie. I wish to acknowledge the incredibly dedicated efforts of undergraduate researchers Jesse Remeis, for her help getting the three-dimensional simulations off the ground; Emma Abele, for her assistance with kickstarting the experiments; and Ian Ho, for being the most reliable engineering whiz and building a syringe pump from scratch for our experimental trials. I would also like to thank Hibbitt Engineering Fellow Jeong- Hyun Kim for first familiarizing me with microfluidic experiments and concepts when I arrived in the lab. I am deeply grateful to Kenny Breuer, who graciously allowed us the use of his inverted microscope and lab space to carry out our experimental trials. The scope of this current thesis would not have been possible without his generosity. Finally, I wish to thank my family for their infinite love and support. They have been by my side at every step of the way, celebrating my successes and holding me up through any difficult times. It is because of them that I have the opportunity to attend graduate school and complete this thesis, for that I could not be more grateful. They will always be my main source of motivation. iv Contents 1 Introduction 1 1.1 Introduction to Taylor Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Microfluidics Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.3 Discovery of an Effective Dispersion Coefficient . . . . . . . . . . . . . . . . . . . . . . . 2 1.4 Governing Equations and Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Benchmarking Theory 6 2.1 Taylor-Aris Dispersivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Dispersion Factor for Rectangular Channels . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Longitudinal Asymmetries in Transitional Time Regimes . . . . . . . . . . . . . . . . . . . 8 3 Optimized Numerical Simulations for Microfluidic Dispersion 12 3.1 Two-Dimensional Channel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.1.1 Problem Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.1.2 COMSOL Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.1.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2 Axisymmetric Channel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2.1 Problem Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2.2 COMSOL Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3 Three-Dimensional Channel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3.1 Problem Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3.2 COMSOL Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4 Optimized Craft-Cutter Fabrication Technique 33 4.1 Introduction to Microchannel Fabrication . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.2 Optimization of Desktop Cutter Technique . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.2.1 Material and Cutter Setting Optimization . . . . . . . . . . . . . . . . . . . . . . . 36 4.2.2 Fabrication Instructions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.3 Experimental Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.3.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5 Discussion 50 A Appendix 53 A.1 Two-Dimensional Velocity Conversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 A.2 Axisymmetric Velocity Conversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 A.3 Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 A.4 Two-Dimensional Simulation Processing Code . . . . . . . . . . . . . . . . . . . . . . . . 56 A.4.1 Defining Parameters and Importing Data . . . . . . . . . . . . . . . . . . . . . . . 56 A.4.2 Theoretical Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 A.4.3 COMSOL Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 A.4.4 Mesh Convergence Plotting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 A.4.5 Moment Benchmarking Plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 A.5 Axisymmetric Simulation Processing Code . . . . . . . . . . . . . . . . . . . . . . . . . . 61 v A.5.1 Defining Parameters and Importing Data . . . . . . . . . . . . . . . . . . . . . . . 61 A.5.2 Theoretical Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 A.5.3 Simulation Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 A.5.4 Moment Benchmark Plotting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 A.5.5 Mesh Convergence Plotting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 A.6 Ballpark Calculations for 3D Experiments and Simulations . . . . . . . . . . . . . . . . . . 66 A.7 Three-Dimensional Simulation Processing Code . . . . . . . . . . . . . . . . . . . . . . . . 67 A.7.1 Initializing Constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 A.7.2 Simulation Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 A.7.3 Plotting M2 Simulation vs. Theory . . . . . . . . . . . . . . . . . . . . . . . . . . 70 A.7.4 Slope Plotting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 A.8 Intensity to Concentration Calibration Code . . . . . . . . . . . . . . . . . . . . . . . . . . 72 A.9 Experimental Processing Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 A.9.1 Generating Concentration Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 A.9.2 Concentration vs Time, calculation and plot . . . . . . . . . . . . . . . . . . . . . . 75 A.9.3 M2 Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 A.9.4 Trial 1: Aspect Ratio (a/b=0.22) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 A.9.5 Trial 2: Aspect Ratio (a/b=0.48) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 vi List of Figures 1 An initial plug of dye in 2D Poiseuille flow, where viscous forces dominate, will experience Taylor Dispersion and spread over time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Dispersion factor as a function of channel aspect ratio for rectangular channels. The value f = 1 corresponds to flow between parallel plates. . . . . . . . . . . . . . . . . . . . . . . . 8 3 High aspect ratio channels create back-loaded solute distributions, while low aspect ratio channels form front-loaded distributions during the transitional time scale.3 . . . . . . . . . 9 4 Condition of initial Gaussian solute distribution defined by ‘cini’ variable, located at x = xm with standard deviation σ0 defined as COMSOL parameters. . . . . . . . . . . . . . . . . . 15 5 A rectangular geometry element is defined within the settings window to function as the flow domain 2D simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 6 Initial condition and boundary conditions on "Transport of Diluted Species" physics for a two-dimensional channel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 7 Sample gradient mesh. This mesh has 10 elements spanning the height, 40 spanning the width, and an element ratio of 15 in an arithmetic sequence, indicating a linear distribution. . 17 8 Mesh convergence analysis results for two-dimensional channel with a gradient mesh for Pe=1000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 9 Concentration profiles for the 2D channel at two time steps, normalized so the area under each curve is one. The expected "front-loaded" behavior representative of negative skewness is visible. Initial standard deviation of solute is σ0 /a = 4, Pe=1000. The simulation images have their height scale increased by a factor of 100 for visibility. . . . . . . . . . . . . . . . 19 10 Non-dimensional M2 (a), M3 (b) and skewness (c) quantities derived from simulation results at Pe=1000 are superimposed with their corresponding mathematical predictions.2, 4 . . . . . 21 11 Initial condition and boundary conditions for "Transport of Diluted Species" physics for an axisymmetric channel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 12 Mesh convergence analysis results for axisymmetric COMSOL simulation at Pe=1000. . . . 25 13 Concentration profiles in an axisymmetric channel are normalized so the area under each curve is 1. Two time steps are shown, τ = 0.15 and τ = 0.3 for Pe=1000. Initial standard deviation of solute is σ0 /a = 4. The simulation images have their height scale increased by a factor of 100 for visibility. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 14 Axisymmetric benchmarking results. Quantities M2 (a), M3 (b) and skewness (c) overlaid with exact analytical expressions.2, 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 15 Three-dimensional velocity distributions for two aspect ratios in a stationary reference frame. The magnitude of the profiles will depend on the desired average inflow velocity, while the shape depends on aspect ratio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 16 COMSOL results for M2 over non-dimensional time for a 3D rectangular channel at Pe=27200 for aspect ratio 1. The results show agreement to the theoretical prediction at long times for a 3D channel.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 17 Laser-cut edge of PDMS showing surface roughness issues along internal channel wall. . . . 33 18 500 µm wide channel fabricated with unrefined desktop craft-cutter technique. Rough sur- face features on the scale of 60 µm can be seen along the channel walls. Figured adapted from Yuen&Goral.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 19 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 20 With refined techniques, craft-cutters can generate tape and transparency film based chan- nels with wall smoothness comparable to industry leading techniques. This particular chan- nel is 210 µm in width. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 21 Silhouette Cameo 3 Cutting Machine with cutting mat and self-adjusting autoblade. . . . . . 35 vii 22 Double-sided adhesive options tested for formulation of the microchannel walls. . . . . . . . 36 23 Microscope image of channel made with Recollections double-sided sheets. This material was not chosen because removal of the channel negative often left pieces of tape remaining along the channel walls. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 24 Channels sandwiched between glass microscope slides often had leakage problems because of poor tolerancing of the slide’s flatness. . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 25 The polyimide tape should be secured with scotch tape to the cutting board to prevent move- ment during cutting. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 26 Cutting settings available in Silhouette Studio software. . . . . . . . . . . . . . . . . . . . . 39 27 Removal of the channel negative after cut is made. . . . . . . . . . . . . . . . . . . . . . . 40 28 Photograph of experimental setup including microchannel, microscope, and syringe pump. . 42 29 Calibration plot illustrating the relationship between intensity and concentration. Experi- mental results must be taken within the linear regime below 0.25 g/L for accurate conversion. 43 30 Microscope images of the two experimental microchannels. Channel (a) is 458.6 ± 1.0 µm wide. Channel (b) is 209.4 ± 0.5 µm in width. Both channels have a height of 101 ± 3 µm. . 44 31 Experimental solute profiles over the channel height at three different times for rectangular channel with (a/b) = 0.22 ± 0.01, Pe=88.6, with videos taken 2 cm downstream from solute inlet. Experimental frames included for each time step, with brightness doubled for clarity. . 45 32 Experimental data showing intensity over time for rectangular channel with aspect ratio of 0.22±0.01 with asymmetric Gaussian fitted curve (Equation (36). For this trial, Pe=88.6. Videos are taken 2 cm downstream from the solute inlet. . . . . . . . . . . . . . . . . . . . 46 33 Experimental data intensity plot at times t > td overlaid with asymmetrical Gaussian fit curve of form Equation (36). For this trial, (a/b) = 0.48 ± 0.03, Pe=531.58, ta =3.85s and videos are taken 6 cm downstream from the solute inlet. . . . . . . . . . . . . . . . . . . . . . . . 48 34 Experimental values of the dispersion factor, f , overlaid with the theoretical predictions de- rived by Section 2.2. The experimental values agree with the theory within the experimental error. Complete error analysis is covered in Appendix A.3. . . . . . . . . . . . . . . . . . . 49 viii List of Tables 1 Summary of governing variables. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 The dispersion factor, f , for various channel geometries. . . . . . . . . . . . . . . . . . . . 6 3 List of parameters for 2D simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4 List of parameters for 3D simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 5 Thicknesses of channel materials and best cutting settings. . . . . . . . . . . . . . . . . . . 37 ix 1 Introduction 1.1 Introduction to Taylor Dispersion The benefits of down sizing standard laboratory procedures are illuminated by recent developments in the field of microfluidics. Studying fluid-chemical interactions on a micron scale allows scientists to mini- mize costs and increase productivity for applications such as diagnostic testing,6 chemical preparation and more.7 Of interest for many of these applications is the behavior of solute spreading in fully developed laminar flow, first examined by Sir Geoffrey Taylor in 1953.1 Taylor discovered that under these flow con- ditions, an initially localized solute being advected by pressure-driven flow will exhibit enhanced spreading compared to solely molecular diffusion. This phenomenon is now known commonly as “Taylor Dispersion”. Recent studies have more thoroughly quantified this behavior at long time scales and even documented de- pendencies of the dispersing solute profile on channel aspect ratio.2, 3, 8 However, unexplored phenomena relating to this problem are still plentiful. Continued exploration of the topic will enhance understanding and facilitate technological developments. This thesis serves to improve accessibility to ongoing Taylor Dispersion research by simplifying and standardizing numerical and experimental techniques. These opti- mized techniques and protocols allow for streamlined studies into the world of microfluidic dispersion in order to accelerate new discoveries. Numerical studies are performed in commercial software and bench- marked using analytical theory. Further experimental comparison to these results is facilitated by optimizing a low-cost, rapid-prototyping procedure for microchannel fabrication. 1.2 Microfluidics Motivation Advantages to minimizing the scale of laboratory processes include decreases in cost, material, and waste. Microfluidic devices also boast faster results, better repeatability, ease of use and transportability compared to their full-scale counterparts. These are especially worthwhile improvements in the field of biomedicine, where a sample might be blood for a diagnostic test. Microfluidic devices have allowed these procedures to become increasingly affordable and accessible in regions without laboratory resources or personnel.6 An additional advantage of microfluidic analysis is the relative simplicity of fluid motion. Fluid systems at the micron scale tend to experience well-behaved flow with little mixing since it is unlikely to reach a turbulent regime with such small characteristic length scales. The flows are thus easier to quantify, under- stand, and predict. Viscous forces dominate inertial forces, resulting in dynamics known as “creeping” or “Stokes flow,” which have accurate and well-characterized mathematical models. The exploration and doc- 1 umentation of fluid phenomena at small scales is an active area of research, guided by the needs of a variety of industries. 1.3 Discovery of an Effective Dispersion Coefficient The first literature documenting the behavior of solute dispersion in shear pipe flow is credited to Tay- lor, who demonstrates through theory and experiment that solute dispersion is governed by an effective diffusivity coefficient dependent on convective effects in small-scale pipe flow.1 Time u(y) Figure 1: An initial plug of dye in 2D Poiseuille flow, where viscous forces dominate, will experience Taylor Dispersion and spread over time. The considered laminar pipe flow may be assumed for the duration of this thesis to be fully developed and have reached steady state, thus the velocity profile does not change along the length or with time. In order to examine dispersion of a chemical in this flow, an initial plug of solute or dye is introduced near the flow inlet. As shown in Figure 1, the solute is entrained in the flow and propagates downstream as time progresses. Molecular diffusion also begins instantaneously after the solute is introduced, working to min- imize concentration gradients between the dye and surrounding fluid. During the transitional time regime (illustrated by the two middle channels in Figure 1), the dye mimics a profile similar to the surrounding fluid velocity profile. At first the species moves more quickly near the center of the channel than near the walls, where the no-slip condition has a strong influence. This creates local transverse concentration gra- dients in the channel at the front and back end of the advecting solute. Molecular diffusion then acts to equalize these gradients. Solute is diffusively driven towards the channel walls at the leading edge, and towards the center at the trailing edge. The transitional time regime ends when all transverse gradients have effectively been eliminated. Then the solute travels as a simple Gaussian distribution which continues to 2 diffuse predominantly in the axial direction, as shown by the last channel in Figure 1. Looking back to the effective dispersion coefficient proposed by Taylor, it can be seen how this non- uniform shear flow enhances overall dispersivity of the solute versus molecular diffusion alone. In the next section, literature is presented in which this effective dispersion coefficient can be quantified for long time scales for various microchannel geometries. 1.4 Governing Equations and Variables The two systems of physics simultaneously governing a Taylor Dispersion problem are fluid motion and solute spreading. These can be evaluated by partial differential equations known as the Navier-Stokes equations and the advection-diffusion equation, respectively. For ease of understanding and reduction of parameters, each equation is non-dimensionalized. All variables and their definitions are summarized in Table 1. A governing non-dimensional quantity for the fluid flow problem is Reynolds number, Re = Ua/ν = ρUa/µ. This quantity is representative of the flow’s ratio of inertial to viscous forces. In a microfluidic set- ting, the Reynold’s number should be very low (often much lower than 1) due to the dominance of viscous forces. As the Reynold’s number approaches zero, the Navier-Stokes equations for steady, incompress- ible, fully-developed flow simplify to what is known as the Stokes equation. This formula will govern the creeping fluid motion present in Taylor Dispersion studies, ˜ p˜ + ∇ 0 = −∇ ˜ 2 u˜ . (1) In order to examine the solute concentration versus time, the advection-diffusion equation is utilized.9 Beginning from the full equation, ∂C = ∇ · (κ∇C) − ∇ · (uC), (2) ∂t a simplification can be made for incompressible flow which has a divergence of zero, ∇ · u = 0. The equation then becomes ∂C + u · ∇C = κ∇2C ∂t 3 Variable Units Definition u m/s fluid velocity t s time ρ kg/m3 fluid density ∇p Pa axial pressure gradient ν m2 /s kinematic viscosity µ N·s/m2 dynamic viscosity a m channel half-height U m/s average fluid velocity C mol/m3 solute concentration C0 mol/m3 initial solute concentration xm m position of center of mass of initial solute σ0 m standard deviation of initial condition in x-direction x m axial coordinate (cartesian) y m transverse coordinate (cartesian) z m depth coordinate (cartesian) κ m2 /s molecular diffusivity of solute u˜ - non-dimensional fluid velocity ˜ p˜ ∇ - non-dimensional axial pressure gradient x˜ - non-dimensional axial coordinate (cartesian) y˜ - non-dimensional transverse coordinate (cartesian) z˜ - non-dimensional depth coordinate (cartesian) τ - time non-dimensionalized by the diffusive time scale Table 1: Summary of governing variables. which for parallel shear flow in the x-direction reduces to: ∂C ∂C + u(y, z) = κ∇2C. (3) ∂t ∂x This equation governs all motion of the initially localized solute in incompressible flow. Again, each variable is converted to a non-dimensional form. Time is non-dimensionalized with a quantity known as the diffusive 4 time scale a2 td = (4) κ which marks the end of the transitional time regime when all localized concentration gradients have been neutralized. A ‘long time scale’, as often referenced in the literature, refers to a time much greater than td . After non-dimensionalization, Equation (5) is derived: y u˜ = Uu x˜ = x a y˜ = a z˜ = z a t˜ = tκ a2 C˜ = CC0 ˜ = ∇a ∇ ∂ C˜ ∂ C˜ ˜ 2C. ˜ + Pe u(y, ˜ z) =∇ (5) ∂ t˜ ∂ x˜ A single non-dimensional quantity is extracted, the Peclet number; Ua Pe = . (6) κ This quantity represents a ratio of advective to diffusive forces in a flow. The flow problem is always steady, fully-developed and laminar. Within the solute dispersion problem, inertial effects are neglected due to scale. The solute is assumed to have no influence on convection since the molecules are much smaller than the height of the channel.1 Additionally, the solute is assumed to comprise a dilute solution, so fluid properties such as density and viscosity are assumed independent of concentration. In order to reach a unique solution for each partial differential equation, boundary conditions and an initial condition must be specified. Along the channel walls, the no-slip condition must be satisfied for the flow problem, and no-flux for the solute dispersion problem. The initial condition is a patch of solute defined at a distance xm from the channel inlet. The patch is modeled specifically by a Gaussian distribution, defined as 2 Ci (x) = C0 e−0.5 (x−xm )/σ0 (7) 5 for cartesian coordinates with standard deviation σ0 . The choice of C0 , the magnitude of initial solute concentration, is arbitrary from the mathematical perspective, since the problem is linear in concentration. 2 Benchmarking Theory 2.1 Taylor-Aris Dispersivity A review paper completed by Debashis Dutta et al. titled "Effect of channel geometry on solute dis- persion in pressure-driven microfluidic systems", summarizes the foundational research on quantifying the aforementioned dispersion coefficient at long time scales.2 A realized dependence on cross-sectional chan- nel geometry elucidated these numerics. The general form is as follows in Equation (8); 2 Ua 2   K 2 = 1+ f = 1+ Pe 2 f (8) κ 105 κ 105 where K denotes the effective dispersion coefficient or ‘Taylor-Aris’ dispersivity. Lastly, f is a geom- etry dependent factor which considers the influence of channel walls, known as the dispersion factor. The dispersion factor may be presented as a constant for particular geometries of interest, and is computed for three-dimensional geometries in Section 2.2.2 Channel Geometry f 2D Channel 1 Circular Cylinder 1.09 Square Channel 1.76 Table 2: The dispersion factor, f , for various channel geometries. Using a known dispersion factor, fluid properties and channel dimensions, an exact expression for Taylor-Aris dispersivity can be obtained. This expression will provide the theoretical values for bench- marking procedures in both the 3D numerical simulations and physical experiments presented in subsequent sections by the relationship that σ 2 (t) K lim 2 = 2τ (9) t→∞ a κ for asymptotic times where σ is the standard deviation of the solute distribution in the axial direction. Experimentally, the dispersion factor, f , is of great interest because so many variations on dispersion type 6 studies remain to be explored. A complete understanding of Taylor-Aris Dispersivity relies on continued characterization in varied microfluidic environments, such as channels with altered wall conditions.2, 9 2.2 Dispersion Factor for Rectangular Channels For numerical and experimental benchmarking purposes, the values of the dispersion factor for rectan- gular geometries are desired. These may be calculated following the Dutta et al review paper2 and utilizing the partial differential equation solver available within COMSOL Multiphysics. The geometry dependent dispersion factor is a function of channel aspect ratio, specifically, 840 u∗ Z f (a/b) = ∗ ∗ g dA∗ (10) A D∗ U A A∗ = a2 where u∗ is non-dimensionalized velocity, µL u∗ = −u 4a2 ∆P U ∗ is unitless average velocity, 1 Z ∗ U = ∗ u∗ dA∗ A D∗ A is the channel’s cross-sectional area, D∗ is the channel’s cross section, δ D∗ is the cross section boundary, and g is the first moment of solute concentration at long times. For a rectangular cross section with a center at zero, the velocity terms are known and may be expressed as:   2  ∞     1 y 32 y 1 cosh(nπx/2a) u∗ = 1− − ∑ 3 3 sin nπ + , 8 a n=1,3,5... n π 2a 2 cosh(nπb/2a) ∞     ∗1 16 a nπb U = − ∑ tanh . 2 n=1,3,5... n5 π 5 b a The function g must be evaluated as a partial differential equation with boundary conditions: u∗ ∇∗2 g = 1 − U∗ 7 Figure 2: Dispersion factor as a function of channel aspect ratio for rectangular channels. The value f = 1 corresponds to flow between parallel plates. constrained by ∗ ∇ g · nˆ =0 δ D∗ where nˆ is the vector normal to the boundary and Z gdA∗ = 0. D∗ Evaluation of this system in COMSOL Multiphysics results in Figure 2 that captures the dispersion factor constants for various aspect ratios, (a/b), of rectangular channels. 2.3 Longitudinal Asymmetries in Transitional Time Regimes Considering the effect of channel geometry on dispersion at long time scales, one may presume that an influence exists during the transitional regime as well. Explorations into this topic are relatively new in the literature, leaving opportunities for continued research. Recently, researchers at the University of North Carolina Chapel Hill uncovered a relationship between aspect ratio and longitudinal asymmetries in the advecting solute profile.3 Aspect ratio is defined as half-height over half-width, a/b, equivalent to height over width. Channels with a high aspect ratio, such as an axisymmetric cylinder, spread the solute with a 8 Figure 3: High aspect ratio channels create back-loaded solute distributions, while low aspect ratio channels form front-loaded distributions during the transitional time scale.3 “back-loaded” profile that has a peak upstream from center. Low aspect ratio channels, such as short and wide rectangular prisms, deliver the solute with a sharp front and long tapering tail i.e. a “front-loaded” distribution. This distinction is illustrated in Figure 3. A mathematical description of longitudinal asymmetries in the translating solute profile can be devel- oped using moments of the concentration distribution. The nth moment is defined by 1 Z ∞ Mn = R ∞ ¯ (x˜ − x˜m − Peτ)n C( ¯ x, ˜ τ)d x. ˜ (11) ˜ τ)d x˜ −∞ C(x, −∞ ¯ is given by The cross sectionally averaged concentration, C, ¯ =1 ZZ C(x) C(x, y, z) dA, (12) A D and x˜m is the initial length-wise location of the solute plug. Qualitatively, the first moment represents the location of the center of mass of the solute, which is known to move with the average velocity. The second moment is the variance i.e. the squared standard deviation of the distribution; a good measurement of the axial width of the spreading distribution.10 At long time scales, the slope of the second moment is proportional to the dispersion coefficient. The third moment measures asymmetries in the distribution. Normalizing the third moment gives skewness; M3 Skewness = 3/2 . (13) M2 9 Typically, positive skewness indicates a back-loaded distribution while negative skewness is the con- verse.3 Francesca Bernardi derived exact analytical expressions for variance and third moment for all time of the solute concentratino distribution.4 Equations (14) & (15) are for a two-dimensional channel and Equations (16) & (17) are for an axisymmetric channel: ∞ −(nπ) τ2 2 4 e M2 (τ) = σ02 + 2τ − Pe2 + Pe2 τ + 36Pe2 ∑ 8 (14) 105 105 n=1 (nπ) 2  27 3 ∞ e−(nπ) τ  172 3 8 3 1488 M3 (τ) = Pe − Pe τ + Pe ∑ 10 144 − 24τ − (15) 716625 5775 8 n=1 (nπ) (nπ)2   ∞ −βn2 τ 1 1 e M2 (τ) = σ02 + 2τ + Pe 2 − 2 + τ + 128Pe ∑ 8 .2in (16) 360 24 n=1 βn   ∞   3 17 1 3 −βn2 τ 240 18 1 M3 (τ) = Pe − + τ + 128Pe ∑ e − 12 + 10 + 8 τ , (17) 53760 480 n=1 βn βn βn where the symbol β denotes zeros of the Bessel function. The initial dye patch is taken as a Gaussian dis- tribution with standard deviation σ0 at instantaneous time zero for ease of translation into the experimental realm. The second and third moment described here are used for 2D numerical simulation benchmarking discussed herewithin, since simulation accuracy may be checked over all time. At long times, the second moment expressions should match the asymptotic predictions from Dutta et al.2 Considering a large limit of τ, Equation (14) for a two-dimensional channel simplifies to: 4 lim M2 (τ) = 2τ + Pe2 τ. (18) τ→∞ 105 From Dutta et al,2 the second moment for a two-dimensional channel in the long time regime is deduced by Equations (8) & (9), which gives the same expression. For a circular pipe, Equation (16) in the long time regime becomes 1 2 lim M2 (τ) = 2τ + Pe τ. (19) τ→∞ 24 10 The asymptotic expression from Dutta et al may be approximated by Equations (8) & (9) while using f = 1.09. However, a more exact expression is defined for elliptical channels that may replace Equation (8), 24 − 24e2 + 5e4   K 1 = 1 + Pe2 κ 48 24 − 12e2 p where e = 1 − a2 /b2 . So for a circular cylinder, e = 0 and the expression reduces to Equation (19) identically. 11 3 Optimized Numerical Simulations for Microfluidic Dispersion More than ever, researchers are turning to numerical simulation software to make light work of com- plex experimental ideas. Recent microfluidic dispersion simulations include electro-osmotic pipe flow with slip-stick wall conditions,11 varying Peclet number flow through Cassie-Baxter state hydrophobic walls,12 and chemical tracing of particles through a micromixer.13 Such simulations often compose the very first documentation of dispersion behavior in new settings, making them a vital tool for researchers. Numerical simulations are also a useful tool for determining the viability of experiments before spending time and money on elaborate experimental setups. Computational Fluid Dynamics (CFD) is one of many features in the user-friendly, FEM-based multi- physics simulation program; COMSOL Multiphysics. This software is excellent for acclimating new users, even with no prior simulation experience. Users have the ability to solve intricate, coupled multiphysics problems across numerous fields including but not limited to: electromagnetics, solid mechanics, heat trans- fer and fluid dynamics. In terms of user accessibility, COMSOL is a good choice for Taylor Dispersion simulations. Succeeding methods further facilitate these studies by minimizing computation time and mem- ory usage within COMSOL so that runs are executable on a single workstation. The impressive accuracy of results derived by Taylor Dispersion type COMSOL models is supported by enclosed mesh convergence analyses and benchmarking results of simulations created to replicate the work in Section 2.3. 3.1 Two-Dimensional Channel 3.1.1 Problem Setup A two-dimensional channel represents flow between parallel flat plates with infinite depth into the page. For the duration of this study, the fluid is assumed incompressible. The problem may be modeled with one-way coupled physics, assuming the solute comprises particles with a small scale compared to the height of the channel and does not affect fluid properties or flow. The Peclet number is selected to be 1000 in this section, representative of common microfluidic experiments. The flow is pressure driven, steady state, fully developed, and thus can be modeled by an exact solution to the Navier-Stokes equations. The calculation begins with the steady state Stokes equation where inertial effects are neglected, 1 0 = − ∆p + ν∇2 u. (20) ρ Rearranging and taking only the x-component gives 12 ∂ 2u 1 ∂p 2 = . ∂y µ ∂x To arrive at an expression for u(y), the expression is integrated twice; ∂u 1 ∂p = y +C1 ∂y µ ∂x 1 ∂p 2 u(y) = y +C1 y +C2 2µ ∂ x and boundary conditions are utilized to solve for the integration constants. There is no-slip at y = −a, a, 1 ∂p u(−a) = 0 = (−a)2 −C1 a +C2 2µ ∂ x 1 ∂p 2 0= a −C1 a +C2 (21) 2µ ∂ x 1 ∂p 2 u(a) = 0 = (a) +C1 a +C2 2µ ∂ x 1 ∂p 2 0= a +C1 a +C2 . (22) 2µ ∂ x Subtracting Equation (21) from Equation (22) yields 0 = 2C1 a so C1 = 0. To find C2 , similar proceedings follow with Equation (21), 1 ∂p 2 0= a +C2 2µ ∂ x a2 ∂ p C2 = − . 2µ ∂ x 13 The expression for velocity is 1 ∂p 2 u(y) = (y − a2 ). (23) 2µ ∂ x Calculating the average of Equation (23) will allow an expression without pressure so the fluid flow may be driven by a mean velocity, U, Z a 1 1 1 ∂p Z U= u(y)dA = (y2 − a2 )dy A A 2a −a 2µ ∂ x a2 ∂ p U= . (24) 3µ ∂ x Thus, the final expression for velocity, considering a negative pressure gradient, is given by 3 y2 u(y) = U(1 − 2 ) . (25) 2 a 3.1.2 COMSOL Setup For readability, COMSOL setup instructions are ordered to follow typical workflow: defining parameters and definitions, geometry, physics, boundary conditions, mesh creation and then study settings. Optimiza- tion options are included throughout. Parameters: Within the ‘Model Builder’ window under ‘Global Definitions’, define parameters per Table 3. Definitions: To define the initial Gaussian distribution of the solute, define a variable cini under ‘Compo- nent - Definitions’ by Equation (7).14 Employing this technique results in a smooth Gaussian distribution at time t = 0. 14 Parameter Value a 5x10−5 m L 0.05 m U 0.005 m/s µ 0.001 Pa·s κ 2.5x10−10 m2 /s xm 2.5x10−4 m σ0 2x10−5 m C0 1 mol/m3 Table 3: List of parameters for 2D simulation. Figure 4: Condition of initial Gaussian solute distribution defined by ‘cini’ variable, located at x = xm with standard deviation σ0 defined as COMSOL parameters. Geometry: For the geometry, define a rectangle of height ‘a’ and length ‘L’ positioned with the top surface along y = 0. This is only the bottom half of the channel, as reflectional symmetry is exploited to save computation time. 15 Figure 5: A rectangular geometry element is defined within the settings window to function as the flow domain 2D simulations. Physics: The physics that are present in a microfluidic dispersion problem are fluid motion and solute dispersion, governed by Navier-Stokes equations and an advection-diffusion equation. Within COMSOL, the corresponding physics titles are "Creeping Flow" and "Transport of Diluted Species". During trial sim- ulations using both of these physics, it is clear that the flow problem takes a bulk of total computation time despite being steady state and identical at each channel cross section. To eliminate this lengthy computation, one may utilize solely "Transport of Diluted Species" physics and input the exact flow solution, Equation (25), as the x-direction velocity field under ‘Transport Properties’. The resulting simulation takes a frac- tion of the run time and may save the user money since purchase of the "Microfluidics" COMSOL package containing creeping flow is not necessary. With the physics specified, initial and boundary conditions may now be defined. The variable ‘cini’ should be added in ‘Initial Values’ settings under ‘Concentration’ to appear at instantaneous time zero. Boundary conditions should be applied as in Figure 6. The symmetry condition cuts computation time in half. 16 Figure 6: Initial condition and boundary conditions on "Transport of Diluted Species" physics for a two-dimensional channel. Mesh: Finite element method necessitates discretization of a problem’s spatial domain, specifying solution nodes along the geometry. Since the output of interest in this problem is solute concentration, the resolution must be sharper upstream, where the solute is tightly packed. Once the solute has been smeared downstream, accurate results are achievable with a coarser mesh. As always, designing a mesh is a balance of generating effective results with minimal mesh elements to save computation time. For Taylor Dispersion problems, this optimization may be achieved by building a gradient mesh with rectangular mesh elements. COMSOL has the capacity for building a gradient mesh using the ‘Distribution’ and ‘Mapped’ mesh features. Figure 7: Sample gradient mesh. This mesh has 10 elements spanning the height, 40 spanning the width, and an element ratio of 15 in an arithmetic sequence, indicating a linear distribution. Executing a mesh convergence analysis (Figure 8) varying element size in both the transverse and lon- gitudinal directions determined the most efficient settings. For the analysis, the third moment of solute concentration (Equation (11)) is considered since it is the most difficult calculation for the simulation to match. The variable ‘N’ dictates the number of elements transversing the channel height, while 100N is the number of elements spanning the channel length. 17 Figure 8: Mesh convergence analysis results for two-dimensional channel with a gradient mesh for Pe=1000. As indicated by the mesh convergence analysis, the most effective mesh is 60 elements across the chan- nel height and 6000 elements across the channel length. The third parameter in COMSOL’s ‘Distribution’ mesh feature is ‘Element Ratio’, which should be set to 50. This sets the last element in the pattern as 50 times longer than the first element in the pattern. In these studies, an arithmetic sequence is used, which re- sults in a linear distribution. A geometric sequence could also be used, which is an exponential distribution. Study: A Time-Dependent study is required to see how the solute disperses over a transitional time scale. The study begins at time zero and runs to an end time of five seconds with 0.025 second incremental steps. COMSOL’s time-dependent solver determines time steps necessary for study convergence based on relative tolerance independently of user-input time steps, which only dictate the displayed output values. For this reason, a time step convergence analysis is not necessary as it may be in other CFD softwares. To explore the transitional time, times must stay within the "diffusive time scale", td , defined in Equation (4). The value 18 of td in this case is 10 seconds, placing the full study within the diffusive scale. Once setup is complete, the total runtime of this simulation is only 27 minutes on a 2018 Mac Mini with a 3.2 GHz Intel Core i7 processor and 64 GB of memory. 3.1.3 Results The optimized numerical simulation succeeds in replicating expected results for a low aspect ratio chan- nel. During the diffusive time scale, longitudinal asymmetries in the solute concentration profile are present, showing a front-loaded distribution. The solute is delivered with a sharp front end and a long tapering tail indicative of negative skewness (Figure 9).3 Within COMSOL, a linear projection function takes a cross- sectional concentration average so that a one-dimensional concentration profile may be examined. Figure 9: Concentration profiles for the 2D channel at two time steps, normalized so the area under each curve is one. The expected "front-loaded" behavior representative of negative skewness is visible. Initial standard deviation of solute is σ0 /a = 4, Pe=1000. The simulation images have their height scale increased by a factor of 100 for visibility. 19 Benchmarking the full spectrum of results from the Aminian et al paper entails a comparison of the second moment, third moment, and skewness of solute concentration.3 To complete this examination, a constant scaling factor is necessary to convert the literature’s characteristic velocity Up to the mean velocity used here, U. The constant factor, which is dependent on channel geometry (calculated in Appendix A.1) is 3 Up = U. (26) 2 This scaling factor is used in all subsequent post-processing procedures. Figure 10 shows finalized benchmarking results with all moments non-dimensionalized by a unit-proportional characteristic height and time by td . Impressive agreement is apparent between the optimized numerical simulation and the mathematical model.4 Skewness is everywhere negative, supporting the evidence of a front-loaded distri- bution profile. This benchmarking agreement illustrates how effective COMSOL is as a tool for solving microfluidic dispersion problems quickly and accurately for a two-dimensional domain. 20 Figure 10: Non-dimensional M2 (a), M3 (b) and skewness (c) quantities derived from simulation results at Pe=1000 are superimposed with their corresponding mathematical predictions.2, 4 21 3.2 Axisymmetric Channel 3.2.1 Problem Setup To benchmark the full spectrum of results from Section 2.3, a high aspect ratio channel is considered next; a circular cylinder of aspect ratio 1. The flow is again modeled as a steady state solution to the Navier-Stokes equations. Beginning from the steady Stokes equation for creeping flow, 1 0 = − ∇p + ν∇2 u (27) ρ in cylindrical coordinates, taking the z-direction, this equation becomes   1 d du 1 dp r = . r dr dr µ dz To arrive at an expression for u(y), the expression is integrated twice, du 1 d p r2 r = +C1 dr µ dz 2 du 1 d p r C1 = + dr µ dz 2 r 1 d p r2 u(r) = +C1 ln(r) +C2 . µ dz 4 Utilizing the boundary conditions, u(r = 0) = constant and u(r = a) = 0, u(0) = constant = C1 ln(0) +C2 , the result of ln(0) is - ∞, thus C1 must be equal to zero. Calculating C2 , 1 d p a2 u(a) = 0 = +C2 µ dz 4 gives d p  a2 C2 = − . dz 4 22 Consequently, the expression for velocity in terms of a pressure gradient is 1 dp 2 u(r) = − (a − R2 ). 4µ dz In order to derive the equation without a pressure gradient term, an average velocity is calculated. Z a 1 1 1 dp 2 Z (a − R2 ) dR  U= u(y)dA = 2πR − A A πa2 0 4µ dz simplified to be 1 a2 d p U =− . 8 µ dz The final velocity expression is: r2   u(r) = 2U 1 − 2 . (28) a 3.2.2 COMSOL Setup When selecting the dimensionality of a new study in COMSOL, users are offered an option for ax- isymmetric studies called "2D-Axisymmetric". This feature solves problems with rotational symmetry as two-dimensional and then automatically revolves the solution for any desired results. This tool should be used whenever possible as it saves significant computation time. Since an axisymmetric study modeled in two dimensions is nearly identical to the previous two-dimensional simulation, only key differences in the COMSOL setup will be highlighted. Parameters: Parameters for this simulation are identical to those described in Table 3 except for ‘xm ’, which is now ‘zm’ in cylindrical coordinates with a value of zm = 5x10−4 meters. The axial coordinate is now ‘z’ and the radial coordinate is represented by ‘r’. Definitions: The expression defining the initial solute’s Gaussian profile is transformed to cylindrical coordinates, 2 cini = C0 e−0.5 (zm−z)/σ0 . (29) Geometry: A single rectangle with length ‘L’ along the z-direction and height ‘a’ along the r-direction 23 should be placed along the axis of rotational symmetry at r = 0. Physics: "Transport of Diluted Species" may be used independently once again, with Equation (28) entered under ‘Transport Properties’ to describe the velocity field. Initial conditions are defined by Equation (29). Inflow and outflow boundary conditions define the channel inlet and outlet, while a no-flux condition im- posed on the internal wall prevents solute from passing through that boundary. Lastly, an ‘Axial Symmetry’ boundary condition is defined on the wall along the line r = 0 as shown in Figure 11. Figure 11: Initial condition and boundary conditions for "Transport of Diluted Species" physics for an axisymmetric channel. Mesh: The convenience of a 2D-Axisymmetric study is especially apparent during mesh creation. The same form of a gradient mesh with rectangular elements from the previous study may be utilized. A con- vergence analysis is completed in this case to ensure accurate results are being obtained, results shown in Figure 12. The optimal mesh has 50 rectangular elements spanning the channel radius and 10,000 spanning the channel length with an element ratio of 50 defining the gradient as an arithmetic sequence. 24 Figure 12: Mesh convergence analysis results for axisymmetric COMSOL simulation at Pe=1000. Study: A time-dependent study is programmed with start time 0 seconds, end time 5 seconds, and 0.025 second increments. 3.2.3 Results For a high aspect ratio channel, the literature predicts a back-loaded distribution indicative of positive skewness during the diffusive time scale. The numerical simulation succeeds in matching this prediction. Concentration profiles for two times within the transitional regime in Figure 13 illustrate the initially gradual delivery of solute followed by a high peak of concentration at the tail. More precisely, the second moment, third moment and skewness of solute concentration profile are compared to exact analytical expressions for the transitional time regime. The M2 plot is also overlaid with the long time asymptotic prediction, which simulation results should approach over time. The results in Figure 14 support the capabilities of COMSOL as an accurate numerical tool for solving solute dispersion problems within fully laminar flows. The conversion used between the characteristic ve- locity used in development of the analytical model and the average velocity used in COMSOL (calculated 25 Figure 13: Concentration profiles in an axisymmetric channel are normalized so the area under each curve is 1. Two time steps are shown, τ = 0.15 and τ = 0.3 for Pe=1000. Initial standard deviation of solute is σ0 /a = 4. The simulation images have their height scale increased by a factor of 100 for visibility. in Appendix A.2) is Up = 2U. (30) 3.3 Three-Dimensional Channel 3.3.1 Problem Setup A three-dimensional simulation is modeled with physically reproducible parameters and conditions so that it may serve to both benchmark theoretical predictions and compare with experimental data. Actual Taylor Dispersion experiments are most often concerned with finding the dispersion coefficient at long time scales, which may be compared with the expected value from Section 2.1 for a given channel geometry. 26 Figure 14: Axisymmetric benchmarking results. Quantities M2 (a), M3 (b) and skewness (c) overlaid with exact analytical expressions.2, 4 27 For the first simulation, a high aspect ratio channel with a square cross section is considered. The flow solution for a rectangular channel of aspect ratio (a/b) is given by an exact analytic expression. To derive this expression, the starting point is the equation for pressure driven flow through a rectangular channel, given by the supplementary materials of the Aminian et al paper.3 This expression considers a reference frame moving with the flow, thus velocity is denoted as uL (y, z) for Lagrangian reference frame. −a2 ∂ p uL (y, z) = m(y, z) (31) µ ∂x where m(y, z) is given by  2 ∞ 4(−1)k        1 y 1 πy 1 πz m(y, z) = − +∑ 3 cos k− cosh (k − +... 3 a k=1 ((k − 1/2)π) cosh((k − 1/2)πb/a) 2 a 2 a (−1)k sinh((k − 1/2)πb/a)  ... + (32) ((k − 1/2)π)2 In order to convert this velocity expression to a stationary reference frame, the average flow velocity is added, −a2 ∂ p u(y, z) = m(y, z) +U. µ ∂x The average velocity, U, is evaluated by −a2 ∂ p Z bZ a 1 U= m(y, z)dydz ab −b −a µ ∂x −a2 ∂ p U= c(a/b) µ ∂x which introduces a function, c(a/b), dependent on the channel aspect ratio. Thus, from a stationary refer- ence frame, 28 −a2 ∂ p u(y, z) = m(y, z) +U µ ∂x and equivalently, −a2 ∂ p   u(y, z) = m(y, z) + c(a/b) . µ ∂x This thesis wishes to examine a flow that is dictated by an average inlet velocity as opposed to a pressure gradient. By rearranging the above expressions, the pressure gradient term may be replaced, resulting in;   m(y, z) u(y, z) = U 1 + . (33) c(a/b) The exact value of c(a/b) is determined by satisfying the no-slip boundary condition at the channel walls, which should hold velocity at zero for a stationary reference frame. Thus,   m(y, z)|wall 0 =U 1+ c(a/b) which rearranges to c(a/b) = −m(y, z) (34) y=a, z=b The evaluation of c(a/b) should generate a constant that is unique to a given aspect ratio, resulting in profiles like Figure 15. 29 (a) (a/b) = 0.22 (b) (a/b) = 0.48 Figure 15: Three-dimensional velocity distributions for two aspect ratios in a stationary reference frame. The magnitude of the profiles will depend on the desired average inflow velocity, while the shape depends on aspect ratio. 3.3.2 COMSOL Setup A three-dimensional model will be used to study Taylor Dispersion through a rectangular channel within COMSOL. The flow travels lengthwise down the channel in the x-direction, the height of the channel is in the y-direction, and the depth of the channel in the z-direction. Parameters: The necessary parameters are identical to the two previous studies, now with different values as defined in Table 4 to mimic plausible experimental quantities. Parameter Value a 5x10−4 m b 5x10−4 m L 10.88 m U 0.031008 m/s µ 8.9x10−4 Pa· s κ 5.7x10−10 m2 /s xm 0.005 m σ0 0.001 m C0 1 mol/m3 Table 4: List of parameters for 3D simulation. Definitions: The variable necessary to define an initially localized Gaussian distribution of solute is iden- tical to the two-dimensional case, defined in Equation (7). 30 Geometry: A block with length L, height a, and width b fully defines the geometry. The corner should be located at (0,0,0). This build generates a quarter of the channel, since symmetry will be exploited later on. The model is theoretically solving for a channel with total height 2a and width 2b. Physics: One “Transport of Diluted Species” physics sufficiently models solute transport when an exact solution to the laminar flow in a rectangular channel (Equation (33)) is input as the transport property. Symmetry conditions should be selected for the surfaces along the y = 0 and z = 0 planes. An inflow and outflow condition should be defined on the ends of the channel, and a no-flux boundary condition must be imposed on the remaining channel walls. Mesh: For the three-dimensional study, full optimization is of even greater importance than before due to the large size of matrix equations which must be solved within the numerical software. The most important step for minimizing computation expense is mesh selection. A gradient mesh with rectangular elements is selected once again, but this time is distributed as a ‘Geometric Sequence’ which generates an exponential distribution. This allows for less mesh elements to be defined along the channel length, while still achieving the fine mesh resolution required near the channel inlet. This type of mesh would also work for the two- dimensional and axisymmetric channels, but is especially important for the three-dimensional channels as it was found to save computation time. The optimal mesh has 15 elements spanning both the height and depth, 2000 elements spanning the channel axially, and an element ratio of 700. Study: This study will begin at time 0 and run for 200 seconds with increments of 5 seconds. The diffusive time scale ends at td = 439 seconds so the study will reach about τ = 0.46, which is enough to benchmark the desired results. 3.3.3 Results The goal of the three-dimensional simulations is to match theory from Section 2.1. The effective dis- persion coefficient predicted in the review paper for a given channel geometry and aspect ratio is related to the second moment of solute concentration at long time scales t  td by Equations (8) & (9). Through the diffusive time scale, M2 will gradually evolve to a linear behavior which should match the theoretical prediction at long times. Figure 16 shows the agreement between theory and the COMSOL simulation. 31 Figure 16: COMSOL results for M2 over non-dimensional time for a 3D rectangular channel at Pe=27200 for aspect ratio 1. The results show agreement to the theoretical prediction at long times for a 3D channel.2 32 4 Optimized Craft-Cutter Fabrication Technique 4.1 Introduction to Microchannel Fabrication Typical experimental methodology for microfluidic systems involves a channel fabrication technique known as “Soft Lithography” which requires expensive equipment, highly trained personnel, and a clean- room environment. Another fabrication technique involves laser cutting thinly spin-coated polydimethyl- siloxane (PDMS) on an acrylic plate and then bonding additional PDMS layers to the top and bottom using a plasma cleaner.15 Once again, specialized equipment and trained personnel are needed for the process. Fur- thermore, laser cutting often generates undesirable wall roughness from burns and melting in the material. Figure 17: Laser-cut edge of PDMS showing surface roughness issues along internal channel wall. Researchers without access to these resources may hire a third party for channel fabrication, resulting in a long lead time and inflated costs. Within this thesis, a fabrication method utilizing an affordable and widely available desktop craft-cutter is presented to bypass exorbitant costs and wait time. This technique has been proposed in the past,5 but lacked the necessary wall smoothness for many microfluidic experiments (Figure 18). With the included optimization techniques, surface roughness is comparable to leading fabrication techniques in channels with a minimum channel width of about 200 µm (Figure 20). 33 Figure 18: 500 µm wide channel fabricated with unrefined desktop craft-cutter technique. Rough surface features on the scale of 60 µm can be seen along the channel walls. Figured adapted from Yuen&Goral.5 Figure 19 (a) Fluorescent image. (b) Bright field image. Figure 20: With refined techniques, craft-cutters can generate tape and transparency film based channels with wall smoothness comparable to industry leading techniques. This particular channel is 210 µm in width. 4.2 Optimization of Desktop Cutter Technique The enclosed methods have been carried out using the Silhouette Cameo Model 3 desktop cutting ma- chine, available for purchase under $200. This model features a self-adjusting autoblade which allows for 34 increased customization of cuts in materials of various thicknesses. Cutting patterns may be designed using the included "Silhouette Studio Design" software, or imported as a DXF file. Silhouette Studio Design v4.1 Basic Edition is used here, which allows dimensions as small as 250 micrometers to be specified within internal drawings. If smaller features are required, he design should be drawn in an external program and then imported. Figure 21: Silhouette Cameo 3 Cutting Machine with cutting mat and self-adjusting autoblade. The cutting procedure begins by adhering the desired material to the cutting mat, aligning the cutting mat along the machine, and loading it. Within the Silhouette software, cutting settings may be specified before choosing the "Send" button which initiates the print. The cutting machine calibrates the blade and makes the cut. Once the cutting mat is unloaded, the cut design is available to the user. For the first iterations of desktop-cutter produced microchannels, straight rectangular channels are fab- ricated. The optimization goals include: internal wall smoothness comparable to soft lithography PDMS devices, consistent channel width, and consistent channel height along the channel length. The final device should be water tight so that fluid doesn’t leak between device layers, and be highly reproducible. Lastly, experiments with these channels should successfully benchmark expected theoretical predictions for the effective dispersion coefficient and dispersion factor at long times from Sections 2.1. 35 4.2.1 Material and Cutter Setting Optimization Appropriate material selection is paramount for fabrication of a high quality channel. The material must be thin and flexible enough to be through cut by the craft cutter. The material must also be adhesive on both sides so that it properly seals an upper and lower layer to form a complete device. A number of low cost, commercially available, double-sided adhesive sheet and tape options are tested to find an optimal selection (Figure 22). Figure 22: Double-sided adhesive options tested for formulation of the microchannel walls. Due to the varying thicknesses and tackiness of each tape, the cutting settings must be adjusted for each to achieve an optimal cut. By comparing the best channels from each material, a final material selection is made. Silhouette Studio offers the cutting settings: "Speed", "Force", "Passes" and "Razor". The razor setting is intended to scale with material thickness, as it dictates the extension distance of the blade from the razor housing. Interestingly, this setting drastically influences the resolution of each cut and plays a key role in cutting a sharp edge in any material. The channel quality is always best with only one cut pass at the lowest speed setting, so these settings remain unchanged while the other two are toggled. The optimal cutting settings for each material are listed in Table 5 along with material thicknesses. The best performing material for a straight microchannel is the polyimide double-sided tape. The chan- nel negative removes cleanly and easily with tweezers and this tape has the added benefit that it is manu- 36 Tape Thickness Razor (1 → 10) Force (1 → 33) Recollections Adhesive Sheets 100 µm 10 1 Silhouette Adhesive Sheets 60 µm 9 1 Polyimide Double-Sided Tape 101.6 µm 9 1 Scotch Double-Sided Tape 67 µm 8 1 Aleene’s Instant Tacky 81 µm 10 2 Recollections Double-Sided Tape 202 µm 10 4 Table 5: Thicknesses of channel materials and best cutting settings. factured to a consistent thickness, 1 mil or 101.6 microns. This quality will ensure a channel with uniform height. This tape is not so tacky that it is challenging to handle and cut, but enough so that it strongly adheres a top and bottom device layer. Both glass microscope slides and polyester film were tested as the top and bottom layers of each device. The polyester film outperformed the glass slides since their flexibility allows for air bubble removal during device assembly, preventing leakage between device layers. Figure 23: Microscope image of channel made with Recollections double-sided sheets. This material was not chosen because removal of the channel negative often left pieces of tape remaining along the channel walls. 37 Figure 24: Channels sandwiched between glass microscope slides often had leakage problems because of poor tolerancing of the slide’s flatness. 4.2.2 Fabrication Instructions 1. Power on craft cutter and connect to computer. Open Silhouette Studio to the desired channel design file. 2. Cut a sufficiently long piece of polyimide double-sided tape of 101.6 µm thickness with scissors. 3. Secure polyimide tape to the cutting board, sticky side up, by taping down both edges with scotch tape. This ensures that the tape does not move at all during the cut. 4. Prepare the top and bottom layer of transparency by cutting two rectangles of similar size to the tape. Set aside. 5. Align the cutting board appropriately with the craft cutter and select ‘Load’ on the cutter screen. 38 Figure 25: The polyimide tape should be secured with scotch tape to the cutting board to prevent movement during cutting. 6. Choose the ‘Send’ tab in Silhouette Studio, select your channel design and choose ‘Cut Outline’ to make one continuous cut. Adjust the settings to: Razor=9, Speed=1, Force=1, and Passes=1. The razor setting is graphically represented as a dial before the other settings, shown in Figure (26). Press ‘Send’ to begin the cut. Figure 26: Cutting settings available in Silhouette Studio software. 7. Once the cut is complete, choose ‘Unload’ from the cutter screen. Carefully remove the channel negative in one piece with tweezers (Figure 27). 39 Figure 27: Removal of the channel negative after cut is made. 8. Place one rectangle of transparency over the tape, ensuring that the entire channel is covered. Smooth out all air bubbles with a smooth but finely tipped tool that will not scratch the transparency. Remove the device from the cutting board. 9. Take the other transparency rectangle, secure on the cutting board and load. Cut the inlet, outlet, and solute holes in the transparency with settings: Razor=9, Speed=1, Force=33, Passes=1 and unload. 10. Peel adhesive covering off of polyimide tape, align inlet and outlet holes on transparency slide to channel and adhere final layer. Press out air bubbles with smoothing tool. The channel is complete. 40 4.3 Experimental Validation 4.3.1 Methodology To visualize channel resolution and dispersion behavior, a solution of 0.581 g/L fluorescein sodium salt in distilled water is used as the solute. Pure distilled water is the working fluid. The dispersivity of the fluorescein solution is 5.7x10−6 cm2 /s and it has peak excitation and emission wavelengths of 460 nm and 515 nm, respectively. The fluorescein is illuminated by a mercury lamp coupled with a blue light filter cube of excitation wavelength 480 nm with a 40 nm bandwidth and emission wavelength of 535 nm with 50 nm bandwidth. A neutral density (ND) filter is activated to negate potential photobleaching effects. The inverted microscope is a Nikon Eclipse TE2000-U with a high contrast fluorescence objective installed of 10x magnification. Images and videos are captured via a monochrome Mako U-503B camera with 5 megapixels of resolution. Resolution of experimental data is 2240 pixels/mm. A Harvard Apparatus syringe pump (Standard Infuse/Withdraw Pump 11 Elite) from Harvard Apparatus is used to provide a known inlet velocity to the channel, connected by capillary tubing that is epoxied to the channel inlet along with a tape gasket to form an air-tight seal. An initial solute distribution is introduced by adding a drop of fluorescein solution to the solute inlet, sealing with tape, and then allowing it to diffuse to approximately 2 mm in width. The solute inlet has a distance of xm = .01 m from the laminar flow inlet. Videos are captured at one fixed location per trial, recording the solute as it travels downstream. Solute intensity calibration is necessary to verify a linear relationship between brightness of illuminated fluorescein and solute concentration values, illustrated in Figure 29. To deduce this relationship, 16 different concentrations of fluorescein and distilled water were prepared. The experimental microchannel was flooded with one solution at a time and then a video is captured at a fixed location along the channel. Each video is processed by a MATLAB code, included in Appendix A.8, that reads the brightness intensity within the channel. The mean and standard deviation of fluorescein intensity is recorded by the code, considering each time frame of the video. In between trials of the different solute mixtures, the microchannel is cleaned by rinsing with distilled water and then flushing with air. The camera and lighting configurations were held constant throughout these trials at an exposure time of 100 ms (specified within MATLAB’s image acquisition tool by setting the exposure time parameter to 100,000) and frame rate of 10 frames per second. The mean intensity values of each concentration are documented in Figure 29, along with errorbars which represent two standard deviations of the intensity. As expected, there is a linear regime for lower fluorescein concentrations, and then the behavior becomes non-linear. At high enough concentrations (above 0.55 g/L), the images are fully saturated and read the 41 Figure 28: Photograph of experimental setup including microchannel, microscope, and syringe pump. 42 highest intensity value of 255. Below 0.25 g/L is the linear regime, thus the solution should stay below this concentration in the video test section. This may be confirmed by ensuring all data from the experiments have intensity values below the corresponding 179 intensity units. A best fit curve is used to extract the slope and offset relating intensity and concentration within the linear regime, shown in Figure 29 and expressed by Intensity Value = 687.5C + 11.83 Figure 29: Calibration plot illustrating the relationship between intensity and concentration. Experimental results must be taken within the linear regime below 0.25 g/L for accurate conversion. 4.3.2 Results Experimental trials are completed for two microchannels of different aspect ratios, (a/b) = 0.22 ± 0.01 and (a/b) = 0.48 ± 0.03. To begin, the aspect ratio of (a/b) = 0.22 ± 0.01 is considered. For this trial, average flow velocity is set to 0.001 m/s, which corresponds to an infusion flow rate of 2779 nl/min on the syringe pump and Re=0.06. The entrance length of the flow conditions is 0.5 µm from the laminar flow inlet. 43 (a) Trial 1: (a/b) = 0.22 ± 0.01 (b) Trial 2: (a/b) = 0.48 ± 0.03 Figure 30: Microscope images of the two experimental microchannels. Channel (a) is 458.6 ± 1.0 µm wide. Channel (b) is 209.4 ± 0.5 µm in width. Both channels have a height of 101 ± 3 µm. It is vital for the experimental results that data is collected outside of the diffusive time scale. For rectangular channels, the dispersivity is known to reach 95% of its asymptotic value at 1 b2 ta = , 5κ so this limit will be used to ensure the solute is well outside of transitional times.2 This microchannel has geometry defined by a = 50.5 ± 1.5 µm and b = 229.3 ± 0.5 µm. Consequently, ta = 18.45 seconds and thus videos captured farther than 1.845 cm downstream from the solute entrance location will be safely within the long time regime. The Peclet number is 88.6 for this trial. The profile of the advecting solute as it passes through the test region is displayed in Figure 31, with videos taken 2 cm downstream from the solute inlet. The plots are generated by taking an intensity average over a cross sectional slice in the channel 100 pixels wide, which minimizes noise. The profiles are consistent with Taylor Dispersion literature and bode well for the functionality of craft-cutter fabricated microchannels. 44 Figure 31: Experimental solute profiles over the channel height at three different times for rectangular channel with (a/b) = 0.22 ± 0.01, Pe=88.6, with videos taken 2 cm downstream from solute inlet. Experimental frames included for each time step, with brightness doubled for clarity. To concretely quantify the success of the microchannel, the effective dispersion coefficient is calculated and compared to the theoretical prediction from Section 2.1. The benchmarking process begins by gener- ating a plot of fluorescein intensity as a function of time (Figure 33) on data collected past the diffusive time scale. For these trials, videos are captured 6 cm downstream from the solute inlet. Since the solute concentration transversing the test section is linearly related to intensity in these trials, either parameter may be used for the duration of analysis. Intensity will be used here. As can be seen in Figure 32, the curve of intensity versus time depicts an asymmetric trend. Videos are 45 collected at a fixed location along the channel. Time elapses as solute passes through the frame, allowing solute at the trailing edge more time to diffuse. This results in an asymmetrical distribution for intensity versus time. The experimental data points from this plot will be used to extract the experimental effective dispersion coefficient by following methods similar to Bontoux et al.16 Figure 32: Experimental data showing intensity over time for rectangular channel with aspect ratio of 0.22±0.01 with asymmetric Gaussian fitted curve (Equation (36). For this trial, Pe=88.6. Videos are taken 2 cm downstream from the solute inlet. The dispersing solute profile should reach a state in the long time regime shaped like a Gaussian centered around position x = Ut in the form −(x −Ut)2   C0 C(x,t) = p exp 4πKexpt 4Kexpt where C0 is an initial concentration factor and Kexp is the experimental effective dispersion coefficient. A variable tmax denotes the time of peak intensity at the location of interest. At one fixed position, x = Utmax , the expression above becomes −U 2 (tmax − t)2   C0 C(t) = p exp (35) 4πKexpt 4Kexpt 46 This distribution is asymmetric consistent with the experimental measurements. A fit is determined for the experimental data with the internal curve fitting application in MATLAB. A custom equation with four coefficients is input for the form of the fitted curve; −(h2 − t)2   h3 C(t) = √ exp + h4 (36) t h1 t where h1 = 4Kexp /U 2 , h2 = tmax , h3 = C0 / p 2πKexp , and h4 is a correction for nonzero minimum intensity values. Working in SI units, the fitted curve determined by MATLAB, shown in Figure 32, has coefficients of: h1 = 1.927 ± 0.136, h2 = 29.39 ± 0.08, h3 = 310.7 ± 9.2, and h4 = 8.15 ± 1.75. These may may be used to evaluate the experimental values of Kexp non-dimensionalized by κ, and fexp . The effective dispersion coefficient Kexp is determined by h1 , resulting in Kexp = 845.6 ± 117.5 κ with full error analysis included in Appendix A.3. The error analysis concludes that the parameter with the highest source of error is measurement of the channel height, followed by the error from confidence bounds of the curve fit parameters. Error from the other parameters are low enough to be neglected. By Equation (8), fexp is determined to be fexp = 5.65 ± 0.78 . The expected theoretical values are: Ktheory = 841 κ ftheory = 5.62 . for the appropriate properties and channel geometry.2 These results, plotted on Figure 34, show excellent agreement, providing evidence that craft-cutter microchannels fabricated by an optimized technique are capable of performing accurately for Taylor Dispersion type studies. To further substantiate this argument, another microchannel with a different aspect ratio is considered. The next trial with a channel of aspect ratio (a/b) = 0.48 ± 0.03 has geometry defined by a = 50.5 ± 1.5 47 µm, b = 104.7 ± 0.2 µm, average inflow velocity U = 0.006 m/s, Pe= 531.58, Re=0.34, and ta = 3.85 s. The syringe pump is set to 7613.8 nanoliters per minute, and videos are collected at a location 6 cm downstream from the solute inlet. The entrance length for flow conditions is 2 µm from the laminar flow inlet. By an identical process to the previous trial, a plot is generated of solute intensity over time and an asymmetric Gaussian fit is determined (Figure 33) using MATLAB’s curve fitting application. Figure 33: Experimental data intensity plot at times t > td overlaid with asymmetrical Gaussian fit curve of form Equation (36). For this trial, (a/b) = 0.48 ± 0.03, Pe=531.58, ta =3.85s and videos are taken 6 cm downstream from the solute inlet. The coefficients for this trial are determined to be: h1 = 1.12 ± 0.03, h2 = 15.23 ± 0.01, h3 = 53.03 ± 0.47, and h4 = 9.48±0.12, which results in the following values for effective dispersion coefficient and dispersion factor: Kexp = 17543.9 ± 2175.4 , κ fexp = 3.29 ± 0.40 . Again, methodology by Section 2.2 is used to evaluate the theoretical value of the dispersion coefficient and dispersion factor for a channel with the appropriate aspect ratio. The error analysis included in Appendix A.3 shows that error on the theoretical values for Kexp may be neglected for both trials. 48 Figure 34: Experimental values of the dispersion factor, f , overlaid with the theoretical predictions derived by Section 2.2. The experimental values agree with the theory within the experimental error. Complete error analysis is covered in Appendix A.3. Ktheory = 18820.8 κ ftheory = 3.50 Summarizing these results is the dispersion factor plot in Figure 34 49 5 Discussion Improved accessibility to research topics across all fields increases the likelihood of advances and dis- coveries. For a field such as microfluidics, ongoing research will directly improve chemical preparation, drug-delivery procedures and lab-on-a-chip devices, to name a few. Using the enclosed procedures, labs with no professional microfluidic equipment, low budgets, or no numerical simulation experience will now be able to study classical Taylor Dispersion problems and numerous variations. 50 References [1] G. I. Taylor. Dispersion of soluble matter in solvent flowing slowly through a tube. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 219(1137):186–203, 1953. [2] D. Dutta, A. Ramachandran, and D. T. Leighton. Effect of channel geometry on solute dispersion in pressure-driven microfluidic systems. Microfluid Nanofluid, 2:275–290, 2006. [3] M. Aminian, F. Bernardi, R. Camassa, D. M. Harris, and R. M. Mclaughlin. How boundaries shape chemical delivery in microfluidics. Science, 354(6317):1252–1256, 2016. [4] F. Bernardi. Space/Time Evolution in the Passive Tracer Problem. PhD thesis, 2018. [5] Po Ki Yuen and Vasiliy N. Goral. Low-cost rapid prototyping of flexible microfluidic devices using a desktop digital craft cutter. Lab Chip, 10(3):384–387, 2010. [6] Frank A. Gomez. The future of microfluidic point-of-care diagnostic devices. Future Science, Bio- analysis, (12):307, 2013. [7] Gina S. Fiorini and Daniel T. Chiu. Disposable microfluidic devices: fabrication, function, and appli- cation. BioTechniques, 38(3):429–446, 2005. [8] P. C. Chatwin and Paul J. Sullivan. The effect of aspect ratio on longitudinal diffusivity in rectangular channels. Journal of Fluid Mechanics, 120(-1):347, 1982. [9] R. Aris. On the dispersion of a solute in a fluid flowing through a tube. Proceedings of the Royal Society of London, (235):67–77, 1956. [10] Manuchehr Aminian, Francesca Bernardi, Roberto Camassa, and Richard M. Mclaughlin. Squaring the circle: Geometric skewness and symmetry breaking for passive scalar transport in ducts and pipes. Physical Review Letters, 115(15), 2015. [11] Qi Zhou and Chiu-On Ng. Electro-osmotic dispersion in a circular tube with slip-stick striped wall. Fluid Dynamics Research, 47(1), 2014. [12] Alessandra Adrover and Stefano Cerbelli. Laminar dispersion at low and high peclet numbers in finite- length patterned microtubes. Physics of Fluids, 29(6), 2017. [13] Particle tracing in a micromixer. COMSOL Application Gallery. 51 [14] Optimizing band dispersion in an electroosmotic flow through a curved microchannel. COMSOL Application Gallery. [15] Hao-Bing Liu and Hai-Qing Gong. Templateless prototyping of polydimethylsiloxane microfluidic structures using a pulsed co2 laser. Journal of Micromechanics and Microengineering, 19(3), 2009. [16] Nathalie Bontoux, Anne PÃl’pin, Yong Chen, Armand Ajdari, and Howard A. Stone. Experimental characterization of hydrodynamic dispersion in shallow microchannels. Lab Chip, 6(7):930–935, 2006. 52 A Appendix A.1 Two-Dimensional Velocity Conversion Beginning with Equations (24) & (25), the Laplacian of physical velocity, u, is ∂ 2u ∂p 1 = ∆u(y) = . ∂ y2 ∂x µ Aminian et al define their characteristic velocity,3 UP , a2 ∂ pP UP = . (37) µ ∂x They define the Laplacian of their velocity as 2 ∂ pP ∆uP (y, z) = . (38) µ ∂x Aminian et al is using a pressure-based velocity rather than the mean velocity.3 Relating the expressions, the constant is 3 UP = U. 2 A.2 Axisymmetric Velocity Conversion Beginning with the Laplacian of physical velocity,u, ∂ 2u 1 dp 2 = ∆u = − , dz 2µ dz the characteristic pressure velocity is a2 ∂ p p Up = . (39) µ ∂z The Laplacian of velocity is 2 ∂ pp ∆u p (y, z) = . (40) µ ∂z Considering the difference in Laplacian and velocity, the factor relating the two will be 53 Up = 2U. A.3 Error Analysis The error analysis begins by breaking down each final equation for Kexp and fexp into fundamental components that introduce error into the system, h1U 2 h1 Q2 N 2 Kexp = = 4 64a2 b2p Kexp (Q, h1 , a, b p , N) where Q is the flow rate with error introduced by the syringe pump, h1 is a fit coefficient with error from the curve fitting application within MATLAB, a is the channel half-height with measurement error using a micrometer, b is the channel half-width which has error from estimating the pixels across the channel (b p error), and error from measurements of the calibration slide itself converting from pixels to microns, denoted by N error. These are related so b = b p /N. As for f , 105 κ 2 κ 2 b2p   Kexp 105 κh1 f (κ, h1 , a, b p , Q) = −1 = − 840 κ 2 U 2 a2 8 a2 N 2 Q2 where each of these parameters has the same definition as before. The error analysis will proceed to find solutions to the following equations: ∂ Kexp ∂ Kexp Kexp (Q, h1 , a, b p , N) = Kexp (Qm , h1m , am , b pm , Nm ) + Qe + h1e + ... ∂Q ∂ h1 where any parameter with subscript m is the mean of that parameter and a subscript of e indicates the distribution of that parameter. The variation of a parameter is considered to be one single standard deviation of the measured values. This expression is assumed to be equivalent to Kexp (Q, h1 , a, b, N) = Kexp (Qm , h1m , am , bm , Nm ) ± ∆Kexp (41) where s 2  2 ∂ Kexp ∂ Kexp ∆Kexp = ∆Q + ∆h1 + ... (42) ∂Q ∂ h1 54 To fill in the necessary information, data is collected on the error in each parameter. For Q, the accuracy of the syringe pump is recorded within the apparatus documentation. The particular syringe pump being used herewithin is the Standard Infuse/Withdraw Pump 11 Elite Programmable Syringe Pump from Harvard Apparatus, which has an accuracy of ±0.5%. For the two trials, Qm and ∆Q may be calculated from the intended flow rate and appropriate geometry, resulting in Q1 and Q2 below for the lower aspect ratio channel and higher aspect ratio channel trials respectively. For the error on the channel height, ten measurements were taken at different locations along the channel with a micrometer. The resulting mean and standard deviation is that 2a = 0.101 ± 0.003018. For conversion of pixels to a unit of length by the calibration slide, ten measurements were taken on various images of the calibration slide from various marker points. This determined that Nm = 4.544 pixels per micron and ∆N = 0.011738 pixels per micron. As for b p , both the experimental channels were considered. First, the lower aspect ratio channel, to find b p,1 . Then the higher aspect ratio channel to find b p,2 , results below. For h1 , the 95% confidence bounds were recorded from the curve fitting application for both channels to determine h1,1 and h1,2 . ∆Q = ±0.5% Q1 = 4.63x10−11 ± 2.32x10−13 m3 /s Q2 = 1.27x10−10 ± 6.34x10−13 m3 /s a = 5.05x10−5 ± 3.02x10−6 m b p,1 = 1037.3 ± 2.38 pixels b p,2 = 473.7 ± 1.1 pixels N = 4.52 ± 0.012 pixels/µm 55 h1,1 = 1.93 ± 0.14 s h1,2 = 1.12 ± 0.028 s By evaluating Equations (42) & (41), the error analysis may be finalized for both fexp and Kexp values on each of the two experimental trials. This allows errorbars to be added to the plots comparing experimental results with theoretical predictions from Dutta et al.2 An error analysis on the aspect ratio (a/b) is also completed to generate horizontal error bars on plot. Code generated in Mathematica evaluates Equations (41) & (42). It is found that the greatest source of error in the experiment is from measurement of the channel height by use of a micrometer. The parameters with the second highest error is h, which is obtained from the confidence bounds on the fitted curve in MATLAB’s curve fitting application. For the lower aspect ratio channel, the error from Q, N, and b p may be neglected as they are many order of magnitude smaller than errors from h. The same terms may be neglected for the higher aspect ratio channel, as they are even smaller. The tape used for the experimental channels is manufactured to have a consistent height, so it is unlikely that the error is truly as high as the measurements suggest. For future trials, a more refined measurement technique should be used on the channel height to decrease error. A.4 Two-Dimensional Simulation Processing Code A.4.1 Defining Parameters and Importing Data % Defining constants, COMSOL parameters a = 5e-5; % [m] channel half-height kappa = 2.5e-10; % [m^2/s], diffusivity of material L = 0.05; % [m], length of channel U = .005; % [m/s], average inflow velocity x_m = 2.5e-4; % [m] starting location of initial solute patch sigma = 2e-5; % [m] standard deviation of initial gaussian plug Pe = (U*a/kappa)*(3/2); % Non-Dimensional Peclet number, converted to Brown from UNC by 3/2 tc = (a^2)/kappa; % non-dimensional time parameter, diffusive time scale t_steps = 200; % number of time steps (defined in comsol) 56 total_time = 5; % [s] grid_steps = 1000; % number of grid steps interpolated over channel length in export % Importing results from text file and formatting into matrices data = load(’gradient_mesh1.txt’); z_column = data(:,1); % defining steps along z-axis sweeps = 5; % number of parametric sweep iterations for t = 3:2+(t_steps+1)*sweeps % Populating concentration matrix with COMSOL data for i=1:grid_steps conc_exp(t-2,i) = data(i,t); end end t_exp = (0:total_time/t_steps:total_time)’; % defining time vector A.4.2 Theoretical Calculations %M2 Theoretical Calculation sumterm = 0; syms t for n=1:20 sumterm = sumterm + (exp(-t*(n*pi)^2))/((n*pi)^8); end M2 = @(t) (sigma^2)+(2*t)-((8*(Pe^2))/4725)+((16*t*(Pe^2))/945)+(16*(Pe^2)*sumterm); % M3 Theoretical Calculation sumterm = 0; syms t for n=1:20 sumterm = sumterm + ((exp(-t*(n*pi)^2))/((n*pi)^10))*(144-(24*t)-(1488/(n*pi)^2)); end M3 = @(t) ((1376*(Pe^3))/19348875)-((64*t*(Pe^3))/155925)+((Pe^3)*sumterm); 57 A.4.3 COMSOL Calculations % M2 COMSOL Calculation row = 1; for row = 1:t_steps+1 prefactor = 0; M2_val = 0; for x = 1:grid_steps M2_val = M2_val + ((z_column(x,1)-(x_m+(U*t_exp(row,1))))^2)*... conc_exp(row+((t_steps+1)*p),x); prefactor = prefactor + conc_exp(row+((t_steps+1)*p),x); end M2_exp5(row,1) = M2_val/prefactor; end % M3 COMSOL Calculations row = 1; for row = 1:t_steps+1 prefactor = 0; val = 0; for x = 1:grid_steps val = val + ((z_column(x,1)-(x_m+(U*t_exp(row,1))))^3)*conc_exp(row+((t_steps+1)*p),x); prefactor = prefactor + conc_exp(row+((t_steps+1)*p),x); end M3_exp5(row,1) = val/prefactor; end % Skewness COMSOL Calculations for sk_row = 1:(t_steps+1) sk_exp5(sk_row,1) = M3_exp5(sk_row,1)/((M2_exp5(sk_row,1))^(3/2)); end 58 A.4.4 Mesh Convergence Plotting figure; plot((t_exp(1:5:end)’/10), M3_exp1(1:5:end)/(a^3),’.-’) hold on plot((t_exp(1:5:end)’/10), M3_exp2(1:5:end)/(a^3),’.-’) hold on plot((t_exp(1:5:end)’/10), M3_exp3(1:5:end)/(a^3),’.-’) hold on plot((t_exp(1:5:end)’/10), M3_exp4(1:5:end)/(a^3),’.-’) hold on plot((t_exp(1:5:end)’/10), M3_exp5(1:5:end)/(a^3),’.-’) xlim([0 0.5]) ylim([-5e5 .5e5]) xlabel(’\tau (Non-Dimensional Time)’); ylabel(’M3’); title(’M3 Mesh Convergence Analysis’) legend(’N = 20’,’N = 30’,’N = 40’,’N = 50’,’N = 60’, ’Location’, ’southwest’) set(gca,’fontsize’, 16) % Zoom in on plot convergence axes(’position’,[.64 .63 .25 .25]) box on % put box around new pair of axes indexOfInterest = (t_exp/10 <= 0.5) & (t_exp/10 > 0.49); plot(t_exp(indexOfInterest),M3_exp1(indexOfInterest)) hold on plot(t_exp(indexOfInterest),M3_exp2(indexOfInterest)) hold on plot(t_exp(indexOfInterest),M3_exp3(indexOfInterest)) hold on plot(t_exp(indexOfInterest),M3_exp4(indexOfInterest)) 59 hold on plot(t_exp(indexOfInterest),M3_exp5(indexOfInterest)) axis tight set(gca,’ytick’,[]) A.4.5 Moment Benchmarking Plots clf close all % M2 Plots subplot(3,1,1); syms t fplot(M2(t),[0,.5],’k’,’LineWidth’, 1.5) % Theoretical M2 hold on plot((t_exp(1:5:end)’/10),M2_exp5(1:5:end)/(a^2),’o’) % Simulation M2 xlim([0 0.5]) ylim([-100 15000]) xlabel(’\tau (Non-Dimensional Time)’); ylabel(’M2’); legend(’Theory’,’COMSOL’, ’Location’, ’northwest’) % M3 Plots subplot(3,1,2); fplot(M3(t),[0,.5],’k’,’LineWidth’, 1.5) % Theoretical M3 hold on plot((t_exp(1:5:end)’/10), M3_exp5(1:5:end)/(a^3),’o’) % Simulation M3 xlim([0 0.5]) ylim([-5e5 .5e5]) xlabel(’\tau (Non-Dimensional Time)’); ylabel(’M3’); legend(’Theory’,’COMSOL’, ’Location’, ’southwest’) 60 % Skewness Plots subplot(3,1,3); fplot(M3(t)/(M2(t)^(3/2)),[0,.5],’k’,’LineWidth’, 1.5) % Theoretical skewness hold on plot((t_exp(1:5:end)’/10), sk_exp5(1:5:end),’o’) % Simulation skewness hold on xlim([0 0.5]) ylim([-.7 0]) xlabel(’\tau (Non-Dimensional Time)’); ylabel(’Skewness’); legend(’Theory’,’COMSOL’,’Location’, ’southeast’ ) A.5 Axisymmetric Simulation Processing Code A.5.1 Defining Parameters and Importing Data Defining constants, COMSOL parameters a = 5e-5; % m channel half-height kappa = 2.5e-10; % [m^2/s], diffusivity of material L = 0.05; % [m], length of channel U = .005; % [m/s], average inflow velocity zm = 5e-4; % starting location of initial dye patch sigma = 2e-5; % [m] standard deviation of initial gaussian plug Pe = (U*a/kappa)*(2); % Non-Dimensional Peclet number, converted to Brown from UNC by 3/2 tc = (a^2)/kappa; % non-dimensional time parameter t_steps = 50; % number of time steps (defined in comsol) total_time = 5; % [s] grid_steps = 1000; % number of steps interpolated over channel length in export % Importing results from text file and formatting into matrices data = load(’axi_gradient_mesh.txt’); % loading text file into matrix z_column = data(:,2); % first column of data matrix are our x values 61 % Populating concentration matrix with data values, each row is a time step p = 4; % Selection of parameter sweep iteration conc_axi = zeros(t_steps+1,grid_steps); for t = 3+(t_steps+1)*p:2+(t_steps+1)*p+(t_steps+1) % Populating concentration matrix for i=1:grid_steps conc_axi(t-2-(t_steps+1)*p,i) = data(i,t); end end t_exp = (0:total_time/t_steps:total_time)’; A.5.2 Theoretical Calculations sumterm = 0; syms t % M2 Theoretical Calculation for n = 1:20 bessel = besselzero(1,n); sumterm = sumterm + ((exp(-t*(bessel(end))^2))/((bessel(end))^8)); end M2 = @(t) ((sigma)^2)+ (2*t) + (Pe^2)*(((-1/1440)+(t/96))+(32*(sumterm))); % M3 Theoretical Calculation sumterm = 0; syms t for n = 1:20 bessel = besselzero(1,n); sumterm = sumterm + ((exp(-t*(bessel(end))^2))*((-240/(bessel(end))^12)+... (18/(bessel(end))^10)+(t/(bessel(end))^8))); end M3 = @(t) (Pe^3)*(((1/480)*((-17/896)+((112*t)/896))) + (16*(sumterm))); 62 A.5.3 Simulation Calculations % Calculating simulation M1 M1_axi_exp = zeros(t_steps+1,1); for row = 1:t_steps+1 prefactor = 0; val = 0; for x = 1:grid_steps val = val + ((z_column(x,1)-zm)*conc_axi(row,x)); prefactor = prefactor + conc_axi(row,x); end M1_axi_exp(row,1) = val/prefactor; end % Calculating simulation M2 M2_axi_exp = zeros(t_steps+1,1); row = 1; for row = 1:t_steps+1 prefactor = 0; val = 0; for x = 1:grid_steps val = val + ((z_column(x,1)-(zm+(U*t_exp(row,1))))^2)*conc_axi(row,x); prefactor = prefactor + conc_axi(row,x); end M2_axi_exp(row,1) = val/prefactor; end % Calculating simulation M3 M3_axi_exp = zeros(t_steps+1,1); row = 1; for row = 1:t_steps+1 prefactor = 0; 63 val = 0; for x = 1:grid_steps val = val + ((z_column(x,1)-(zm+(U*t_exp(row,1))))^3)*conc_axi(row,x); prefactor = prefactor + conc_axi(row,x); end M3_axi_exp(row,1) = val/prefactor; end % Calculating simulation skewness sk_axi = zeros(t_steps+1,1); row=1; for row = 1:t_steps+1 sk_axi(row,1) = M3_axi_exp(row,1)/((M2_axi_exp(row,1))^(3/2)); end A.5.4 Moment Benchmark Plotting clf close all % Plotting M2 subplot(3,1,1); fplot(M2(t),[0,0.5],’k’) % Theoretical M2 hold on plot((t_exp’/10),M2_axi_exp/(a^2), ’.r’) % Simulation M2 xlabel(’\tau (Non-Dimensional Time)’); ylabel(’M_2’); legend(’Mathematical Model’,’COMSOL’, ’Location’, ’northwest’) set(gca,’fontsize’, 14) % Plotting M3 subplot(3,1,2); fplot(M3(t),[0,0.5],’k’) % Theoretical M3 64 hold on plot((t_exp’/10),M3_axi_exp/(a^3), ’.r’) % Simulation M3 xlabel(’\tau (Non-Dimensional Time)’); ylabel(’M_3’); legend(’Mathematical Model’,’COMSOL’, ’Location’, ’northwest’) set(gca,’fontsize’, 14) % Plotting Skewness subplot(3,1,3); fplot(M3(t)/(M2(t)^(3/2)),[0,.5],’k’,’LineWidth’, 1.5) % Theoretical skewness hold on plot((t_exp’/10), sk_axi,’.r’) % Simulation skewness xlabel(’\tau (Non-Dimensional Time)’); ylabel(’Skewness’); legend(’Mathematical Model’,’COMSOL’,’Location’, ’southeast’ ) ylim([0 .4]) set(gca,’fontsize’, 14) A.5.5 Mesh Convergence Plotting clf close all fplot(M3(t),[0,0.5],’k’) % Theoretical M2 hold on plot((t_exp’/10),M3_axi_exp0/(a^3),’.-’) hold on plot((t_exp’/10),M3_axi_exp1/(a^3),’.-’) hold on plot((t_exp’/10),M3_axi_exp2/(a^3),’.-’) hold on plot((t_exp’/10),M3_axi_exp3/(a^3),’.-’) hold on 65 plot((t_exp’/10),M3_axi_exp4/(a^3),’.-’) xlabel(’\tau (Non-Dimensional Time)’); ylabel(’M3’); title(’Axisymmetric Mesh Convergence Analysis’) legend(’N = 10’,’N = 20’,’N = 30’, ’N = 40’,’N = 50’,’Location’,’Northwest’) set(gca,’fontsize’, 14) % Zooming in on convergence axes(’position’,[.63 .16 .25 .25]) box on % put box around new pair of axes indexOfInterest = (t_exp/10 <= 0.5) & (t_exp/10 > 0.485); % Range to zoom in on plot(t_exp(indexOfInterest),M3_axi_exp0(indexOfInterest)) hold on plot(t_exp(indexOfInterest),M3_axi_exp1(indexOfInterest)) hold on plot(t_exp(indexOfInterest),M3_axi_exp2(indexOfInterest)) hold on plot(t_exp(indexOfInterest),M3_axi_exp3(indexOfInterest)) hold on plot(t_exp(indexOfInterest),M3_axi_exp4(indexOfInterest)) axis tight set(gca,’ytick’,[]) A.6 Ballpark Calculations for 3D Experiments and Simulations a = 2.6e-5; % [meters] half-height of channel W = 21e-5; % [meters] full width of channel kappa = 5.7e-10; % [meters^2/s] diffusivity of fluorescein U = .005; % [m/s] average inflow velocity flowrate = U*((2*a)*W*1000) % [liters/second] infusion rate td = (1/20)*(W^2)/kappa; % Diffusive time scale 66 fprintf(’T = %5.5f seconds \n\n’,td); % Print diffusive time scale Pe = (U*a)/kappa; % Peclet number Pe_term = (1/210)*Pe^2; % Peclet term from dutta, should have Pe_term>>1 fprintf(’Pe_term = %5.5f >> 1\n\n’,Pe_term) % Print Peclet term syringe_flow_rate = flowrate*1e9*60 % [nanoliters/minute] flowrate if flowrate > 4.16e-10 && flowrate < 4.16e-4 % check if syringe pump capable of flowrate fprintf(’Good Flowrate \n\n’) else fprintf(’Bad Flowrate \n\n’) end Channel_Length = 2 * (td)*U*100; % Estimate a sufficient channel length fprintf(’Channel Length = %2.3f cm \n\n’, Channel_Length) % Print channel length A.7 Three-Dimensional Simulation Processing Code A.7.1 Initializing Constants % Defining constants, COMSOL parameters a = 5e-4; % [meters] channel half-height kappa = 5.7e-10; % [m^2/s], diffusivity of material L = 5.44*2; %.04*10; % [meters], length of channel U = .0310078; %.005; % [m/s], average inflow velocity x_m = 0.005; % starting location of initial dye patch tc = (a^2)/kappa; % non-dimensional time parameter t_steps = 10; % number of time steps (defined in comsol) total_time = 200; % [seconds] t_exp = (0:total_time/t_steps:total_time)’; % Dutta Theoretical Calculation 67 f = 1.76; % Dispersion factor, based on geometry and aspect ratio K = kappa + f*(kappa/210)*(U*2*a/kappa)^2; syms t M2_dutta = @(t) 2*K*t; A.7.2 Simulation Calculations clf close all % Calculating M2 from simulations using M0 and M2 integral from COMSOL M0 = [0.0000 6.2666E-10 % Values from COMSOL table, volume integral evaluating M0 5.0000 6.2680E-10 10.000 6.2694E-10 15.000 6.2709E-10 20.000 6.2723E-10 25.000 6.2738E-10 30.000 6.2752E-10 35.000 6.2766E-10 40.000 6.2781E-10 45.000 6.2795E-10 50.000 6.2809E-10 55.000 6.2824E-10 60.000 6.2838E-10 65.000 6.2853E-10 70.000 6.2867E-10 75.000 6.2882E-10 80.000 6.2896E-10 85.000 6.2910E-10 90.000 6.2925E-10 95.000 6.2939E-10 100.00 6.2954E-10 105.00 6.2968E-10 68 110.00 6.2983E-10 115.00 6.2997E-10 120.00 6.3012E-10 125.00 6.3026E-10 130.00 6.3040E-10 135.00 6.3055E-10 140.00 6.3069E-10 145.00 6.3084E-10 150.00 6.3098E-10 155.00 6.3112E-10 160.00 6.3126E-10 165.00 6.3138E-10 170.00 6.3144E-10 175.00 6.3136E-10 180.00 6.3101E-10 185.00 6.3019E-10 190.00 6.2871E-10 195.00 6.2642E-10 200.00 6.2330E-10]; M2_int = [0.0000 6.2736E-16 % From COMSOL evaluation of M2 integrand 5.0000 5.4785E-12 10.000 2.0782E-11 15.000 4.4497E-11 20.000 7.5664E-11 25.000 1.1318E-10 30.000 1.5632E-10 35.000 2.0440E-10 40.000 2.5690E-10 45.000 3.1336E-10 50.000 3.7337E-10 55.000 4.3730E-10 69 60.000 5.0367E-10 65.000 5.7223E-10 70.000 6.4308E-10 75.000 7.1599E-10 80.000 7.9074E-10 85.000 8.6716E-10 90.000 9.4509E-10 95.000 1.0244E-9 100.00 1.1049E-9 105.00 1.1866E-9 110.00 1.2692E-9 115.00 1.3529E-9 120.00 1.4373E-9 125.00 1.5226E-9 130.00 1.6087E-9 135.00 1.6954E-9 140.00 1.7827E-9 145.00 1.8707E-9 150.00 1.9591E-9 155.00 2.0481E-9 160.00 2.1374E-9 165.00 2.2265E-9 170.00 2.3141E-9 175.00 2.3979E-9 180.00 2.4743E-9 185.00 2.5381E-9 190.00 2.5857E-9 195.00 2.6145E-9 200.00 2.6255E-9]; A.7.3 Plotting M2 Simulation vs. Theory syms t 70 fplot((M2_dutta(t)/kappa),[0.1 1],’--’) hold on t = linspace(0,200,41); M2 = M2_int(:,2:end)./M0(:,2:end); plot(t(1:37)/tc,M2(1:37)/a^2,’-o’) hold on xlabel(’\tau’) ylabel(’M_2’) xlim([0 .5]) legend(’Dutta Theory’,’COMSOL Simulation’,’Location’,’Northwest’) set(gca,’fontsize’, 16) ylim([0 2e7]) xlim([0 0.5]) A.7.4 Slope Plotting clear M2_slopes % Initializing Euler Forward figure h = 5; % step’s size N = 39; % number of steps tau = 0:h/tc: h*9/tc % Non dimensional time and time steps % Calculating Simulation M2 slope for n=1:N M2_slopes(n+1)= (M2(n+1)- M2(n))/(h); end % Calculating Dutta Slope (Theoretical) M2_slope_dutta(1)=1; for n=1:N M2_slope_dutta(n+1)= (M2_dutta((n+1)*h)-M2_dutta(n*h))/h; 71 tau(n+1)=n*h/tc; end plot(tau(2:34),M2_slope_dutta(2:34),’--’) hold on plot(tau(2:33),M2_slopes(1:end,2:33),’-o’) xlabel(’\tau’); ylabel(’M_2 Slope’); legend(’Theory’,’COMSOL Simulation’,’Location’,’southeast’) set(gca,’fontsize’, 16) A.8 Intensity to Concentration Calibration Code function datamats = intensity_check(file); % Example input: %data_vec=intensity_check(’file’); close all MOVIEfile=VideoReader(file); Pend=MOVIEfile.NumberOfFrames; I=read(MOVIEfile,1); [~,RECT] = imcrop(I); %crop close all datamats=[]; for n=1:Pend I=read(MOVIEfile,n); 72 I=imcrop(I,RECT); datamats=cat(3,datamats,I); end [i1,i2] = size(datamats(:,:,1)); close all for n=1:size(datamats,3) % Dealing with outliers data_vec(:,n) = double(reshape(datamats(:,:,n),[i1*i2,1])); outlier_check = isoutlier(data_vec(:,n)); data_vec = data_vec(:,n).*double(~outlier_check); end data_vec(data_vec==0) = []; Im=max(data_vec); avg = mean(data_vec); stdev = std(data_vec); fprintf(’avg = %4.4f \n\n’,avg) fprintf(’stdev = %4.4f \n\n’,stdev) end A.9 Experimental Processing Code % Defining experimental constants x_m = 0.01; % meters, location of initial solute from laminar inlet U = .0033333; % m/s, average inflow velocity of fluid kappa = 5.7e-10; % meters^2/s, diffusivity of fluorescein 73 a = 5e-5; % [meters], half-height of channel td = a^2/kappa; % diffusive time scale % file = ’Users/AbbyTaylor/Dropbox (Brown)/Research/Thesis/Experiments/Vids/12cm.mov’ MOVIEfile=VideoReader(file); Pend=MOVIEfile.NumberOfFrames; % I=read(MOVIEfile,200); % choose a time saturated with fluorescein for accurate cropping % [~,RECT] = imcrop(I); %crop % Cropping video file for all frames and forming datamats tensor full of intensities datamats=[]; for n=1:Pend I=read(MOVIEfile,n); I=imcrop(I,RECT); datamats=cat(3,datamats,I); end [i1,i2] = size(datamats(:,:,1)); %imshow(datamats(:,:,n)) to show image of nth time step A.9.1 Generating Concentration Matrix for i = 1:size(datamats,3) conc_exp(i,:) = mean(datamats(100:600,:,i)); % mean of domain to eliminate noise end % Concentration Plots x = linspace(0,210,size(conc_exp,2)); plot(x,conc_exp(500,:)) xlim([0 210]) xlabel(’y [um]’) 74 ylabel(’Fluorescein Intensity’) set(gca,’fontsize’, 15) A.9.2 Concentration vs Time, calculation and plot cmin = 10; clear I_avg I_avg = mean(conc_exp’)’-cmin; C_avg = I_avg/364.1; time = linspace(40,size(I_avg,1)/10,size(I_avg,1)); plot(time, C_avg) xlabel(’Time [s]’) ylabel(’Concentration [g/L]’) %xlim([35 55]) ylim([0 .07]) set(gca,’fontsize’, 15) A.9.3 M2 Calculations % Theoretical Calculation f = 3.9982; % dispersion factor, get from COMSOL PDE solving program Deff_dutta = (2*f/105)*(U^2*a^2/kappa); Deff_exp = (24.4629*U^2)/(4*48.42); A.9.4 Trial 1: Aspect Ratio (a/b=0.22) Evaluating f at the mean of each parameter κ = 5.7 ∗ 10−10 ; (*m∧ 2/s*) h = 1.927; (*s*) a = 5.05 ∗ 10−5 ; (*m*) bp = 2074.6/2; (*pixels*) n = 4.524 ∗ 106 ; (*pixels/meter*) Q = 4.6319 ∗ 10−11 ; (*m∧ 3/s*) 75 Remove[ f ]   κ 2 bp2 f = 105 κh 8 a2 − 840 n2 Q2 5.64624 Evaluating K at the mean of each parameter 2 2 hQ n K = 64a2 bp2 4.818077378078232`*∧ -7 Calculating ∆f Clear[“Global`*”] ∆h = 0.136; (*seconds*) ∆a = 3.018 ∗ 10−6 ; (*meters*) ∆bp = 2.3828; (*pixels*) ∆n = 0.011738 ∗ 106 ; (*pixels/meter*) ∆Q = 2.3159 ∗ 10−13 ; (*m∧ 3/s*)  2 2 κ bp f = 105 κh 8 a2 − 840 n2 Q2 ; p ∆f = (D[ f , h]∆h)2 + (D[ f , a]∆a)2 + (D[ f , bp]∆bp)2 + (D[ f , n]∆n)2 + (D[ f , Q]∆Q)2 ; ans1 = N (D[ f , h]∆h)2 ; (*calculating which parameter introduces most error*)   ans2 = N (D[ f , a]∆a)2 ;   ans3 = N (D[ f , bp]∆bp)2 ;   ans4 = N (D[ f , n]∆n)2 ;   ans5 = N (D[ f , Q]∆Q)2 ;   κ = 5.7 ∗ 10−10 ; (*m∧ 2/s*) h = 1.927; (*s*) a = 5.05 ∗ 10−5 ; (*m*) bp = 2074.6/2; (*pixels*) n = 4.524 ∗ 106 ; (*pixels/meter*) 76 Q = 4.6319 ∗ 10−11 ; (*m∧ 3/s*) N[∆f] N[ans1] N[ans2] N[ans3] N[ans4] N[ans5] 0.78466 0.15917 0.456522 9.440061880336055`*∧ -10 1.2043464364258696`*∧ -9 4.4722870461102055`*∧ -9 Calculating ∆K Clear[“Global`*”] ∆h = 0.136; (*seconds*) ∆a = 3.018 ∗ 10−6 ; (*meters*) ∆bp = 2.3828; (*pixels*) ∆n = 0.011738 ∗ 106 ; (*pixels/meter*) ∆Q = 2.3159 ∗ 10−13 ; (*m∧ 3/s*) hQ n 2 2 K = 64a2 bp2 ; p ∆K = (D[K, h]∆h)2 + (D[K, a]∆a)2 + (D[K, bp]∆bp)2 + (D[K, n]∆n)2 + (D[K, Q]∆Q)2 ; κ = 5.7 ∗ 10−10 ; (*m∧ 2/s*) h = 1.927; (*s*) a = 5.05 ∗ 10−5 ; (*m*) bp = 2074.6/2; (*pixels*) n = 4.524 ∗ 106 ; (*pixels/meter*) 77 Q = 4.6319 ∗ 10−11 ; (*m∧ 3/s*) N[∆K] 6.71342892830287`*∧ -8 Calculating Error on Aspect Ratio Clear[“Global`*”] ∆a = 3.018 ∗ 10−6 ; (*meters*) ∆bp = 2.3828; (*pixels*) ∆n = 0.011738 ∗ 106 ; (*pixels/meter*) an ar = bp ; p ∆ar = (D[ar, a]∆a)2 + (D[ar, bp]∆bp)2 + (D[ar, n]∆n)2 ; a = 5.05 ∗ 10−5 ; (*m*) bp = 2074.6/2; (*pixels*) n = 4.524 ∗ 106 ; (*pixels/meter*) N[∆ar] 0.0131846 Calculating Error on Theoretical K Clear[“Global`*”] κ = 5.7 ∗ 10−10 ; (*m∧ 2/s*) h = 1.927; (*s*) a = 5.05 ∗ 10−5 ; (*m*) bp = 2074.6/2; (*pixels*) n = 4.524 ∗ 106 ; (*pixels/meter*) Q = 4.6319 ∗ 10−11 ; (*m∧ 3/s*) F = 5.6178; Remove[ f ] 78  2 2 Qn Ktheory = 1 + 105 4bpκ F 841.025 Clear[“Global`*”] ∆bp = 2.3828; (*pixels*) ∆n = 0.011738 ∗ 106 ; (*pixels/meter*) ∆Q = 2.3159 ∗ 10−13 ; (*m∧ 3/s*) p ∆Ktheory = (D[Ktheory, bp]∆bp)2 + (D[Ktheory, n]∆n)2 + (D[Ktheory, Q]∆Q)2 ; F = 5.6178; κ = 5.7 ∗ 10−10 ; (*m∧ 2/s*) bp = 2074.6/2; (*pixels*) n = 4.524 ∗ 106 ; (*pixels/meter*) Q = 4.6319 ∗ 10−11 ; (*m∧ 3/s*) N[∆Ktheory] 0. A.9.5 Trial 2: Aspect Ratio (a/b=0.48) Evaluating f at the mean of each parameter Clear[“Global`*”] κ = 5.7 ∗ 10−10 ; (*m∧ 2/s*) h = 1.121; (*s*) a = 5.05 ∗ 10−5 ; (*m*) bp = 473.7; (*pixels*) n = 4.524 ∗ 106 ; (*pixels/meter*) Q = 1.269 ∗ 10−10 ; (*m∧ 3/s*) Remove[ f ]   κ 2 bp2 f = 105 κh 8 a2 − 840 n2 Q2 3.28831 79 Evaluating K at the mean of each parameter hQ n 2 2 K = 64a2 bp2 0.0000101039 Calculating ∆f Clear[“Global`*”] ∆h = 0.028; (*seconds*) ∆a = 3.018 ∗ 10−6 ; (*meters*) ∆bp = 1.107; (*pixels*) ∆n = 0.011738 ∗ 106 ; (*pixels/meter*) ∆Q = 6.3448 ∗ 10−13 ; (*m∧ 3/s*)  2 2 κ bp f = 105 κh 8 a2 − 840 n2 Q2 ; p ∆f = (D[ f , h]∆h)2 + (D[ f , a]∆a)2 + (D[ f , bp]∆bp)2 + (D[ f , n]∆n)2 + (D[ f , Q]∆Q)2 ; ans1 = N (D[ f , h]∆h)2 ; (*calculating which parameter introduces most error*)   ans2 = N (D[ f , a]∆a)2 ;   ans3 = N (D[ f , bp]∆bp)2 ;     ans4 = N (D[ f , n]∆n)2 ; ans5 = N (D[ f , Q]∆Q)2 ;   κ = 5.7 ∗ 10−10 ; (*m∧ 2/s*) h = 1.121; (*s*) a = 5.05 ∗ 10−5 ; (*m*) bp = 473.7; (*pixels*) n = 4.524 ∗ 106 ; (*pixels/meter*) Q = 1.269 ∗ 10−10 ; (*m∧ 3/s*) N[∆f] N[ans1] N[ans2] N[ans3] 80 N[ans4] N[ans5] 0.401547 0.00674682 0.154493 7.541940884179388`*∧ -13 9.296889508814442`*∧ -13 3.4522900869458783`*∧ -12 Calculating ∆K Clear[“Global`*”] ∆h = 0.028; (*seconds*) ∆a = 3.018 ∗ 10−6 ; (*meters*) ∆bp = 1.107; (*pixels*) ∆n = 0.011738 ∗ 106 ; (*pixels/meter*) ∆Q = 6.3448 ∗ 10−13 ; (*m∧ 3/s*) hQ n 2 2 K = 64a2 bp2 ; p ∆K = (D[K, h]∆h)2 + (D[K, a]∆a)2 + (D[K, bp]∆bp)2 + (D[K, n]∆n)2 + (D[K, Q]∆Q)2 ; κ = 5.7 ∗ 10−10 ; (*m∧ 2/s*) h = 1.121; (*s*) a = 5.05 ∗ 10−5 ; (*m*) bp = 473.7; (*pixels*) n = 4.524 ∗ 106 ; (*pixels/meter*) Q = 1.269 ∗ 10−10 ; (*m∧ 3/s*) N[∆K] 1.237941141726084`*∧ -6 81 Calculating Error on Aspect Ratio Clear[“Global`*”] ∆a = 3.018 ∗ 10−6 ; (*meters*) ∆bp = 1.107; (*pixels*) ∆n = 0.011738 ∗ 106 ; (*pixels/meter*) ar = an bp ; p ∆ar = (D[ar, a]∆a)2 + (D[ar, bp]∆bp)2 + (D[ar, n]∆n)2 ; a = 5.05 ∗ 10−5 ; (*m*) bp = bp = 473.7; (*pixels*) n = 4.524 ∗ 106 ; (*pixels/meter*) N[∆ar] 0.0288721 Calculating Error on Theoretical K Clear[“Global`*”] κ = 5.7 ∗ 10−10 ; (*m∧ 2/s*) bp = 473.7; (*pixels*) n = 4.524 ∗ 106 ; (*pixels/meter*) Q = 1.269 ∗ 10−10 ; (*m∧ 3/s*) F = 3.4969; Remove[ f ]  2 2 Qn Ktheory = 1 + 105 4bpκ F 18820.8 Clear[“Global`*”] ∆bp = 2.3828; (*pixels*) ∆n = 0.011738 ∗ 106 ; (*pixels/meter*) ∆Q = 2.3159 ∗ 10−13 ; (*m∧ 3/s*) p ∆Ktheory = (D[Ktheory, bp]∆bp)2 + (D[Ktheory, n]∆n)2 + (D[Ktheory, Q]∆Q)2 ; 82 F = 5.6178; κ = 5.7 ∗ 10−10 ; (*m∧ 2/s*) bp = 2074.6/2; (*pixels*) n = 4.524 ∗ 106 ; (*pixels/meter*) Q = 4.6319 ∗ 10−11 ; (*m∧ 3/s*) N[∆Ktheory] 0. 83