Search in book...
Toggle Font Controls
Create new playlist

Name your new playlist

Playlist description (optional)
Sign In

Email address

Password

Forgot Password?

or

Continue with Facebook

Continue with Google
Sign Up

Full Name

Email address

Confirm Email Address

Password

or

Continue with Facebook

Continue with Google

Chapter 19

Taylor-Aris Dispersion

19.1 Introduction

As we have learned, most microfluidic flows are strictly laminar. For many applications this is an advantage, as the flow can be controlled very precisely. However, applications, such as, mixing are difficult in microfluidics because mass transport is mostly based on molecular diffusion. If two laminar flows are co-flowing, the mass transport across the interface is purely based on diffusion. However, if a flow is segmented, i.e., a plug of one liquid is (repeatedly) inserted into a second stream, the boundary is normal to the flow direction and thus experiences convective flow. In these cases we find that mixing is significantly enhanced. This is due to the fact that the molecular mass transport by diffusion is supported by convective mass transport. In general, we refer to the transport of mass (irrespective of whether it is due to convection or diffusion) as dispersion. The most important model for dispersion that we come across in microfluidics is the Taylor-Aris dispersion that was first introduced by Taylor^¹ in 1953 [1] and later refined by Aris^² [2].

19.2 Dispersion

19.2.1 Dispersion on a Static Fluid Plug

If a static plug of a fluid in inserted into a second fluid, it will diffuse over time. This happens also in a microfluidic channel where such a liquid plug may be inserted, e.g., via channel intersections (see section 19.9). Due to diffusion, the sharp plug interfaces, i.e., the sharp concentration gradients will smear out (see Fig. 19.1).

f19-01-9781455731411 — Fig. 19.1 Static plug of a fluid or a solute inserted in a second fluid in microfluidic channel. Over time the plug will slowly disperse into the surrounding fluid due to diffusion.

19.2.2 Dispersion in a Moving Fluid Plug

If the bulk of the liquid is moving, e.g., due to an applied pressure drop across the microfluidic channel, the interface will be curved due to the characteristic flow profiles that will form. This leads to strong concentration gradients and therefore to oriented diffusion (see Fig. 19.2a). In the case of diffusion only along the axis of the channel, the curved interfaces will be continually expanding, keeping their characteristic shape (see Fig. 19.2b).

f19-02-9781455731411 — Fig. 19.2 Moving plug of a fluid or a solute inserted in a second fluid in a microfluidic channel. a) Due to bulk liquid movement, e.g., as a consequence of an applied pressure drop, the plug will move forming a characteristic flow profile. This creates concentration gradients along the curved interfaces. b) If the substance of which the plug consists shows strong axial diffusion, the shape of the plug will be smeared into a concentration profile along the length of the channel. c) If axial diffusion is strong the plug profile will straighten out eventually forming a near-rectangular shape again.

If only radial diffusion is considered, the curved surfaces will blur and the plug will eventually obtain the shape of a rectangular pulse traveling along the channel (see Fig. 19.2c). In practical applications, there will always be a mixture of both cases as diffusion can mostly be considered as being isotropic. However, due to the significantly shorter length scales in the radial compared to the axial dimension, radial diffusion happens significantly faster.

In the following we will discuss the case of convection-diffusion in cylindrical cross-sections. We have already studied the case of laminar flows in circular cross-sections that is referred to as Hagen-Poiseuille flow in section 16.3 where we have found parabolic flow profiles described by Eq. 16.36

$v (r) = 2 v_{a v} (1 - {(\frac{r}{R})}^{2})$ $v (r) = 2 v_{a v} (1 - {(\frac{r}{R})}^{2})$

si3_e (Eq. 19.1)

19.3 Convection-Diffusion Equation for Cylindrical Cross-Sections

The second equation we require is the convection-diffusion equation (see Eq. 10.20) which we write (for convenience) for a single substance only, i.e., dropping the indices to find

$\begin{array}{l} \frac{\partial ρ}{\partial t} + \vec{\nabla} \cdot (ρ \vec{v}) & = D Δ ρ \\ \frac{\partial ρ}{\partial t} + ρ (\vec{\nabla} \vec{v}) + \vec{v} (\vec{\nabla} ρ) & = D Δ ρ \\ \frac{\partial ρ}{\partial t} + \vec{v} (\vec{\nabla} ρ) & = D Δ ρ \end{array}$ $\begin{array}{l} \frac{\partial ρ}{\partial t} + \vec{\nabla} \cdot (ρ \vec{v}) & = D Δ ρ \\ \frac{\partial ρ}{\partial t} + ρ (\vec{\nabla} \vec{v}) + \vec{v} (\vec{\nabla} ρ) & = D Δ ρ \\ \frac{\partial ρ}{\partial t} + \vec{v} (\vec{\nabla} ρ) & = D Δ ρ \end{array}$

si4_e (Eq. 19.2)

where we applied the continuity equation to find $\vec{\nabla} \vec{v} = 0$ $\vec{\nabla} \vec{v} = 0$ (see Eq. 10.16). In cylindrical coordinates, we can rewrite Eq. 19.2 using the fact that there is no velocity component in the azimuthal and radial direction. We use the gradient $\vec{\nabla} ρ$ $\vec{\nabla} ρ$ in cylindrical coordinates (see Eq. 7.90) as well as the scalar Laplacian in cylindrical coordinates (see Eq. 7.99) which yields

$\frac{\partial ρ}{\partial t} + v_{z} \frac{\partial ρ}{\partial z} = D (\frac{1}{r} \frac{\partial ρ}{\partial r} + \frac{\partial^{2} ρ}{\partial r^{2}} + \frac{\partial^{2} ρ}{\partial z^{2}})$ $\frac{\partial ρ}{\partial t} + v_{z} \frac{\partial ρ}{\partial z} = D (\frac{1}{r} \frac{\partial ρ}{\partial r} + \frac{\partial^{2} ρ}{\partial r^{2}} + \frac{\partial^{2} ρ}{\partial z^{2}})$

si7_e (Eq. 19.3)

19.4 Mass Concentration Function

The problem with Eq. 19.3 is that the concentration ρ is a function of time t, radius r, and z. However, as we have originally stated, we expect the diffusion in radial direction to proceed fairly quickly and thus, quickly smear out the radial concentration profiles. In this case we would be left only with ρ being a function of time t and z.

Obviously, this is only an assumption that we still have to proof. However, for the moment, we will assume that we will achieve in doing so and thus, we rewrite ρ as two functions:

• Radial average mass concentration $\bar{ρ} (t, z) = \frac{1}{π R^{2}} \int_{0}^{R} ρ (r, t, z) 2 πrdr$ $\bar{ρ} (t, z) = \frac{1}{π R^{2}} \int_{0}^{R} ρ (r, t, z) 2 πrdr$ which is defined as the average concentration of the concentration across the radius. It is a function only of time t and z.

• Radial mass concentration $\hat{ρ} (r, t, z)$ $\hat{ρ} (r, t, z)$ with $\frac{1}{π R^{2}} \int_{0}^{R} \hat{ρ} (r, t, z) 2 πrdr = 0$ $\frac{1}{π R^{2}} \int_{0}^{R} \hat{ρ} (r, t, z) 2 πrdr = 0$ as the function which “overlays” the average concentration to account for potential changes of the concentration across the radius. This then is a function of radius r, time t, and z.

As stated, we expect to proof $\hat{ρ}$ $\hat{ρ}$ to be negligibly small which will allow us to drop it during our calculations. Using these two functions, we now define the mass concentration function ρ as

$ρ = \bar{ρ} (t, z) + \hat{ρ} (r, t, z)$ $ρ = \bar{ρ} (t, z) + \hat{ρ} (r, t, z)$

(Eq. 19.4)

19.5 Convection-Diffusion Equation

Using Eq. 19.4 we can now rewrite Eq. 19.3 to

$\frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \hat{ρ}}{\partial t} + v_{z} (\frac{\partial \bar{ρ}}{\partial z} + \frac{\partial \hat{ρ}}{\partial z}) = D (\frac{1}{r} (\frac{\partial \bar{ρ}}{\partial r} + \frac{\partial \hat{ρ}}{\partial r}) + \frac{\partial^{2} \bar{ρ}}{\partial r^{2}} + \frac{\partial^{2} \hat{ρ}}{\partial r^{2}} + \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} + \frac{\partial^{2} \hat{ρ}}{\partial z^{2}})$ $\frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \hat{ρ}}{\partial t} + v_{z} (\frac{\partial \bar{ρ}}{\partial z} + \frac{\partial \hat{ρ}}{\partial z}) = D (\frac{1}{r} (\frac{\partial \bar{ρ}}{\partial r} + \frac{\partial \hat{ρ}}{\partial r}) + \frac{\partial^{2} \bar{ρ}}{\partial r^{2}} + \frac{\partial^{2} \hat{ρ}}{\partial r^{2}} + \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} + \frac{\partial^{2} \hat{ρ}}{\partial z^{2}})$

si13_e (Eq. 19.5)

We note that $\frac{\partial \bar{ρ}}{\partial r} = \frac{\partial^{2} \bar{ρ}}{\partial r^{2}} = 0$ $\frac{\partial \bar{ρ}}{\partial r} = \frac{\partial^{2} \bar{ρ}}{\partial r^{2}} = 0$ as we have defined $\bar{ρ}$ $\bar{ρ}$ as the radial average mass concentration and thus independent r. In this case Eq. 19.5 becomes

$\begin{array}{l} \frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \hat{ρ}}{\partial t} + v_{z} (\frac{\partial \bar{ρ}}{\partial z} + \frac{\partial \hat{ρ}}{\partial z}) = D (\frac{1}{r} \frac{\partial \hat{ρ}}{\partial r} + \frac{\partial^{2} \hat{ρ}}{\partial r^{2}} + \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} + \frac{\partial^{2} \hat{ρ}}{\partial z^{2}}) \\ \frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \hat{ρ}}{\partial t} + v_{z} (\frac{\partial \bar{ρ}}{\partial z} + \frac{\partial \hat{ρ}}{\partial z}) = D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r}) + \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} + \frac{\partial^{2} \hat{ρ}}{\partial z^{2}}) \end{array}$ $\begin{array}{l} \frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \hat{ρ}}{\partial t} + v_{z} (\frac{\partial \bar{ρ}}{\partial z} + \frac{\partial \hat{ρ}}{\partial z}) = D (\frac{1}{r} \frac{\partial \hat{ρ}}{\partial r} + \frac{\partial^{2} \hat{ρ}}{\partial r^{2}} + \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} + \frac{\partial^{2} \hat{ρ}}{\partial z^{2}}) \\ \frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \hat{ρ}}{\partial t} + v_{z} (\frac{\partial \bar{ρ}}{\partial z} + \frac{\partial \hat{ρ}}{\partial z}) = D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r}) + \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} + \frac{\partial^{2} \hat{ρ}}{\partial z^{2}}) \end{array}$

si16_e (Eq. 19.6)

Eq. 19.6 is the first of the two important equations we need for solving the Taylor-Aris dispersion problem. It contains all the physics involved in the problem. However, the equation is very complicated and difficult to solve. We will solve it by splitting it into two solutions: A solution for the radial mass concentration $\hat{ρ}$ $\hat{ρ}$ and a solution for the (radially) averaged mass concentration $\bar{ρ}$ $\bar{ρ}$ . For the latter we will perform averaging of Eq. 19.6.

19.5.1 Averaging Over the Cross-Section

We now average both sides of Eq. 19.6 across the radius which allows us to drop some terms

$\frac{1}{π R^{2}} \int_{0}^{R} (\frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \hat{ρ}}{\partial t} + v_{z} (\frac{\partial \bar{ρ}}{\partial z} + \frac{\partial \hat{ρ}}{\partial z})) 2 πrdr = \frac{1}{π R^{2}} \int_{0}^{R} D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r}) + \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} + \frac{\partial^{2} \hat{ρ}}{\partial z^{2}}) 2 πrdr$ $\frac{1}{π R^{2}} \int_{0}^{R} (\frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \hat{ρ}}{\partial t} + v_{z} (\frac{\partial \bar{ρ}}{\partial z} + \frac{\partial \hat{ρ}}{\partial z})) 2 πrdr = \frac{1}{π R^{2}} \int_{0}^{R} D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r}) + \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} + \frac{\partial^{2} \hat{ρ}}{\partial z^{2}}) 2 πrdr$

si19_e

Where we rewrite to

$\frac{2}{R^{2}} \int_{0}^{R} (\frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \hat{ρ}}{\partial t} + v_{z} (\frac{\partial \bar{ρ}}{\partial z} + \frac{\partial \hat{ρ}}{\partial z})) r d r = \frac{2 D}{R^{2}} \int_{0}^{R} D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r}) + \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} + \frac{\partial^{2} \hat{ρ}}{\partial z^{2}}) r d r$ $\frac{2}{R^{2}} \int_{0}^{R} (\frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \hat{ρ}}{\partial t} + v_{z} (\frac{\partial \bar{ρ}}{\partial z} + \frac{\partial \hat{ρ}}{\partial z})) r d r = \frac{2 D}{R^{2}} \int_{0}^{R} D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r}) + \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} + \frac{\partial^{2} \hat{ρ}}{\partial z^{2}}) r d r$

si20_e

where we split the expression into individual integrals finding

$\frac{2}{R^{2}} (\int_{0}^{R} \frac{\partial \bar{ρ}}{\partial t} r d r + \int_{0}^{R} \frac{\partial \hat{ρ}}{\partial t} r d r + \int_{0}^{R} v_{2} (\frac{\partial \bar{ρ}}{\partial z} + \frac{\partial \hat{ρ}}{\partial z}) r d r) = \frac{2 D}{R^{2}} (\int_{0}^{R} (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r}) r d r) + \int_{0}^{R} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} r d r + \int_{0}^{R} \frac{\partial^{2} \hat{ρ}}{\partial z^{2}} r d r)$ $\frac{2}{R^{2}} (\int_{0}^{R} \frac{\partial \bar{ρ}}{\partial t} r d r + \int_{0}^{R} \frac{\partial \hat{ρ}}{\partial t} r d r + \int_{0}^{R} v_{2} (\frac{\partial \bar{ρ}}{\partial z} + \frac{\partial \hat{ρ}}{\partial z}) r d r) = \frac{2 D}{R^{2}} (\int_{0}^{R} (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r}) r d r) + \int_{0}^{R} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} r d r + \int_{0}^{R} \frac{\partial^{2} \hat{ρ}}{\partial z^{2}} r d r)$

si21_e (Eq. 19.7)

We now use the fact that we can interchange the order of the integral and the differential and thus we can rewrite Eq. 19.7 to

$\frac{2}{R^{2}} (\frac{\partial}{\partial t} (\int_{0}^{R} \bar{ρ} r d r) + \frac{\partial}{\partial t} (\int_{0 .}^{R} \hat{ρ} r d r) + \int_{0}^{R} v_{z} \frac{\partial \bar{ρ}}{\partial z} r d r + \int_{0}^{R} v_{z} \frac{\partial \hat{ρ}}{\partial z} r d r) = \frac{2 D}{R^{2}} (\frac{\partial}{\partial r} (\int_{0}^{R} r \frac{\partial \hat{ρ}}{\partial r} d r) + \frac{\partial^{2}}{\partial z^{2}} (\int_{0}^{R} \bar{ρ} r d r) + \frac{\partial^{2}}{\partial z^{2}} (\int_{0}^{R} \hat{ρ} r d r))$ $\frac{2}{R^{2}} (\frac{\partial}{\partial t} (\int_{0}^{R} \bar{ρ} r d r) + \frac{\partial}{\partial t} (\int_{0 .}^{R} \hat{ρ} r d r) + \int_{0}^{R} v_{z} \frac{\partial \bar{ρ}}{\partial z} r d r + \int_{0}^{R} v_{z} \frac{\partial \hat{ρ}}{\partial z} r d r) = \frac{2 D}{R^{2}} (\frac{\partial}{\partial r} (\int_{0}^{R} r \frac{\partial \hat{ρ}}{\partial r} d r) + \frac{\partial^{2}}{\partial z^{2}} (\int_{0}^{R} \bar{ρ} r d r) + \frac{\partial^{2}}{\partial z^{2}} (\int_{0}^{R} \hat{ρ} r d r))$

si22_e (Eq. 19.8)

We will look at the individual terms one by one. We first note that the average value for the concentration across the cross-section is defined as

$\frac{2}{R^{2}} \int_{0}^{R} \bar{ρ} r d r = \bar{ρ}$ $\frac{2}{R^{2}} \int_{0}^{R} \bar{ρ} r d r = \bar{ρ}$

which is nothing other than the “average velocity of the average velocity”. This turns the first term on the left-hand side of Eq. 19.8 to

$\begin{array}{l} \frac{2}{R^{2}} \frac{\partial}{\partial t} (\int_{0}^{R} \bar{ρ} r d r) & = \frac{\partial}{\partial t} (\frac{2}{R^{2}} \int_{0}^{R} \bar{ρ} r d r) \\ = \frac{\partial \bar{ρ}}{\partial t} \end{array}$ $\begin{array}{l} \frac{2}{R^{2}} \frac{\partial}{\partial t} (\int_{0}^{R} \bar{ρ} r d r) & = \frac{\partial}{\partial t} (\frac{2}{R^{2}} \int_{0}^{R} \bar{ρ} r d r) \\ = \frac{\partial \bar{ρ}}{\partial t} \end{array}$

si24_e

Furthermore, we use the average velocity v_z across the cross-section defined as

$\frac{2}{R^{2}} \int_{0}^{R} v_{z} r d r = v_{a v}$ $\frac{2}{R^{2}} \int_{0}^{R} v_{z} r d r = v_{a v}$

This allows us to note that the second term on the left-hand side of Eq. 19.8 can be removed because

$\frac{2}{R^{2}} \int_{0}^{R} \hat{ρ} r d r = 0$ $\frac{2}{R^{2}} \int_{0}^{R} \hat{ρ} r d r = 0$

(Eq. 19.9)

which is the property by which we defined the function $\hat{ρ}$ $\hat{ρ}$ . The third term on the left-hand side of Eq. 19.8 can be rewritten as an integration by parts (see section 3.1.6.3) according to

$\begin{array}{l} \frac{2}{R^{2}} (\int_{0}^{R} v_{z} \frac{\partial \bar{ρ}}{\partial z} r d r) & = \frac{2}{R^{2}} (\frac{\partial \bar{ρ}}{\partial z} \int_{0}^{R} v_{z} r d r - \int_{0}^{R} v_{z} \frac{d}{d r} (\frac{\partial \bar{ρ}}{\partial z}) d r) \\ = \frac{\partial \bar{ρ}}{\partial z} v_{a v} - \frac{2}{R^{2}} \int_{0}^{R} v_{z} \frac{\partial}{\partial z} (\frac{d \bar{ρ}}{d r}) d r \\ = \frac{\partial \bar{ρ}}{\partial z} v_{a v} \end{array}$ $\begin{array}{l} \frac{2}{R^{2}} (\int_{0}^{R} v_{z} \frac{\partial \bar{ρ}}{\partial z} r d r) & = \frac{2}{R^{2}} (\frac{\partial \bar{ρ}}{\partial z} \int_{0}^{R} v_{z} r d r - \int_{0}^{R} v_{z} \frac{d}{d r} (\frac{\partial \bar{ρ}}{\partial z}) d r) \\ = \frac{\partial \bar{ρ}}{\partial z} v_{a v} - \frac{2}{R^{2}} \int_{0}^{R} v_{z} \frac{\partial}{\partial z} (\frac{d \bar{ρ}}{d r}) d r \\ = \frac{\partial \bar{ρ}}{\partial z} v_{a v} \end{array}$

si28_e

where we use the fact that $\bar{ρ}$ $\bar{ρ}$ is not a function of r and therefore $\frac{d \bar{ρ}}{d r} = 0$ $\frac{d \bar{ρ}}{d r} = 0$ . The last term on the left-hand side of Eq. 19.8 cannot be further simplified. However, we can drop the first term on the right-hand side of Eq. 19.8 if we carry out the following integration by parts

$\begin{array}{l} \int_{0}^{R} r \frac{\partial \hat{ρ}}{\partial r} d r & = r \int_{0}^{R} \frac{\partial \hat{ρ}}{\partial r} d r - \int \int \frac{\partial \hat{ρ}}{\partial r} d r d r \\ = r \int_{0}^{R} \frac{\partial \hat{ρ}}{\partial r} d r - r \int_{0}^{R} \frac{\partial \hat{ρ}}{\partial r} d r \\ = 0 \end{array}$ $\begin{array}{l} \int_{0}^{R} r \frac{\partial \hat{ρ}}{\partial r} d r & = r \int_{0}^{R} \frac{\partial \hat{ρ}}{\partial r} d r - \int \int \frac{\partial \hat{ρ}}{\partial r} d r d r \\ = r \int_{0}^{R} \frac{\partial \hat{ρ}}{\partial r} d r - r \int_{0}^{R} \frac{\partial \hat{ρ}}{\partial r} d r \\ = 0 \end{array}$

si31_e

The second term of the right-hand side of Eq. 19.8 also cannot be simplified further. The last term, on the other hand, can again be dropped because of Eq. 19.9. Using these simplifications, Eq. 19.8 thus becomes

$\frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \bar{ρ}}{\partial z} v_{a v} + \frac{2}{R^{2}} \int_{0}^{R} v_{z} \frac{\partial \hat{ρ}}{\partial z} r d r = D \frac{\partial^{2} \bar{ρ}}{\partial z^{2}}$ $\frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \bar{ρ}}{\partial z} v_{a v} + \frac{2}{R^{2}} \int_{0}^{R} v_{z} \frac{\partial \hat{ρ}}{\partial z} r d r = D \frac{\partial^{2} \bar{ρ}}{\partial z^{2}}$

(Eq. 19.10)

Eq. 19.10 is the second equation we require for solving the Taylor-Aris dispersion. As in Eq. 19.6 it is exact as we did not drop any terms. We merely performed local averaging which means that Eq. 19.10 will not give us information about the radial mass concentration profile, but merely about the averaged values. As you can see, the only term still involving $\hat{ρ}$ $\hat{ρ}$ is the integral that we were not quite able to drop yet. We will see in a moment how to cope with this.

19.5.2 Subtracting the Two Equations

We now subtract Eq. 19.10 from Eq. 19.6. By doing so, we are effectively left with an equation that only describes the radial mass concentration profile $\hat{ρ}$ $\hat{ρ}$ . This is due to the fact that by subtracting Eq. 19.6, we remove the average values. We then obtain

$\frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \hat{ρ}}{\partial t} + v_{z} (\frac{\partial \bar{ρ}}{\partial z} + \frac{\partial \hat{ρ}}{\partial z}) - \frac{\partial \bar{ρ}}{\partial t} - \frac{\partial \bar{ρ}}{\partial z} v_{a v} - \frac{2}{R^{2}} \int_{0}^{R} v_{z} \frac{\partial \hat{ρ}}{\partial z} r d r = D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r}) + \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} + \frac{\partial^{2} \hat{ρ}}{\partial z^{2}}) - D \frac{\partial^{2} \bar{ρ}}{\partial z^{2}}$ $\frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \hat{ρ}}{\partial t} + v_{z} (\frac{\partial \bar{ρ}}{\partial z} + \frac{\partial \hat{ρ}}{\partial z}) - \frac{\partial \bar{ρ}}{\partial t} - \frac{\partial \bar{ρ}}{\partial z} v_{a v} - \frac{2}{R^{2}} \int_{0}^{R} v_{z} \frac{\partial \hat{ρ}}{\partial z} r d r = D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r}) + \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} + \frac{\partial^{2} \hat{ρ}}{\partial z^{2}}) - D \frac{\partial^{2} \bar{ρ}}{\partial z^{2}}$

si35_e

from which we find

$\frac{\partial \hat{p}}{\partial t} + \frac{\partial \bar{p}}{\partial z} (v_{z} - v_{a v}) + v_{z} \frac{\partial \hat{p}}{\partial z} - \frac{2}{R^{2}} \int_{0}^{R} v_{z} \frac{\partial \hat{p}}{\partial z} r d r = D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{p}}{\partial r}) + \frac{\partial^{2} \hat{p}}{\partial z^{2}})$ $\frac{\partial \hat{p}}{\partial t} + \frac{\partial \bar{p}}{\partial z} (v_{z} - v_{a v}) + v_{z} \frac{\partial \hat{p}}{\partial z} - \frac{2}{R^{2}} \int_{0}^{R} v_{z} \frac{\partial \hat{p}}{\partial z} r d r = D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{p}}{\partial r}) + \frac{\partial^{2} \hat{p}}{\partial z^{2}})$

si36_e (Eq. 19.11)

This equation is exact, as it does not contain any simplifications that would have involved dropping terms. It will give us the concentration $\hat{ρ}$ $\hat{ρ}$ as a function of the radius and the coordinate z. As you can see, Eq. 19.11 still contains one term that refers to the average concentration, i.e., $\frac{\partial \bar{ρ}}{\partial z}$ $\frac{\partial \bar{ρ}}{\partial z}$ . Therefore, in order to solve Eq. 19.11, we need to make simplifications.

19.6 Solving for $\hat{ρ}$ $\hat{ρ}$

19.6.1 Coordinate System Transformation

Before solving the equation, we first define a new velocity function that combines the two independent functions t and z. We term this velocity function

$v_{z}^{'} = v_{z} - v_{a v} = v_{a v} (1 - 2 {(\frac{r}{R})}^{2})$ $v_{z}^{'} = v_{z} - v_{a v} = v_{a v} (1 - 2 {(\frac{r}{R})}^{2})$

si40_e (Eq. 19.12)

This effectively corresponds to a coordinate system transformation where we use the independent variable

$z' = z - v_{a v} t$ $z' = z - v_{a v} t$

This coordinate system travels with the control volume at the constant velocity v_av.

19.6.2 Gradients Along z-axis

We first note that we expect all gradients along the z-axis to be significantly smaller than gradients along the radial directions. Doing so will allow us to drop all terms containing $\frac{\partial}{\partial z}$ $\frac{\partial}{\partial z}$ from Eq. 19.11 as we can set $\frac{\partial}{\partial z} = \frac{\partial^{2}}{\partial z^{2}} = 0$ $\frac{\partial}{\partial z} = \frac{\partial^{2}}{\partial z^{2}} = 0$ . We are then left to solve

$\frac{\partial \hat{ρ}}{\partial t} + \frac{\partial \bar{ρ}}{\partial z} v_{z}^{'} = D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r}))$ $\frac{\partial \hat{ρ}}{\partial t} + \frac{\partial \bar{ρ}}{\partial z} v_{z}^{'} = D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r}))$

si44_e (Eq. 19.13)

Obviously, this assumption can be made as long as the length L of the microfluidic channel investigated is very much larger than the diameter R, i.e., L .R.

19.6.3 Time-Dependency

The next simplification we make is by assuming that the diffusion in radial direction is fast. If this assumption holds, then we may be able to drop the time-dependent term $\frac{\partial \hat{ρ}}{\partial t}$ $\frac{\partial \hat{ρ}}{\partial t}$ . By doing this, we lose the time-dependency of our solution, i.e., we make our solution independent of time. If a function of time is required, then this simplification cannot be made. However, by dimensionless analysis, we will quickly discuss which timescales we are dealing with. Obviously we would require

$\frac{\partial \hat{ρ}}{\partial t} « D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r}))$ $\frac{\partial \hat{ρ}}{\partial t} « D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r}))$

si46_e

By making this expression dimensionless we find

$\frac{1}{t_{char}} \frac{\partial \hat{ρ}}{\partial t} « \frac{D}{R^{2}} (\frac{1}{R} \frac{\partial}{\partial R} (R \frac{\partial \hat{ρ}}{\partial R}))$ $\frac{1}{t_{char}} \frac{\partial \hat{ρ}}{\partial t} « \frac{D}{R^{2}} (\frac{1}{R} \frac{\partial}{\partial R} (R \frac{\partial \hat{ρ}}{\partial R}))$

si47_e

from which we find

$t_{char} » \frac{R^{2}}{D}$ $t_{char} » \frac{R^{2}}{D}$

Just to give a rough idea what the timescale actually is: for a molecule of ethanol in water (see Tab. 6.9) in a capillary with a diameter of 100 μm, the characteristic timescale is around 0.4 μs. If the diffusion is to be determined around this precision, then it is not possible to drop the time dependency. As you can see, this is a reasonably small timescale, and in most cases it is sufficient to study diffusion at longer timescales, thus the simplification is appropriate. We are thus able to simplify Eq. 19.13 to

$\begin{array}{l} v_{z}^{'} \frac{\partial \bar{ρ}}{\partial z} & = D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r})) \\ v_{z}^{'} \frac{\partial \bar{ρ}}{\partial z} |_{z'} & = D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r})) \end{array}$ $\begin{array}{l} v_{z}^{'} \frac{\partial \bar{ρ}}{\partial z} & = D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r})) \\ v_{z}^{'} \frac{\partial \bar{ρ}}{\partial z} |_{z'} & = D (\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial \hat{ρ}}{\partial r})) \end{array}$

si49_e (Eq. 19.14)

This is where the coordinate system transformation we made really helps. In our new frame of reference $\hat{ρ}$ $\hat{ρ}$ is a function only of the radius. The gradient of $\hat{ρ}$ $\hat{ρ}$ along the z-axis in this frame of reference is constant, which is why we can simply write $\frac{\partial \hat{ρ}}{\partial z} |_{z'}$ $\frac{\partial \hat{ρ}}{\partial z} |_{z'}$ which states “treat $\frac{\partial \hat{ρ}}{\partial z}$ $\frac{\partial \hat{ρ}}{\partial z}$ as a constant with respect to z′ and simply evaluate it at the location z′”.. We will now solve this equation for $\hat{ρ}$ $\hat{ρ}$ with the radius r as the independent variable. We also note that $\bar{ρ}$ $\bar{ρ}$ is not a function of the radius, as it is the average value along the radius. Therefore we can rewrite Eq. 19.14 to an ODE which reads

$\begin{array}{l} \frac{\partial \bar{ρ}}{\partial z} |_{z'} v_{a v} (1 - \frac{2 r^{2}}{R^{2}}) & = D (\frac{1}{r} \frac{d}{d r} (r \frac{d \hat{ρ}}{d r})) \\ \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (r - \frac{2 r^{3}}{R^{2}}) d r & = d (r \frac{d \hat{ρ}}{d r}) \\ \int \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (r - \frac{2 r^{3}}{R^{2}}) d r & = \int d (r \frac{d \hat{ρ}}{d r}) \\ \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{r^{2}}{2} - \frac{r^{4}}{2 R^{2}}) + c_{1} & = r \frac{d \hat{ρ}}{d r} \end{array}$ $\begin{array}{l} \frac{\partial \bar{ρ}}{\partial z} |_{z'} v_{a v} (1 - \frac{2 r^{2}}{R^{2}}) & = D (\frac{1}{r} \frac{d}{d r} (r \frac{d \hat{ρ}}{d r})) \\ \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (r - \frac{2 r^{3}}{R^{2}}) d r & = d (r \frac{d \hat{ρ}}{d r}) \\ \int \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (r - \frac{2 r^{3}}{R^{2}}) d r & = \int d (r \frac{d \hat{ρ}}{d r}) \\ \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{r^{2}}{2} - \frac{r^{4}}{2 R^{2}}) + c_{1} & = r \frac{d \hat{ρ}}{d r} \end{array}$

si56_e

where we have used Eq. 19.12 for the velocity. We also have to find the first integration constant that we obtain by noting that the profile we expect is symmetrical and therefore $\frac{d \hat{ρ}}{d r} |_{0} = 0$ $\frac{d \hat{ρ}}{d r} |_{0} = 0$ , in which case we find c₁ = 0. The second integration yields

$\begin{array}{l} \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{r}{2} - \frac{r^{3}}{2 R^{2}}) d r & = d \hat{ρ} \\ \int \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{r}{2} - \frac{r^{3}}{2 R^{2}}) d r & = \int d \hat{ρ} \\ \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{r^{2}}{4} - \frac{r^{4}}{8 R^{2}}) + c_{2} & = \hat{ρ} \end{array}$ $\begin{array}{l} \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{r}{2} - \frac{r^{3}}{2 R^{2}}) d r & = d \hat{ρ} \\ \int \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{r}{2} - \frac{r^{3}}{2 R^{2}}) d r & = \int d \hat{ρ} \\ \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{r^{2}}{4} - \frac{r^{4}}{8 R^{2}}) + c_{2} & = \hat{ρ} \end{array}$

si58_e (Eq. 19.15)

The second integration constant is found using the fact that we defined $\hat{ρ}$ $\hat{ρ}$ with the property

$\begin{array}{l} \frac{1}{π R^{2}} \int_{0}^{R} \hat{ρ} 2 πrdr & = 0 \\ \frac{2}{R^{2}} \int_{0}^{R} \hat{ρ} r d r & = 0 \end{array}$ $\begin{array}{l} \frac{1}{π R^{2}} \int_{0}^{R} \hat{ρ} 2 πrdr & = 0 \\ \frac{2}{R^{2}} \int_{0}^{R} \hat{ρ} r d r & = 0 \end{array}$

si60_e (Eq. 19.16)

in which case we find

$\begin{array}{l} \frac{2}{R^{2}} \int_{0}^{R} \hat{ρ} r d r & = \frac{2}{R^{2}} \int_{0}^{R} (\frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{r^{3}}{4} - \frac{r^{5}}{8 R^{2}}) + c_{2} r) d r \\ = \frac{2}{R^{2}} {[\frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{r^{4}}{16} - \frac{r^{6}}{48 R^{2}}) + c_{2} \frac{r^{2}}{2}]}_{0}^{R} \\ = \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{R^{2}}{8} - \frac{R^{4}}{24 R^{2}}) + c_{2} \\ \overset{!}{=} 0 \\ c_{2} & = - \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{R^{2}}{8} - \frac{R^{2}}{24}) \\ = - \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} \frac{R^{2}}{12} \end{array}$ $\begin{array}{l} \frac{2}{R^{2}} \int_{0}^{R} \hat{ρ} r d r & = \frac{2}{R^{2}} \int_{0}^{R} (\frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{r^{3}}{4} - \frac{r^{5}}{8 R^{2}}) + c_{2} r) d r \\ = \frac{2}{R^{2}} {[\frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{r^{4}}{16} - \frac{r^{6}}{48 R^{2}}) + c_{2} \frac{r^{2}}{2}]}_{0}^{R} \\ = \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{R^{2}}{8} - \frac{R^{4}}{24 R^{2}}) + c_{2} \\ \overset{!}{=} 0 \\ c_{2} & = - \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{R^{2}}{8} - \frac{R^{2}}{24}) \\ = - \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} \frac{R^{2}}{12} \end{array}$

si61_e

in which case we find the solution for the radial concentration $\hat{ρ}$ $\hat{ρ}$ from Eq. 19.6 as

$\begin{array}{l} \hat{ρ} & = \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{r^{2}}{4} - \frac{r^{4}}{8 R^{2}} - \frac{R^{2}}{12}) \\ = \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{24 D} (6 r^{2} - 3 \frac{r^{4}}{R^{2}} - 2 R^{2}) \end{array}$ $\begin{array}{l} \hat{ρ} & = \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{D} (\frac{r^{2}}{4} - \frac{r^{4}}{8 R^{2}} - \frac{R^{2}}{12}) \\ = \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{24 D} (6 r^{2} - 3 \frac{r^{4}}{R^{2}} - 2 R^{2}) \end{array}$

si63_e (Eq. 19.17)

$= \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{4 D} (r^{2} - \frac{1}{2} \frac{r^{4}}{R^{2}} - \frac{1}{3} R^{2})$ $= \frac{\partial \bar{ρ}}{\partial z} |_{z'} \frac{v_{a v}}{4 D} (r^{2} - \frac{1}{2} \frac{r^{4}}{R^{2}} - \frac{1}{3} R^{2})$

si64_e (Eq. 19.18)

Eq. 19.17 gives us an analytical expression to calculate the concentration of the radial concentration in our moving reference system. This equation is by itself not very useful, but we require it to evaluate the integral term we dropped earlier from Eq. 19.14.

19.6.4 Validity of the “fast radial diffusion” Assumption

We have now found the equation that describes the diffusion in radial direction. Our original assumption was that we can neglect this diffusion given a “sufficiently rough” timescale. From Eq. 19.18, we can derive when this assumption is correct. It is correct if $\hat{ρ}$ $\hat{ρ}$ is very small. In order to satisfy this for any position along the radius r, we therefore require the prefactor to be very small thus

$\frac{R^{2} v_{a v}}{4 D} \frac{\partial \hat{ρ}}{\partial z} |_{z'} « 1$ $\frac{R^{2} v_{a v}}{4 D} \frac{\partial \hat{ρ}}{\partial z} |_{z'} « 1$

si66_e

where for the term $\frac{\partial \hat{ρ}}{\partial z} |_{z'}$ $\frac{\partial \hat{ρ}}{\partial z} |_{z'}$ really only the total length L of the channel we are investigating is of importance, as this defines the steepness of the gradient. We can therefore note that

$\begin{array}{l} \frac{R^{2} v_{a v}}{4 D} \frac{1}{L} « 1 \\ P e \frac{R}{4 L} « 1 \\ P e « 4 \frac{L}{R} \end{array}$ $\begin{array}{l} \frac{R^{2} v_{a v}}{4 D} \frac{1}{L} « 1 \\ P e \frac{R}{4 L} « 1 \\ P e « 4 \frac{L}{R} \end{array}$

si68_e (Eq. 19.19)

Eq. 19.19 tells us that in our assumption that we may be able to neglect the time-dependency in radial diffusion only if the Péclet number is significantly smaller than $4 \frac{L}{R}$ $4 \frac{L}{R}$ .

19.7 Solving for $\bar{ρ}$ $\bar{ρ}$

We have found the solution to the radial mass concentration profiles with Eq. 19.17. We therefore have a solution to one of our original equations, i.e., Eq. 19.6. However, we still need to find a solution for the average concentration $\bar{ρ}$ $\bar{ρ}$ , i.e., to Eq. 19.10. For this we first need to evaluate the integral in Eq. 19.10 that we could not do before, as we require $\hat{ρ}$ $\hat{ρ}$ for this. We just found this expression (see Eq. 19.18) so now we can proceed to find

$\frac{2}{R^{2}} \int_{0}^{R} v_{z} \frac{\partial \hat{ρ}}{\partial z} r d r$ $\frac{2}{R^{2}} \int_{0}^{R} v_{z} \frac{\partial \hat{ρ}}{\partial z} r d r$

(Eq. 19.20)

For this we first need to find $\frac{\partial \hat{ρ}}{\partial z}$ $\frac{\partial \hat{ρ}}{\partial z}$ using Eq. 19.17 which is

$\frac{\partial \hat{ρ}}{\partial z} = \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} \frac{v_{a v}}{24 D} (6 r^{2} - 3 \frac{r^{4}}{R^{2}} - 2 R^{2})$ $\frac{\partial \hat{ρ}}{\partial z} = \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} \frac{v_{a v}}{24 D} (6 r^{2} - 3 \frac{r^{4}}{R^{2}} - 2 R^{2})$

si75_e (Eq. 19.21)

Using Eq. 19.12 and Eq. 19.21 we can evaluate Eq. 19.20 to

$\begin{array}{l} \frac{2}{R^{2}} \int_{0}^{R} v_{z} \frac{\partial \hat{ρ}}{\partial z} r d r & = \frac{4}{R^{2}} \int_{0}^{R} \frac{v_{a v}^{2}}{24 D} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} (r - {(\frac{r^{2}}{R})}^{2}) (6 r^{2} - 3 \frac{r^{4}}{R^{2}} - 2 R^{2}) d r \\ = \frac{v_{a v}^{2}}{6 R^{2} D} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} \int_{0}^{R} (6 r^{3} - 3 \frac{r^{5}}{R^{2}} - 2 r R^{2} + - \frac{6}{R^{2}} r^{5} + 3 \frac{r^{7}}{R^{4}} + 2 r^{3}) d r \\ = \frac{v_{a v}^{2}}{6 R^{2} D} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} {[\frac{3}{2} r^{4} - 3 \frac{r^{6}}{6 R^{2}} - R^{2} r^{2} + - \frac{6}{6 R^{2}} r^{6} + 3 \frac{r^{8}}{8 R^{4}} + \frac{2}{4} r^{4}]}_{0}^{R} \\ = \frac{v_{a v}^{2}}{6 R^{2} D} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} (\frac{3}{2} R^{4} - \frac{1}{2} R^{4} - R^{4} - R^{4} + \frac{3}{8} R^{4} + \frac{1}{2} R^{4}) \\ = \frac{R^{2} v_{a v}^{2}}{6 D} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} (\frac{3}{2} - \frac{1}{2} - 2 + \frac{3}{8} + \frac{1}{2}) \\ = - \frac{1}{48} \frac{R^{2} v_{a v}^{2}}{D} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} \end{array}$ $\begin{array}{l} \frac{2}{R^{2}} \int_{0}^{R} v_{z} \frac{\partial \hat{ρ}}{\partial z} r d r & = \frac{4}{R^{2}} \int_{0}^{R} \frac{v_{a v}^{2}}{24 D} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} (r - {(\frac{r^{2}}{R})}^{2}) (6 r^{2} - 3 \frac{r^{4}}{R^{2}} - 2 R^{2}) d r \\ = \frac{v_{a v}^{2}}{6 R^{2} D} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} \int_{0}^{R} (6 r^{3} - 3 \frac{r^{5}}{R^{2}} - 2 r R^{2} + - \frac{6}{R^{2}} r^{5} + 3 \frac{r^{7}}{R^{4}} + 2 r^{3}) d r \\ = \frac{v_{a v}^{2}}{6 R^{2} D} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} {[\frac{3}{2} r^{4} - 3 \frac{r^{6}}{6 R^{2}} - R^{2} r^{2} + - \frac{6}{6 R^{2}} r^{6} + 3 \frac{r^{8}}{8 R^{4}} + \frac{2}{4} r^{4}]}_{0}^{R} \\ = \frac{v_{a v}^{2}}{6 R^{2} D} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} (\frac{3}{2} R^{4} - \frac{1}{2} R^{4} - R^{4} - R^{4} + \frac{3}{8} R^{4} + \frac{1}{2} R^{4}) \\ = \frac{R^{2} v_{a v}^{2}}{6 D} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} (\frac{3}{2} - \frac{1}{2} - 2 + \frac{3}{8} + \frac{1}{2}) \\ = - \frac{1}{48} \frac{R^{2} v_{a v}^{2}}{D} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} \end{array}$

si76_e (Eq. 19.22)

We can see that this integral can be evaluated to a pretty simple expression. This we can set into Eq. 19.10 to obtain the solution to $\hat{ρ}$ $\hat{ρ}$ which is

$\begin{array}{l} \frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \bar{ρ}}{\partial z} v_{a v} - \frac{1}{48} \frac{R^{2} v_{a v}^{2}}{D} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} & = D \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} \\ \frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \bar{ρ}}{\partial z} v_{a v} & = (D + \frac{R^{2} v_{a v}^{2}}{48 D}) \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} \end{array}$ $\begin{array}{l} \frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \bar{ρ}}{\partial z} v_{a v} - \frac{1}{48} \frac{R^{2} v_{a v}^{2}}{D} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} & = D \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} \\ \frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \bar{ρ}}{\partial z} v_{a v} & = (D + \frac{R^{2} v_{a v}^{2}}{48 D}) \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} \end{array}$

si78_e (Eq. 19.23)

$\begin{array}{c} = D (1 + \frac{1}{48} {P e}^{2}) \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} \end{array}$ $\begin{array}{c} = D (1 + \frac{1}{48} {P e}^{2}) \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} \end{array}$

(Eq. 19.24)

$\begin{array}{c} \frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \bar{ρ}}{\partial z} v_{a v} & = D_{e ff .} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} \end{array}$ $\begin{array}{c} \frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \bar{ρ}}{\partial z} v_{a v} & = D_{e ff .} \frac{\partial^{2} \bar{ρ}}{\partial z^{2}} \end{array}$

(Eq. 19.25)

Eq. 19.24 is the PDE of the Taylor-Aris dispersion. It describes the t- and z-dependent changes of the average concentration v_av. This equation is very interesting because it effectively is (again) a convection-diffusion equation as defined in section 10.6, but with a slightly changed effective diffusion coefficient which is dependent on the velocity and the geometry of the flow channel. If this effective diffusion D_eff. coefficient is significantly bigger than the diffusion coefficient D alone, then molecular diffusion can be neglected and the convective dispersion is dominant. As can be seen from Eq. 19.24 this happens if

$\begin{array}{l} \frac{R^{2} v_{a v}^{2}}{48 D} » D \\ \frac{R^{2} v_{a v}^{2}}{D^{2}} » 48 \\ P e » \sqrt{48} \end{array}$ $\begin{array}{l} \frac{R^{2} v_{a v}^{2}}{48 D} » D \\ \frac{R^{2} v_{a v}^{2}}{D^{2}} » 48 \\ P e » \sqrt{48} \end{array}$

si81_e (Eq. 19.26)

Interpretation of the Effective Diffusion Coefficient. The effective diffusion coefficient D_eff. is dependent on the Péclet number Pe. If we consider a given fixed velocity v_av and radius R, the effective diffusion coefficient is at a minimum if the first derivative of Eq. 19.23 with respect to the D is 0

$\begin{array}{l} \frac{d}{d D} D_{eff .} & = \frac{d}{d D} (D + \frac{R^{2} v_{a v}^{2}}{48 D}) \\ = 1 - \frac{R^{2} v_{a v}^{2}}{48 D^{2}} \overset{!}{=} 0 \end{array}$ $\begin{array}{l} \frac{d}{d D} D_{eff .} & = \frac{d}{d D} (D + \frac{R^{2} v_{a v}^{2}}{48 D}) \\ = 1 - \frac{R^{2} v_{a v}^{2}}{48 D^{2}} \overset{!}{=} 0 \end{array}$

si82_e

This is fulfilled if

$D^{2} = \frac{R^{2} v_{a v}^{2}}{412} \to D_{min} = \frac{R v_{a v}}{\sqrt{48}}$ $D^{2} = \frac{R^{2} v_{a v}^{2}}{412} \to D_{min} = \frac{R v_{a v}}{\sqrt{48}}$

si83_e

At the value D = D_min the effective diffusion coefficient is minimal. At values of D > D_min the effective diffusion increase, which is intuitive from Eq. 19.23. This effect is purely due to molecular diffusion, which increases as the diffusion coefficient increase. However, for values D < D_min the effective diffusion coefficient also increases. This increase of effective diffusion at low diffusion coefficients is due to Taylor-Aris dispersion. As the molecular diffusion is low, dispersion has time to stretch the plug according to the Hagen-Poiseuille flow profile. This increases the interaction area and keeps the concentration gradient high. This, in turn, will guarantee effective diffusive mixing even if the diffusion coefficients are low.

19.8 Validity of the Solution

During our derivation, we have found two conditions that need to be fulfilled in order to have valid results with the derivations we took in Eq. 19.19 and Eq. 19.26. Therefore we can state that the Taylor-Aris dispersion is dominant in the range of

$4 \frac{L}{R} » P e » \sqrt{48}$ $4 \frac{L}{R} » P e » \sqrt{48}$

(Eq. 19.27)

In order to assess whether or not Taylor-Aris dispersion actually plays a significant role in a convection- diffusion problem at hand, Eq. 19.27 must be checked first.

19.9 Example

In the next step we will make an example solving Eq. 19.24 for a practical example. For this we assume a microfluidic channel cross-section with a smaller injection channel crossing a larger detection channel (see Fig. 19.3). An example of such a channel geometry would be, e.g., a plug injector which allows inserting a small volume of a fluid into the detection channel by first applying a pressure drop across the injection channel. At a given time, the pressure drop is removed and a pressure drop applied along the detection channel, thus moving the fraction of the liquid within the channel intersection into the detection channel. We will assume that the injected liquid is a sample of, e.g., an analyte in a buffer. A given amount of substance is to be transferred into the detection channel and then moved, e.g., toward a detector. As the longer channel is filled with the same buffer, we expect to see Taylor-Aris dispersion along the channel. For the sake of simplicity, we will assume the channel cross-sections to be circular. Thus we have to solve Eq. 19.24.

f19-03-9781455731411 — Fig. 19.3 Example of Taylor-Aris dispersion using a microfluidic plug injector for an analytical application. a/b) A microfluidic channel with a short injection and a long detection channel is first filled by applying a pressure drop ∆p along the injection channel. c/d) By switching the pressure drop to the detection channel, a small volume of sample is transferred into the detection channel. While being transported in the detection channel, it is diluted due to Taylor-Aris dispersion.

19.9.1 Coordinate System Transformation

We need to find a solution to Eq. 19.25 for which we will apply a coordinate system transformation that we have used previously. As can be seen from Eq. 19.25, the convection term $\frac{\partial \bar{p}}{\partial z} v_{a v}$ $\frac{\partial \bar{p}}{\partial z} v_{a v}$ contains the constant velocity vav. If we use a coordinate system that moves along with the center of the plug, we will have a very neat simplification. For this we note

$z^{'} = z - υ_{a v} t \to z = z^{'} + υ_{a v} t$ $z^{'} = z - υ_{a v} t \to z = z^{'} + υ_{a v} t$

(Eq. 19.28)

Using the chain rule of differentiation (see Eq. 3.14) we can now obtain the following differential which we require for Eq. 19.25

$\frac{\partial \bar{ρ}}{\partial z} = \frac{\partial \bar{ρ}}{\partial z^{'}} \frac{\partial z^{'}}{\partial z} = \frac{\partial \bar{ρ}}{\partial z^{'}}$ $\frac{\partial \bar{ρ}}{\partial z} = \frac{\partial \bar{ρ}}{\partial z^{'}} \frac{\partial z^{'}}{\partial z} = \frac{\partial \bar{ρ}}{\partial z^{'}}$

(Eq. 19.29)

Using Eq. 19.29 we can rewrite Eq. 19.25 to

$\begin{array}{l} \frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \bar{ρ}}{\partial z^{'}} (υ_{a v} - υ_{a v}) & = & D_{eff .} \frac{\partial^{2} \bar{ρ}}{\partial {z^{'}}^{2}} \\ \frac{\partial \bar{ρ}}{\partial t} & = & D_{eff .} \frac{\partial^{2} \bar{ρ}}{\partial {z^{'}}^{2}} \end{array}$ $\begin{array}{l} \frac{\partial \bar{ρ}}{\partial t} + \frac{\partial \bar{ρ}}{\partial z^{'}} (υ_{a v} - υ_{a v}) & = & D_{eff .} \frac{\partial^{2} \bar{ρ}}{\partial {z^{'}}^{2}} \\ \frac{\partial \bar{ρ}}{\partial t} & = & D_{eff .} \frac{\partial^{2} \bar{ρ}}{\partial {z^{'}}^{2}} \end{array}$

si88_e (Eq. 19.30)

which is nothing other than the one-dimensional diffusion equation (see Eq. 6.90).

In section 8.3.6 and section 8.3.8, we have already found several solutions to the diffusion equations including the cases for point source diffusion (see section 8.3.6.6) and limited plane diffusion (see section 8.3.6.7). All of these solutions would be directly applicable to Eq. 19.30 and would yield solutions in the moving coordinate system, i.e., as functions of z', which would be retransformed into the static coordinate system, i.e., to a function of z using Eq. 19.28.

In the following we will solve Eq. 19.30 for the Taylor-Aris dispersion problem shown in Fig. 19.3. As we will see, this is surprisingly simple given the solutions we already derived for more general PDEs.

19.9.2 Substitution and Separation of Variables

We will take a substitution and separation of variables approach similar to the approach taken in section 8.3.8. We therefore assume $\bar{p} (t, z) = T (t) Z (z')$ $\bar{p} (t, z) = T (t) Z (z')$ in which case Eq. 19.25 becomes

$\begin{array}{l} z \frac{\partial T}{\partial t} & = & D_{eff .} T \frac{\partial^{2} z}{\partial z^{2}} \\ \frac{1}{D_{eff}} \frac{1}{T} \frac{\partial T}{\partial t} & = & \frac{1}{z} \frac{\partial^{2} z}{\partial z^{2}} = - λ \end{array}$ $\begin{array}{l} z \frac{\partial T}{\partial t} & = & D_{eff .} T \frac{\partial^{2} z}{\partial z^{2}} \\ \frac{1}{D_{eff}} \frac{1}{T} \frac{\partial T}{\partial t} & = & \frac{1}{z} \frac{\partial^{2} z}{\partial z^{2}} = - λ \end{array}$

si90_e (Eq. 19.31)

where we require both sides of the equation to be constant. For details on this approach, please refer to section 8.3.3.3. Using Tab. 8.2 we can immediately derive the solutions of the two ODEs which we derive from Eq. 19.31. They are given by

$T (t) = c_{0} e^{- λ D_{e ff .} t}$ $T (t) = c_{0} e^{- λ D_{e ff .} t}$

(Eq. 19.32)

$z (z^{'}) = c_{1} sin (\sqrt{λ} z^{'}) + c_{2} cos (\sqrt{λ} z^{'})$ $z (z^{'}) = c_{1} sin (\sqrt{λ} z^{'}) + c_{2} cos (\sqrt{λ} z^{'})$

(Eq. 19.33)

$\bar{ρ} = (t, z^{'}) = e^{- {λ D}_{eff .} t} (c_{3} sin (\sqrt{λ} z^{'}) + c_{4} cos (\sqrt{λ} z^{'}))$ $\bar{ρ} = (t, z^{'}) = e^{- {λ D}_{eff .} t} (c_{3} sin (\sqrt{λ} z^{'}) + c_{4} cos (\sqrt{λ} z^{'}))$

(Eq. 19.34)

Obviously we need boundary conditions or initial conditions in order to determine A and the constants c3 and c4. For this we refer to Fig. 19.3. In Fig. 19.3b the injection channel is filled and the pressure drop is removed in order to be applied along the detection channel (Fig. 19.3c). We now study the propagation of the plug along the detection channel with the scenario in Fig. 19.3b serving as the distribution at t = 0. At this moment, there is a sharply defined rectangular concentration profile at the origin of our coordinate system that we set to the intersection. We choose the width of both channels to be w in which case we find

$\bar{ρ} (0, z^{'}) = \{\begin{array}{c} \bar{ρ} & - \frac{ω}{2} \leq z^{'} \leq + \frac{ω}{2} \\ 0 & otherwise \end{array}$ $\bar{ρ} (0, z^{'}) = \{\begin{array}{c} \bar{ρ} & - \frac{ω}{2} \leq z^{'} \leq + \frac{ω}{2} \\ 0 & otherwise \end{array}$

si94_e

This is the initial condition which we will use to solve for A and our constants. The problem at hand does not have any useful boundary conditions, as there is no boundary in this system which would give us, e.g., a zero-concentration condition. However, as we will see, the initial values are sufficient.

19.9.3 Balanced Rectangular Function

In order to correctly model the initial values of the system, we require a suitable function to represent the initial conditions. As Eq. 19.33 is a function of sine and cosine, it seems suitable to develop the initial distribution to a sine or cosine Fourier series. In the following we choose a cosine expansion, but a sine expansion would work just the same. The procedure to develop such a series expansion is explained in section 4.3.3. We will require Eq. 4.24, Eq. 4.26, Eq. 4.27, and Eq. 4.28 for this expansion. In the following we will derive a more general case of a function, which is a rectangular pulse around the origin that can be developed into a "nice" cosine Fourier series. The following definition may seem a little arbitrary, but it will make sense in just a moment. We define the rectangular function as

$\begin{array}{c} f_{b a l . rect .} (x) & = \end{array} \{\begin{cases} - 1 & - s - \frac{ω}{2} \leq x \leq - s \\ 1 & - \frac{ω}{2} \leq x \leq \frac{ω}{2} \\ - 1 & s - \frac{ω}{2} \leq x \leq s \\ 0 & otherwise \end{cases}$ $\begin{array}{c} f_{b a l . rect .} (x) & = \end{array} \{\begin{cases} - 1 & - s - \frac{ω}{2} \leq x \leq - s \\ 1 & - \frac{ω}{2} \leq x \leq \frac{ω}{2} \\ - 1 & s - \frac{ω}{2} \leq x \leq s \\ 0 & otherwise \end{cases}$

si95_e

where we have introduced two variables, i.e., the spacing s and the width of the rectangular pulse w. We term this function balanced rectangular function. It is depicted exemplarily in Fig. 19.4. It consists of a series of three rectangular pulses symmetric to the origin with the outer two pulses having the height —1 and the width w. The central pulse has the height 1 and the width w. The outer pulses are spaced by the distance s. In the following we will develop this function into a Fourier series. As we will see, the shape and the location of the pulses will account for a very practical series.

f19-04-9781455731411 — Fig. 19.4 Balanced rectangular function f_bai.rect. defined as a series of two pulses with height − 1 of width $\frac{w}{2} s$ $\frac{w}{2} s$ ws spaced by the spacing distance s. The central pulse (which is the one we are interested in) has the width w and the height + 1. This function can be developed to a very simple Fourier cosine series.

19.9.4 Fourier Series Expansion of the Balanced Rectangular Function

Referring to Eq. 4.26, Eq. 4.27, and Eq. 4.28 we find the Fourier coefficients for the balanced rectangular function to be

$\begin{array}{l} a_{0} & = & \frac{2}{s} \int_{- \frac{s}{2}}^{\frac{s}{2}} f_{b a l . rect .} (x) d x \\ = & \frac{2}{s} (\int_{- \frac{s}{2}}^{- \frac{s}{2} + \frac{w}{2}} (- 1) d x + \int_{- \frac{w}{2}}^{\frac{w}{2}} 1 d x + \int_{\frac{s}{2} - \frac{w}{2}}^{\frac{s}{2}} (- 1) d x) \\ = & \frac{2}{s} (- (- \frac{s}{2} + \frac{w}{2} + \frac{s}{2}) + (\frac{w}{2} + \frac{w}{2}) - (\frac{s}{2} - \frac{s}{2} + \frac{w}{2})) \\ = & 0 \end{array}$ $\begin{array}{l} a_{0} & = & \frac{2}{s} \int_{- \frac{s}{2}}^{\frac{s}{2}} f_{b a l . rect .} (x) d x \\ = & \frac{2}{s} (\int_{- \frac{s}{2}}^{- \frac{s}{2} + \frac{w}{2}} (- 1) d x + \int_{- \frac{w}{2}}^{\frac{w}{2}} 1 d x + \int_{\frac{s}{2} - \frac{w}{2}}^{\frac{s}{2}} (- 1) d x) \\ = & \frac{2}{s} (- (- \frac{s}{2} + \frac{w}{2} + \frac{s}{2}) + (\frac{w}{2} + \frac{w}{2}) - (\frac{s}{2} - \frac{s}{2} + \frac{w}{2})) \\ = & 0 \end{array}$

si96_e

$\begin{array}{l} a_{n} & = \frac{2}{s} (\int_{- \frac{s}{2}}^{\frac{s}{2}} f_{b a l . rect .} (x) cos (\frac{2 π n}{s} x) d x) \\ = \frac{2}{s} (\int_{- \frac{s}{2}}^{- \frac{s}{2} + \frac{w}{2}} (- cos (\frac{2 π n}{s} x)) d x + \int_{- \frac{w}{2}}^{\frac{w}{2}} cos (\frac{2 π n}{s} x) d x + \int_{\frac{s}{2} - \frac{w}{2}}^{\frac{s}{2}} (- cos (\frac{2 π n}{s} x))) d x \\ = \frac{2}{s} \frac{s}{2 π n} ({[- sin (\frac{2 π n}{s} x)]}_{- \frac{s}{2}}^{- \frac{s}{2} + \frac{w}{2}} + {[sin (\frac{2 π n}{s} x)]}_{- \frac{w}{2}}^{\frac{w}{2}} + \int - sin (\frac{2 π n}{s} x) \frac{s}{2} - \frac{w}{2} \frac{s}{2}) \\ = \frac{1}{π n} (- sin (\frac{2 π n}{s} (- \frac{s}{2} + \frac{w}{2})) + sin (\frac{2 π n}{s} (- \frac{s}{2})) + sin (\frac{2 π n}{s} \frac{w}{2}) - sin (\frac{2 π n}{s} (- \frac{w}{2})) - sin (\frac{2 π n}{s} \frac{s}{2}) + sin (\frac{2 π n}{s} (\frac{s}{2} - \frac{w}{2}))) \\ = \frac{1}{π n} (sin (\frac{2 π n}{s} (\frac{s}{2} - \frac{w}{2})) - sin (\frac{2 π n}{s} (\frac{s}{2})) + sin (\frac{2 π n}{s} \frac{w}{2}) + sin (\frac{2 π n}{s} (\frac{w}{2})) - sin (\frac{2 π n}{s} \frac{s}{2}) + sin (\frac{2 π n}{s} (\frac{s}{2} - \frac{w}{2}))) \\ = \frac{2}{π n} (sin (n π (1 - \frac{w}{2})) - sin (n π) + sin (n π \frac{w}{2})) \\ = \frac{2}{π n} (sin (n π (1 - \frac{w}{2})) + sin (n π \frac{w}{2})) \\ b_{n} & = \frac{2}{s} (\int_{- \frac{s}{2}}^{\frac{s}{2}} f_{b a l . rect . (x)} sin (\frac{2 π n}{s} x) d x) \\ = \frac{2}{s} (\int_{- \frac{s}{2}}^{- \frac{s}{2} + \frac{w}{2}} (- sin (\frac{2 π n}{s} x)) d x + \int_{- \frac{w}{2}}^{\frac{w}{2}} sin (\frac{2 π n}{s} x) d x + \int_{\frac{s}{2} - \frac{w}{2}}^{\frac{s}{2}} (- sin (\frac{2 π n}{s} x))) d x \\ = \frac{2}{s} \frac{s}{2 π n} ({[cos (\frac{2 π n}{s} x)]}_{- \frac{s}{2}}^{- \frac{s}{2} + \frac{w}{2}} + {[- cos (\frac{2 π n}{s} x)]}_{- \frac{w}{2}}^{\frac{w}{2}} + \int cos (\frac{2 π n}{s} x) \frac{s}{2} - \frac{w}{2} \frac{s}{2}) \\ = \frac{1}{π n} (cos (\frac{2 π n}{s} (- \frac{s}{2} + \frac{w}{2})) - cos (\frac{2 π n}{s} (- \frac{s}{2})) - cos (\frac{2 π n}{s} \frac{w}{2}) + cos (\frac{2 π n}{s} (- \frac{w}{s})) + cos (\frac{2 π n}{s} \frac{s}{2}) - cos (\frac{2 π n}{s} (\frac{s}{2} - \frac{w}{2}))) \\ = \frac{1}{π n} (cos (\frac{2 π n}{s} (\frac{s}{2} - \frac{w}{2})) - cos (\frac{2 π n}{s} (\frac{s}{2})) - cos (\frac{2 π n}{s} \frac{w}{2}) + cos (\frac{2 π n}{s} (\frac{w}{s})) + cos (\frac{2 π n}{s} \frac{s}{2}) - cos (\frac{2 π n}{s} (\frac{s}{2} - \frac{w}{2}))) \\ = 0 \end{array}$ $\begin{array}{l} a_{n} & = \frac{2}{s} (\int_{- \frac{s}{2}}^{\frac{s}{2}} f_{b a l . rect .} (x) cos (\frac{2 π n}{s} x) d x) \\ = \frac{2}{s} (\int_{- \frac{s}{2}}^{- \frac{s}{2} + \frac{w}{2}} (- cos (\frac{2 π n}{s} x)) d x + \int_{- \frac{w}{2}}^{\frac{w}{2}} cos (\frac{2 π n}{s} x) d x + \int_{\frac{s}{2} - \frac{w}{2}}^{\frac{s}{2}} (- cos (\frac{2 π n}{s} x))) d x \\ = \frac{2}{s} \frac{s}{2 π n} ({[- sin (\frac{2 π n}{s} x)]}_{- \frac{s}{2}}^{- \frac{s}{2} + \frac{w}{2}} + {[sin (\frac{2 π n}{s} x)]}_{- \frac{w}{2}}^{\frac{w}{2}} + \int - sin (\frac{2 π n}{s} x) \frac{s}{2} - \frac{w}{2} \frac{s}{2}) \\ = \frac{1}{π n} (- sin (\frac{2 π n}{s} (- \frac{s}{2} + \frac{w}{2})) + sin (\frac{2 π n}{s} (- \frac{s}{2})) + sin (\frac{2 π n}{s} \frac{w}{2}) - sin (\frac{2 π n}{s} (- \frac{w}{2})) - sin (\frac{2 π n}{s} \frac{s}{2}) + sin (\frac{2 π n}{s} (\frac{s}{2} - \frac{w}{2}))) \\ = \frac{1}{π n} (sin (\frac{2 π n}{s} (\frac{s}{2} - \frac{w}{2})) - sin (\frac{2 π n}{s} (\frac{s}{2})) + sin (\frac{2 π n}{s} \frac{w}{2}) + sin (\frac{2 π n}{s} (\frac{w}{2})) - sin (\frac{2 π n}{s} \frac{s}{2}) + sin (\frac{2 π n}{s} (\frac{s}{2} - \frac{w}{2}))) \\ = \frac{2}{π n} (sin (n π (1 - \frac{w}{2})) - sin (n π) + sin (n π \frac{w}{2})) \\ = \frac{2}{π n} (sin (n π (1 - \frac{w}{2})) + sin (n π \frac{w}{2})) \\ b_{n} & = \frac{2}{s} (\int_{- \frac{s}{2}}^{\frac{s}{2}} f_{b a l . rect . (x)} sin (\frac{2 π n}{s} x) d x) \\ = \frac{2}{s} (\int_{- \frac{s}{2}}^{- \frac{s}{2} + \frac{w}{2}} (- sin (\frac{2 π n}{s} x)) d x + \int_{- \frac{w}{2}}^{\frac{w}{2}} sin (\frac{2 π n}{s} x) d x + \int_{\frac{s}{2} - \frac{w}{2}}^{\frac{s}{2}} (- sin (\frac{2 π n}{s} x))) d x \\ = \frac{2}{s} \frac{s}{2 π n} ({[cos (\frac{2 π n}{s} x)]}_{- \frac{s}{2}}^{- \frac{s}{2} + \frac{w}{2}} + {[- cos (\frac{2 π n}{s} x)]}_{- \frac{w}{2}}^{\frac{w}{2}} + \int cos (\frac{2 π n}{s} x) \frac{s}{2} - \frac{w}{2} \frac{s}{2}) \\ = \frac{1}{π n} (cos (\frac{2 π n}{s} (- \frac{s}{2} + \frac{w}{2})) - cos (\frac{2 π n}{s} (- \frac{s}{2})) - cos (\frac{2 π n}{s} \frac{w}{2}) + cos (\frac{2 π n}{s} (- \frac{w}{s})) + cos (\frac{2 π n}{s} \frac{s}{2}) - cos (\frac{2 π n}{s} (\frac{s}{2} - \frac{w}{2}))) \\ = \frac{1}{π n} (cos (\frac{2 π n}{s} (\frac{s}{2} - \frac{w}{2})) - cos (\frac{2 π n}{s} (\frac{s}{2})) - cos (\frac{2 π n}{s} \frac{w}{2}) + cos (\frac{2 π n}{s} (\frac{w}{s})) + cos (\frac{2 π n}{s} \frac{s}{2}) - cos (\frac{2 π n}{s} (\frac{s}{2} - \frac{w}{2}))) \\ = 0 \end{array}$

si97_e

where we have exploited the symmetry properties, both of the sine and the cosine functions (see Eq. 3.16 and Eq. 3.19). We note a couple of things. First as we balance the total area of the function (we have two negative pulses each with only half of the width of the positive pulse), we end up with a total integral of 0 over the domain. This in turn will make a0 = 0. In addition, the pulses are arranged mirror-symmetrically to the origin, which results in bn = 0. We are therefore only left with the coefficients for the cosine term coefficients an. Upon closer inspection, we can simplify these even more as we find for all odd values of n, i.e., for 2n + 1

$sin ((2 n + 1) π (1 - \frac{w}{s})) = - sin ((2 n + 1) π (- \frac{w}{s})) = sin ((2 n + 1) π (- \frac{w}{s}))$ $sin ((2 n + 1) π (1 - \frac{w}{s})) = - sin ((2 n + 1) π (- \frac{w}{s})) = sin ((2 n + 1) π (- \frac{w}{s}))$

and therefore

$a_{n, odd} = \frac{4}{π (2 n + 1)} sin ((2 n + 1) π \frac{w}{s})$ $a_{n, odd} = \frac{4}{π (2 n + 1)} sin ((2 n + 1) π \frac{w}{s})$

On the other hand, we find for all even values of n, i.e., for 2n

$sin (2 n π (1 - \frac{w}{s})) = - sin (2 n π \frac{w}{s})$ $sin (2 n π (1 - \frac{w}{s})) = - sin (2 n π \frac{w}{s})$

and therefore

$a_{n, even} = 0$ $a_{n, even} = 0$

From which we can find the Fourier cosine series of the balanced rectangular function to be

$f_{b a l . rect .} (x) = \sum_{n = 0}^{\infty} \frac{4}{π (2 n + 1)} sin ((2 n + 1) π \frac{w}{s}) cos (2 π (2 n + 1) \frac{x}{s})$ $f_{b a l . rect .} (x) = \sum_{n = 0}^{\infty} \frac{4}{π (2 n + 1)} sin ((2 n + 1) π \frac{w}{s}) cos (2 π (2 n + 1) \frac{x}{s})$

si102_e (Eq. 19.35)

As you can see, the balanced rectangular function amounts to a very simple cosine Fourier series. Do not get confused because there is also a sine term in Eq. 19.35. If you look very carefully, you will see that this sine term is a constant, as it is not dependent on the independent variable x. The fact that we used a total of three pulses with matching orientation and integral allows us to drop the terms a₀ and b_n. We can now set arbitrary values for the central pulse width w (which is the only pulse we are interested in) and the spacing s of the two negative pulses. One thing to keep in mind is that the influence of the two negative pulses on the central pulse must be minimized. This is due to the fact that these pulses are only inserted to balance the rectangular function and make the problem analytically easier to solve. They have no physical meaning and are pure artifacts of the symmetry trick we used. Therefore we have to make sure to choose s sufficiently large in order to avoid all influence on the central pulse. Fig. 19.5 shows an example of the Fourier series of the balanced rectangular pulse with a central pulse width of w = 100 and a spacing s = 2000. The expansion order is n_max = 500. It is important to keep in mind that the higher the ratio $\frac{s}{w},$ $\frac{s}{w},$ the higher the expansion should be, otherwise the central pulse will not be a sharp rectangle.

f19-05-9781455731411 — Fig. 19.5 Fourier series of the balanced rectangular pulse function f_bai.rect. (x) expanded to an expansion order of nmax = 500. The pulse width was chosen as w = 100 the spacing as s = 2000.

Having derived the Fourier series for the balanced rectangular, we can now return to the Taylor-Aris dispersion problem outlined in Fig. 19.3. We derived the preliminary solution in Eq. 19.34 and can now turn to finding the missing constants c₃, c₄, and λ.

19.9.5 Applying the Boundary Condition

For t = 0 we find Eq. 19.34 to be

$\begin{array}{l} \bar{ρ} (0, z^{'}) & = c_{3} sin (\sqrt{λ} z^{'}) + c_{4} cos (\sqrt{λ} z^{'}) \\ \overset{!}{=} \sum_{n = 0}^{\infty} \frac{4}{π (2 n + 1)} (sin (2 n + 1) π \frac{2 R}{s}) cos (2 π (2 n + 1) \frac{z^{'}}{s}) \end{array}$ $\begin{array}{l} \bar{ρ} (0, z^{'}) & = c_{3} sin (\sqrt{λ} z^{'}) + c_{4} cos (\sqrt{λ} z^{'}) \\ \overset{!}{=} \sum_{n = 0}^{\infty} \frac{4}{π (2 n + 1)} (sin (2 n + 1) π \frac{2 R}{s}) cos (2 π (2 n + 1) \frac{z^{'}}{s}) \end{array}$

si104_e

where we have already applied Eq. 19.35, i.e., our solution to the initial rectangular concentration profile for t = 0. We inserted w = 2R and must only choose the spacing s to be sufficiently large. By comparison, we can now determine the missing constants to be

$\begin{array}{l} \sqrt{λ} = \frac{2 π (2 n + 1)}{s} \to λ = {(\frac{2 π (2 n + 1)}{s})}^{2} \\ c_{3} = 0 \\ c_{4} = \frac{4}{π (2 n + 1)} sin ((2 n + 1) π \frac{2 R}{s}) \end{array}$ $\begin{array}{l} \sqrt{λ} = \frac{2 π (2 n + 1)}{s} \to λ = {(\frac{2 π (2 n + 1)}{s})}^{2} \\ c_{3} = 0 \\ c_{4} = \frac{4}{π (2 n + 1)} sin ((2 n + 1) π \frac{2 R}{s}) \end{array}$

si105_e

in which case the solution to the Taylor-Aris dispersion problem can be written down from Eq. 19.34 as

$\begin{array}{c} \bar{ρ} (t, z^{'}) = \sum_{n = 0}^{\infty} e^{- {(\frac{2 π (2 n + 1)}{s})}^{2}} D_{eff .} t \frac{4}{π (2 n + 1)} sin ((2 n + 1) π \frac{2 R}{s}) cos (2 π (2 n + 1) \frac{z^{'}}{s}) \\ \bar{ρ} (t, z) = \sum_{n = 0}^{\infty} e^{- {(\frac{2 π (2 n + 1)}{s})}^{2}} D_{eff .} t \frac{4}{π (2 n + 1)} sin ((2 n + 1) π \frac{2 R}{s}) cos (2 π (2 n + 1) \frac{z - υ_{a v} t}{s}) \end{array}$ $\begin{array}{c} \bar{ρ} (t, z^{'}) = \sum_{n = 0}^{\infty} e^{- {(\frac{2 π (2 n + 1)}{s})}^{2}} D_{eff .} t \frac{4}{π (2 n + 1)} sin ((2 n + 1) π \frac{2 R}{s}) cos (2 π (2 n + 1) \frac{z^{'}}{s}) \\ \bar{ρ} (t, z) = \sum_{n = 0}^{\infty} e^{- {(\frac{2 π (2 n + 1)}{s})}^{2}} D_{eff .} t \frac{4}{π (2 n + 1)} sin ((2 n + 1) π \frac{2 R}{s}) cos (2 π (2 n + 1) \frac{z - υ_{a v} t}{s}) \end{array}$

si106_e (Eq. 19.36)

where we performed the back-transformation from z' to z according to Eq. 19.28.

19.9.6 Visualization

We will now visualize the profile as it moves along the detection channel as depicted in Fig. 19.3c/d. For this we will assume both the injection and the detection channel with a radius of 50 pm. We assume that the injected sample is a solution of ethanol in water, for which we find a diffusion constant of 1.24 × 10^–5 cm² s^− 1 according to Tab. 6.9. Obviously we expect a Hagen-Poiseuille flow profile (see section 16.3) along the detection channel with an average velocity v_av derived from Eq. 16.31 and Eq. 16.35 as

$_{υ_{a v} = - \frac{R^{2}}{8 η} \frac{d p}{d z}}$ $_{υ_{a v} = - \frac{R^{2}}{8 η} \frac{d p}{d z}}$

(Eq. 19.37)

We will calculate the Taylor-Aris dispersion for four different pressure drops along the z-axis: − 0.000 01 mbar mm^− 1, − 0.01 mbar mm^− 1, − 0.055 mbar mm^− 1, and − 0.1 mbarmm^− 1. For all cases we find Reynolds numbers below 1. For the last case (which is obviously the one with the highest velocity) the Reynolds number is Re = 0.156. The Peclet numbers are 0.001 26, 1.26, 6.93, and 12.6, respectively. As you can see from the increasing Peclet number, the convective mixing will increasingly contribute to the purely molecular diffusion, and we expect to see more effective mixing. This we can also see when plotting the effective diffusion constant D_eff. over the molecular diffusion D as a function of the pressure drop (see Fig. 19.6). As you can see, we expect to see a fourfold increase of the effective diffusion at the highest pressure drop in our example. Fig. 19.7 shows the profiles we find along the detection channel at different points in time. As can be seen, the concentration profiles broaden faster at higher Peclet numbers, i.e., if diffusive mixing is supported by convective mixing. As an example, compare the overall height and width of the concentration profile for time point t= 0.2 s at a pressure drop of − 0.00001 mbarmm^− 1 with the profile obtained for a pressure drop of − 0.1 mbarmm^− 1. The first profile reaches an overall relative height of about 0.5, whereas the latter reaches a height of only 0.3 as a consequence of the increased diffusion.

f19-06-9781455731411 — Fig. 19.6 Effective diffusion D_eff. as a function of the pressure drop in Taylor-Aris dispersion. As you can see, the effective diffusion increases substantially when taking the convective mixing into account.

f19-07-9781455731411 — Fig. 19.7 Taylor-Aris dispersion at a channel intersection with circular cross-sections and a radius of R = 50μm and a spacing value of s = 5000 μm. The profiles show the concentration profile $\bar{p}$ $\bar{p}$ (z) of an ethanol plug as it moves along the detection channel at different points in time 0 s, 0.001 s, 0.01 s, and 0.1 s. As can be seen, the concentration profiles flatten out significantly slower at low Péclet numbers, i.e., when the contribution from convective mixing is small. As the average velocity v_avincreases the mixing increases and the plug profiles flatten out faster.

19.9.7 Maple Worksheet

The plots in Fig. 19.7 were created using the Maple worksheet shown in listing 19.1 which we will briefly discuss. In line 1 we include Core.txt and define the velocity velzav in line 4 according to Eq. 19.37 and the effective diffusion Deff according to Eq. 19.24. The solution function to the Taylor-Aris dispersion is defined in line 10 which is given by Eq. 19.36. We define this function as a function of time t, z, radius R and the average velocity v_av, as well as the expansion order of the Fourier series n_max and the spacing variable s. In line 13 we correct the function for unit, expecting the independent variables time t in s, z and radius R in pm, diffusion coefficient in 10^–5 cm² s^− 1, viscosity η in mPas, and pressure drop in mbar mm^− 1. Furthermore we set the expansion order n_max = 200 and the spacing s = 5000 μm. As already discussed, it is necessary to set this length to a high value

Listing 19.1

[Maple] Listing for Taylor-Aris dispersion. A digital version of the listing can be found under the name TaylorAris.mw in the supplementary material.

1 restart:read "Core.txt":
2
3 #this is the average velocity along the detection channel
4 velzav:=(R,eta,dpdz)->(− R~2/(8*eta)*dpdz);
5
6 #this calculates the effective diffusion coefficient
7 Deff:=(D,R,velzav)->D+R~2*velzav~2/(48*D);
8
9 #Taylor-Aris dispersion
10 CmAvNorm:=(t,z,R,D,velzav,nmax,s)->sum(4/(Pi*(2*n+1))*sin(Pi*(2*n+1)*(2*R)/s)*exp(−(2*Pi*(2*n+1)/s) ~ 2*Deff(D,R,velzav)*t)*cos(2*Pi*(2*n+1)/s*(z - velzav*t)),n=0..nmax);
11
12 #with units
13 CmAvNorm_units:=unapply(evalf(simplify(CmAvNorm(_t*Unit(s),_z*Unit(um),_R*Unit(um),_D*10E-5*Unit(cm~2/s),velzav(_R*Unit(um),_eta*Unit(mPa*s),_dpdz*Unit(mbar/mm)),200,5000*Unit(um)))),_t,_z,_R ,_dpdz,_eta,_D):
14
15 #Taylor-Aris diffusion of ethanol in water
16 quickplot([z->CmAvNorm_units(0,z,50,-0.00001,1,1.24),z->CmAvNorm_units(0.01,z,50,-0.01,1,1.24),z-> CmAvNorm_units(0.1,z,50,-0.00001,1,1.24),z->CmAvNorm_units(0.2,z,50,-0.00001,1,1.24) ],-200..1000,"z [um]",typeset(conjugate(rho)/conjugate(rho)[0]~empty),["t=0 s","t=0.01 s","t=0.1 s","t=0.2 s"]);
17 quickplot([z->CmAvNorm_units(0,z,50,-0.01,1,1.24),z->CmAvNorm_units(0.01,z,50,-0.01,1,1.24),z-> CmAvNorm_units(0.1,z,50,-0.01,1,1.24),z->CmAvNorm_units(0.2,z,50,-0.01,1,1.24)],-200..1000,"z [um]",typeset(conjugate(rho)/conjugate(rho)[0]~empty),["t=0 s","t=0.01 s","t=0.1 s","t=0.2 s"]);
18 quickplot([z->CmAvNorm_units(0,z,50,-0.055,1,1.24),z->CmAvNorm_units(0.01,z,50,-0.055,1,1.24),z-> CmAvNorm_units(0.1,z,50,-0.055,1,1.24),z->CmAvNorm_units(0.2,z,50,-0.055,1,1.24)],-200..1000," z [um]",typeset(conjugate(rho)/conjugate(rho)[0]~empty),["t=0 s","t=0.01 s","t=0.1 s","t=0.2 s"]);
19 quickplot([z->CmAvNorm_units(0,z,50,-0.1,1,1.24),z->CmAvNorm_units(0.01,z,50,-0.1,1,1.24),z-> CmAvNorm_units(0.1,z,50,-0.1,1,1.24),z->CmAvNorm_units(0.2,z,50,-0.1,1,1.24)],-200..1000,"z [ um]",typeset(conjugate(rho)/conjugate(rho)[0]~empty),["t=0 s","t=0.01 s","t=0.1 s","t=0.2 s"]);

in order to avoid all artifacts resulting from the balancing of the rectangular function to interfere with the effect on the central pulse. Using the function unapply, we turn the series into a function of the independent variables time t, z, radius R, pressure drop, viscosity η, and diffusion coefficient D. Line 16 calculates the profile for a pressure drop of − 0.00001 mbarmm^− 1, line 17 for a pressure drop of − 0.01 mbarmm^− 1, line 18 for a pressure drop of − 0.055 mbar mm^− 1, and line 19 for a pressure drop of 0.1 mbar mm^− 1. The output of these lines are the plots shown in Fig. 19.7.

Using this worksheet, it is very easy to adapt the Taylor-Aris dispersion to different substances and flow conditions.

19.9.8 Comment on the Timescale

In section 19.6.3 we calculated that the timescale on which we can neglect radial diffusion is in the order of 0.4 ps. As you can see in our example, we are significantly above these timescales, therefore the assumption that we can neglect radial diffusion is fully realistic.

19.9.9 Reducing Taylor-Aris Dispersion

In general, dispersion cannot be avoided. However, it can be reduced by limiting the lengths over which a liquid diffusing sample has to be transported [3]. In addition, several approaches have been described in the literature which make use of, e.g., the insertion of diffusion barriers such as immiscible oil plugs. This concept has been explored, e.g., in indirect microfluidic systems which make use of tetradecane separators to completely block dispersion effects [4].

19.10 SUMMARY

Taylor-Aris diffusion is one of the most important and commonly encountered phenomena in microfluidics. As discussed, mixing in microfluidic systems is often problematic due to the low Reynolds numbers and therefore, the low convective contributions to the overall mixing process. However, given a certain channel cross-section, one can see that convection can increase the diffusive mixing substantially at higher Péclet numbers. These can be achieved in microfluidic channel geometries, thereby making dispersion an effect to consider. As we have seen, the mathematical derivation of the Taylor-Aris dispersion gives rise to a slightly modified diffusion equation with an effective diffusion coefficient which depends on the flow conditions. The most important parameter contributing to this effective diffusion is the average velocity. At higher flow rates, i.e., higher pressure drops the effective diffusion can be increased several-fold resulting in effective mixing.

..................Content has been hidden....................

You can't read the all page of ebook, please click here login for view all page.

Table of Contents for Chapter 19: Taylor-Aris Dispersion

Create new playlist

Sign In

Sign Up