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, mixedparabolichyperbolic 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:
where,
is the velocity of the fluid and ρ(x, t) is the density of the fluid. The
incompressibility of the fluid means that the density does not change with the
change of pressure.
From the above two equation we get:
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:
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:
The Eq. 3, 4 and 5 describe
a model of passive scalar transport in a micro fluidic device.
NON DIMENSIONAL MODEL
The nondimensional model: Using dimensionless quantities, in the model equations, the new nondimensional forms of the model equations are obtained which can be written in component form:
where, R_{e} is the Reynolds number and P_{e} is the Peclet number. These nondimensional numbers play a significant in numerical simulations. If R_{e}<<1 the flow is dominated by advection or inertia force and if R_{e}>>l, then the flow is dominated by viscous force. Similar arguments hold for Peclet number P_{e. }
PRESSURE POISSON FORMULATION
The Poisson equation for the pressure field has been derived by taking divergence
on Eq. 78, and applying the incompressibility
constraint Eq. 6.
The governing model equations can be written in fluxconservative form excluding incompressibility constraint.
where,
is given by Eq. 9. Thus, Eq. 10 and 11
can be rewritten as below:
where:
Eq. 13 and 14 are described as the pressurePoisson
formulation of twodimensional passive scalar advection model.
NUMERICAL SOLUTION FOR SIMULATIONS
Numerical simulation parameters: The dimensional and the nondimensional 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} m^{2} sec^{1}, Kinematic viscosity v = 10^{10} m^{2} sec^{1}, Reynolds number R_{e} = 0.1, Maximum concentration C_{max} = 10^{6} mol, Peclet number P_{e} = 1000.
The nondimensional parameters have been evaluated from characteristic length and velocity scales based on typical laboratory experiments.
Computational domain: Simulation of passive scalar advection in a microfluidic
device requires pumping of fluid through a network of microchips (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(ih_{x}, jh_{y}) = F_{i,j}, were ih_{x} and jh_{y }are the distances between the centroids of two consecutive cells in xand ydirections, 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 semidiscrete set of equation can be written as:
where, U_{i,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 noncompact Laplacian operators respectively. These discrete operators are defined as:
The numerical discretization of the NavierStokes equation has become complicated
by the presence of a pressuregradient term on the right hand side. A compact
discretization has been used for the viscous and diffusion term but a noncompact
discretization for the pressurePoisson equation. This can be verified by simply
replacing F_{i,j} with F_{i,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.
Step 3: Solve the diffusion problem
where, U_{i,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
where:
This consistency of the scheme defined by Eq. 22 can be checked easily by substituting Eq. 23 into Eq. 22, which gives:
Eq. 24 is the first order temporal approximation of the
semidiscrete 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 P_{e}≥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 GaussSeidel 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 xcomponent 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 nondimensional Peclet number increases; i.e., P_{e}≥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 P_{e}. When P_{e}≥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 nondimensional 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 wigglessfree
although a second ordercenteredinspace discretization is being used for advection
derivatives. The presence of wiggles in the numerical solution of an advectiondiffusion
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 nonrealistic 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 noslip boundary conditions in the xdirection 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 xcomponent 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 xcomponent 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 xcomponent of the velocity field computed at the middle of the channel along a line parallel to the yaxes 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: 
xcomponent of the velocity field 

Fig. 9: 
xcomponent 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, mixedparabolichyperbolic 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.