Abstract: The dynamics of passive scalar advection has been modeled with a set of partial differential equations. The study of the passive scalar advection in a micro fluidic device has been observed. Numerical simulation of the scalar advection has been established from the model equation. The advection of the diluted fluid inside the micro-device occurred due to the pressure gradient of the fluid. A set of nonlinear partial differential equations has been solved numerically for the advection process. Splitting based technique has been developed to simulate the passively advected scalar field. The technique has been observed that saves computational time, requires less computer memory and workable with large Peclet number.
INTRODUCTION
Transfer of fluid with scalar quantity known as passive scalar advection. A mathematical model of passive scalar advection in a micro fluidic device has been presented. Hydrodynamically driven flow has been focused in a micro chip. The width of a typical chip varies from 10 to 100 μm (Seiler et al., 1994; Crabtree et al., 2001). This model describes the flow of fluid with a dissolved scalar that consists of a set of elliptic, mixed-parabolic-hyperbolic partial differential equations. The numerical treatment of the model is based on Splitting Based Technique (Alam and Bowman, 2002a), decoupled the pressure field from the momentum equation by means of a Poisson equation.
MODEL EQUATION
Assumption: It is assumed that the solution inside the buffer of the microchip is incompressible and Newtonian with uniform physical properties. The solution is a mixture of a substance called scalar and a dilute fluid. The concentration of the second substance is defined by the amount of the substances per unit volume which is a passive scalar quantity. It is further assumed that the density of the scalar does not change even if the concentration of the fluid changes. It is mentioned that the pressure of the scalar does not change because the fluid is Newtonian.
Thus the mass inside a typical micro chip is conserved and can be expressed in the form of a partial differential equation:
| (1) |
where,
| (2) |
From the above two equation we get:
| (3) |
which is known as incompressibility constraint. It is to be noted that Eq. 3 is independent of density.
Model equations: The assumptions have been made to model the dynamics of the fluid scalar advection in micro fluidic devices. Following the assumptions it can be said that fluid inside the micro device satisfies the equation:
| (4) |
where, P(x, t) is the pressure of the fluid, v is the Kinematic viscosity of the fluid and f(x, t) is external force.
The concentration c(x, t) of the dissolved samples satisfies the Fick’s law, i.e.,
where, Γ(c) = 0 is the flux of the c and we can define
such that D is the molecular diffusion coefficient, Γf flux due to the external force. Since, there is no external force in the model, we assume that Γf = 0. So, the concentration equation becomes:
| (5) |
The Eq. 3, 4 and 5 describe a model of passive scalar transport in a micro fluidic device.
NON DIMENSIONAL MODEL
The non-dimensional model: Using dimensionless quantities, in the model equations, the new non-dimensional forms of the model equations are obtained which can be written in component form:
| (6) |
| (7) |
| (8) |
where, Re is the Reynolds number and Pe is the Peclet number. These non-dimensional numbers play a significant in numerical simulations. If Re<<1 the flow is dominated by advection or inertia force and if Re>>l, then the flow is dominated by viscous force. Similar arguments hold for Peclet number Pe.
PRESSURE POISSON FORMULATION
The Poisson equation for the pressure field has been derived by taking divergence on Eq. 7-8, and applying the incompressibility constraint Eq. 6.
| (9) |
| (10) |
The governing model equations can be written in flux-conservative form excluding incompressibility constraint.
| (11) |
| (12) |
| (13) |
where,
| (14) |
where:
Eq. 13 and 14 are described as the pressure-Poisson formulation of two-dimensional passive scalar advection model.
NUMERICAL SOLUTION FOR SIMULATIONS
Numerical simulation parameters: The dimensional and the non-dimensional parameters that have been used in the numerical treatment are Length scale h = 5x10-4 m, Velocity scale U = 10-3 m, Diffusion coefficient D = 10-10 m2 sec-1, Kinematic viscosity v = 10-10 m2 sec-1, Reynolds number Re = 0.1, Maximum concentration Cmax = 10-6 mol, Peclet number Pe = 1000.
The non-dimensional parameters have been evaluated from characteristic length and velocity scales based on typical laboratory experiments.
Computational domain: Simulation of passive scalar advection in a micro-fluidic device requires pumping of fluid through a network of micro-chips (Seiler et al., 1994). One of the widely used instruments is the cross channel device as shown Fig. 1. The experimental micro fluidic device has an input boundary and an output boundary and two wall boundaries.
Boundary condition: The boundary conditions for the velocity and the concentration field are Dirichlet on the channel walls and Neumann on the input and output open boundaries. The Dirichlet boundary conditions for the pressure field on the input and output boundaries impose a pressure gradient along the channel and Neumann boundary conditions are used on the walls. Specific boundary conditions are needed that no flux of physical quantity is coming in or going out through the boundaries.
Numerical discretization: The numerical grid has been demonstrated in Fig. 2. A continuous function F(x, y) is mapped to approximate the grid cell as F(x, y)→F(ihx, jhy) = Fi,j, were ihx and jhy are the distances between the centroids of two consecutive cells in x-and y-directions, respectively,
| Fig. 1: | A schematic diagram of a cross channel device with the portion of the injection channel that is considered in the simulation as computational domain |
| Fig. 2: | The numerical grid, describing the ghost cell and grid cell |
For convenience, the semi-discrete set of equation can be written as:
| (15) |
| (16) |
where, Ui,j is the approximation of U at (i,j) cell and ∇1, ∇2 and ∇2 are discretized divergence (and hence gradient) operator, compact Laplacian and non-compact Laplacian operators respectively. These discrete operators are defined as:
| (17) |
| (18) |
| (19) |
The numerical discretization of the Navier-Stokes equation has become complicated by the presence of a pressure-gradient term on the right hand side. A compact discretization has been used for the viscous and diffusion term but a non-compact discretization for the pressure-Poisson equation. This can be verified by simply replacing Fi,j with Fi,j in the Eq. 17. The numerical solution obtained from this scheme should automatically satisfy the incompressibility constraint Eq. 3. In a fluid simulation, care must be taken in order to obtain an incompressible velocity field at each time step (Freziger and Peric, 2002).
Solving procedure
Step 1: Compute the pressure field form Eq. 16,
the equation has been designed in such a way that the computed pressure field
preserves the solenoidal nature of the velocity field that satisfies the Eq.
3 (Alam and Bowman, 2002b).
Step 2: Solve the pure advection initial value problem with known initial conditions by Euler explicit method.
| (20) |
Step 3: Solve the diffusion problem
| (21) |
where, Ui,j* is the solution found in step 2.
One of the significant points that has to be noted is that the calculation of the pressure field should satisfies the solenoidal nature of the velocity field. In the SIMPLE method of Patankar and Hu (1998), the velocity and the pressure have to be corrected at each time step until the incompressible velocity field is obtained.
By applying Euler explicit method for step 2, this numerical treatment can be accomplished as
| (22) |
where:
| (23) |
This consistency of the scheme defined by Eq. 22 can be checked easily by substituting Eq. 23 into Eq. 22, which gives:
| (24) |
Eq. 24 is the first order temporal approximation of the semi-discrete Eq. 15. This proves that the Eq. 22 is consistent (Alam and Bowman, 2002b).
Numerical code: The Splitting Based Technique has been used for the Numerical Simulation of Passive Scalar Advection in Micro Fluidic Device. The SBT solver preserves the solenoidal nature of the velocity field. The numerical result of the flow characteristics has compared with analytical and experimental results. The numerical results found to be in good qualitative agreement with theoretical and experimental observations.
RESULTS AND DISCUSSION
The principal aim of this study is to develop a numerical scheme that is capable of simulating a passively advected scalar field. Splitting based technique in short SBT solver has been developed, tested and found the following results. SBT Solver save the computational time requires less computer memory and is Convergent. SBT Solver is found to work with large Peclet number Pe≥1. At the point artificial damping the profile generated by SBT solver, has higher amplitudes at the same point of time as that of upwind solver and takes longer time to be damped. The velocity field takes a parabolic shape which is in good qualitative agreement with the realistic situation. The SBT solver developed is based on the time splitting of the viscous term from the inertia term from the model equations. Once the pressure field is computed by the pressure solver, a guess velocity field is obtained with any suitable numerical treatment of the nonlinear terms. This guess velocity field is used to integrate the viscous part with a trapezoidal method. A trapezoidal integration then required to solve a Helmoholtz type linear system of equations. Thus a Helmoholtz solver was developed based on a Gauss-Seidel relaxation scheme. Since, the previously computed velocity field is used an initial guess, an enormous amount of computational time and memory is saved. In comparison with the Helmholtz like linear system the SBT solver needs less computational time and requires less computer memory. In order to test the convergence of the SBT solver, the norm of the x-component of the velocity field was computed, which is shown in Fig. 7.
The numerical simulation of passive scalar advection is one of the most complicated and computationally expensive tasks that remain as an unsolved puzzle in computational fluid dynamics. An Eulerian numerical scheme suffers from numerical instability as the non-dimensional Peclet number increases; i.e., Pe≥1. To stabilize the result one may need to introduce artificial diffusion. Artificial diffusion dampens the original solution and such a solution is not very useful. For a passive scalar advection problem like capillary electrophoresis or DNA sequencing, Eulerian schemes do not provide qualitative result and do not preserve the boundedness and conservative properties. At this point SBT Solver has been found workable with large Pe. When Pe≥1 the nonlinear term dominates and causes the transfer of the scalar field from one place to another place. Thus an initial pulse of a scalar quantity travels with the flow speed. The scalar field c is distributed in the computational region so that the non-dimensional value of c becomes zero every where except in a small strip near the input boundary where it becomes 1. It is shown in Fig. 3, white is 1 and black is zero. The initial distribution of c is shown in Fig. 4. Numerical experiments are performed to examine the distribution of the scalar field at a later time. The distribution of c is than computed with the SBT solver at three locations (called location 1, location 2 and location 3) on the center line of the channel. The first location is at a quarter ways upstream from the input boundary, the second is at a distance of half of the channel length in upstream and the third is at a distance of three quarter of the channel length in upstream. The profiles of c at these three different locations are computed and the results are presented in Fig. 5. The prediction of the pulse of a traveling c is thus found wiggless-free although a second order-centered-in-space discretization is being used for advection derivatives. The presence of wiggles in the numerical solution of an advection-diffusion equation is because of the numerically accumulated energy in the system. Thus many researchers have been trying to develop elegant and accurate algorithms to obtain wiggless solution.
| Fig. 3: | Initial distribution of scalar field |
| Fig. 4: | The initial profile of the scalar field |
| Fig. 5: | Distribution of the traveling scalar pulse at three different locations computed by the SBT solver |
The first solution to this problem is the method of Neumann and Richtmyer (1950) who introduced numerical damping. Artificial diffusion dampens the numerically accumulated energy which helps the disappearance of the wiggles, but the inherent computational damping destroys the original profile from the computation. At this point, the profile of C given in Fig. 4 and 5 show that the initial traveling square wave is distorted; the sharp edges are smoothed over and takes a shape of a Gaussian distribution. To view a clear picture of non-realistic damping of the artificial viscosity method, the profile c at location 3, computed by a traditional upwind method one compared with that produced by SBT solver. The distribution of c in Fig. 6 shows the artificial damping of the upwind method in comparison with the SBT solver.
The two dimensional flow in a rectangular computational domain with no-slip boundary conditions in the x-direction and Neumann boundary conditions in the other directions leads to a parallel share flow. A pressure drop along the channel causes a fluid flow in the direction of channel length. With a moderate Reynolds number, the x-component of the velocity field takes a parabolic shape once the flow field is fully developed. A series of numerical experiments have been invoked to examine the performance of the SBT solver by computing the x-component of the velocity field with different computational domain and different grid. This solver has the capability of computing the velocity field with any computational domain of rectangular size and with any numerical grid. The x-component of the velocity field computed at the middle of the channel along a line parallel to the y-axes is shown in the Fig. 8. The parabolic profile is in good qualitative agreement with the realistic situation. In order to show the distribution of this velocity field, a snapshot of the field at a typical time step when the velocity is fully developed is produced. Using a color spectrum with black for the lowest value and white for the highest value, this computation is presented in Fig. 9.
| Fig. 6: | Distribution of the traveling scalar pulse computed at location 3 with SBT solver and a first upwind solver |
| Fig. 7: | The norm of the velocity error |
| Fig. 8: | x-component of the velocity field |
| Fig. 9: | x-component of the velocity field |
CONCLUSION
A mathematical model of passive scalar transport in a microchip has been presented. The model describes the flow of fluid carrying a second substance and consists of a set of elliptic, mixed-parabolic-hyperbolic partial differential equations. The equation describing the concentration of the second substance is dominated by the advective terms and hence is of hyperbolic nature. The formulation of the model decouples the pressure field from the momentum equation by means of a Poisson equation. The accumulation of the numerical energy is treated so that physical viscosity can be dissipating the extra energy in the system.
ACKNOWLEDGMENT
The authors are grateful to Prof. Dr. Dilder Hossain, for valuable suggestions.