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 17

Hydraulic Resistance

17.1 Introduction

In microfluidics fluids are probed through channels with small cross-sections. Depending on the material the channel is made of the initial filling may either be facilitated or hindered by capillary effects. This is dependent on the wetting angle of the fluid on the material of the microchannel (see section 21). However, once the channel is filled with the fluid forces must be applied onto the fluid in order to move it. We know intuitively that it requires a certain amount of pressure in order to squeeze a liquid through a tiny orifice such as a microfluidic channel. A good example are squirt guns where the distance the water will travel depends on the velocity of the liquid upon exiting the gun. In order to make the gun shoot farther we need to apply more pressure in order to accelerate the water further. However, the effective width that the squirt gun can shoot additionally depends on the size of the orifice through which the water is squeezed. If a smaller orifice is used (some squirt guns allow changing the nozzle at the tip of the gun) the water will exit at higher speeds, as we can see from Eq. 10.8, and thus travel farther. However, we will also find that we need significantly higher pressures in order to squeeze the water through the nozzle.

This effect is due to the internal friction of the fluid. Internal friction will cause the fluid to lose kinetic energy. This energy is converted into thermal energy, therefore heating up the fluid. One would expect that higher flow rates also result in higher internal friction. As we will see in a moment this intuitive correlation is not correct as fluid friction increases with increasing gradients in velocity within the fluid. Obviously if the flow rates in, e.g., a circular channel are higher, we will also have higher gradients in the velocity profile if we assume a no-slip condition at the channel walls.

The conversion of kinetic energy into thermal energy is ubiquitous in fluid systems. This effect is generally referred to as viscous dissipation. Technically the terms pressure loss and pressure drop are often used to refer to the same effect. The fluid flowing through a system at a given flowrate will lose kinetic energy due to viscous dissipation that results in a reduction of the pressure according to the Bernoulli equation, Eq. 14.11. This pressure drop is dependent on the flowrate and the geometry of the fluidic system.

In the following, we will derive the equations for quantifying viscous dissipation.

17.2 Viscous Dissipation

As stated, viscous dissipation will reduce the kinetic energy of a fluid over time. Obviously, the total kinetic energy of a moving control volume is given by

Ekin=12mv→2 $E_{k i n} = \frac{1}{2} m {\vec{v}}^{2}$

Unfortunately, the velocity is not uniform across our control volume, which is why we need to use an integral notation. This integral needs to be formed over the entire domain Ω, which is given by

Ekin==∫Ω12ρv→2dV12∫Ωρv→2dV $\begin{array}{l} E_{k i n} & = & \int_{Ω} \frac{1}{2} ρ {\vec{v}}^{2} d V \\ = & \frac{1}{2} \int_{Ω} ρ {\vec{v}}^{2} d V \end{array}$

si8_e

Change of Kinetic Energy. The kinetic energy will change over time as a consequence of the viscous dissipation, i.e., the transformation of kinetic energy to heat as a consequence of the friction inside of the fluid. The change of energy over time is described by

∂Ekin∂t====∂∂t(12∫Ωρv→2dV)12∫Ωρ∂v→2∂tdV12∫Ωρ2v→∂v→∂tdV∫Ωρ∂v→∂tv→dV $\begin{array}{l} \frac{\partial E_{kin}}{\partial t} & = & \frac{\partial}{\partial t} (\frac{1}{2} \int_{Ω} ρ {\vec{v}}^{2} d V) \\ = & \frac{1}{2} \int_{Ω} ρ \frac{\partial {\vec{v}}^{2}}{\partial t} d V \\ = & \frac{1}{2} \int_{Ω} ρ 2 \vec{v} \frac{\partial \vec{v}}{\partial t} d V \\ = & \int_{Ω} ρ \frac{\partial \vec{v}}{\partial t} \vec{v} d V \end{array}$

si9_e (Eq. 17.1)

where we have used the fact that we can interchange the integral and the differential and applied the chain rule of differentiation (see Eq. 3.14).

Navier-Stokes Equation. We now need a way to evaluate the integral given by Eq. 17.1. If we restrict our derivation to incompressible Newtonian fluids we can use the Navier-Stokes equation (see Eq. 11.40), which we can rewrite to:

ρ∂v→∂t=k→−∇→p+ηΔv→−ρ(v→⋅∇→)v→ $ρ \frac{\partial \vec{v}}{\partial t} = \vec{k} - \vec{\nabla} p + η Δ \vec{v} - ρ (\vec{v} \cdot \vec{\nabla}) \vec{v}$

si10_e (Eq. 17.2)

which we insert into Eq. 17.1 to find

∂Ekin∂t==ρ∫Ωv→(k→−∇→p+ηΔv→−ρ(v→⋅∇→)v→)dVρ∫Ωv→(k→−∇→p)dV+ρ∫Ωv→(ηΔv→−(v→⋅∇→)v→)dV $\begin{array}{l} \frac{\partial E_{k i n}}{\partial t} & = & ρ \int_{Ω} \vec{v} (\vec{k} - \vec{\nabla} p + η Δ \vec{v} - ρ (\vec{v} \cdot \vec{\nabla}) \vec{v}) d V \\ = & ρ \int_{Ω} \vec{v} (\vec{k} - \vec{\nabla} p) d V + ρ \int_{Ω} \vec{v} (η Δ \vec{v} - (\vec{v} \cdot \vec{\nabla}) \vec{v}) d V \end{array}$

si11_e (Eq. 17.3)

In the following, we will transform the individual terms of the integrals in Eq. 17.3 into gradient forms. We will do this because we want to transform this volume integral into a surface integral using Gauss’s theorem. We will start with the first integral and concentrate on the pressure term. The volume force term cannot be converted into a gradient because it has no partial derivatives to it.

Converting the Pressure Term Into a Gradient. The pressure term of the first integral of Eq. 17.3 is given by

−∫Ωv→⋅(∇→p)dV $- \int_{Ω} \vec{v} \cdot (\vec{\nabla} p) d V$

si12_e

which can be converted from vector notation to read

−v→⋅(∇→p)=−(vx∂p∂x+vy∂p∂y+vz∂p∂z) $- \vec{v} \cdot (\vec{\nabla} p) = - (v_{x} \frac{\partial p}{\partial x} + v_{y} \frac{\partial p}{\partial y} + v_{z} \frac{\partial p}{\partial z})$

si13_e (Eq. 17.4)

Applying the chain rule of differentiation (see Eq. 3.14) we note the following equivalences for v_x:

∂(vxp)∂xvx∂p∂x==∂vx∂xp+vx∂p∂x∂(vxp)∂x−∂vx∂xp $\begin{array}{l} \frac{\partial (v_{x} p)}{\partial x} & = & \frac{\partial v_{x}}{\partial x} p + v_{x} \frac{\partial p}{\partial x} \\ v_{x} \frac{\partial p}{\partial x} & = & \frac{\partial (v_{x} p)}{\partial x} - \frac{\partial v_{x}}{\partial x} p \end{array}$

si14_e (Eq. 17.5)

for v_y:

∂(vyp)∂yvy∂p∂y==∂vy∂yp+vy∂p∂y∂(vyp)∂y−∂vy∂yp $\begin{array}{l} \frac{\partial (v_{y} p)}{\partial y} & = & \frac{\partial v_{y}}{\partial y} p + v_{y} \frac{\partial p}{\partial y} \\ v_{y} \frac{\partial p}{\partial y} & = & \frac{\partial (v_{y} p)}{\partial y} - \frac{\partial v_{y}}{\partial y} p \end{array}$

si15_e (Eq. 17.6)

for v_z:

∂(vzp)∂zvz∂p∂z==∂vz∂zp+vz∂p∂z∂(vzp)∂z−∂vz∂zp $\begin{array}{l} \frac{\partial (v_{z} p)}{\partial z} & = & \frac{\partial v_{z}}{\partial z} p + v_{z} \frac{\partial p}{\partial z} \\ v_{z} \frac{\partial p}{\partial z} & = & \frac{\partial (v_{z} p)}{\partial z} - \frac{\partial v_{z}}{\partial z} p \end{array}$

si16_e (Eq. 17.7)

Inserting Eq. 17.5, Eq. 17.6, and Eq. 17.7 into Eq. 17.4 yields

v→⋅(∇→p)===−∂(vxp)∂x−∂vx∂xp+∂(vyp)∂y−∂vy∂yp+∂(vzp)∂z−∂vz∂zp(−∇→⋅(v→p)−p(∂vx∂x+∂vy∂y+∂vz∂z))−∇→⋅(v→p) $\begin{array}{l} \vec{v} \cdot (\vec{\nabla} p) & = & - \frac{\partial (v_{x} p)}{\partial x} - \frac{\partial v_{x}}{\partial x} p + \frac{\partial (v_{y} p)}{\partial y} - \frac{\partial v_{y}}{\partial y} p + \frac{\partial (v_{z} p)}{\partial z} - \frac{\partial v_{z}}{\partial z} p \\ = & (- \vec{\nabla} \cdot (\vec{v} p) - p (\frac{\partial v_{x}}{\partial x} + \frac{\partial v_{y}}{\partial y} + \frac{\partial v_{z}}{\partial z})) \\ = & - \vec{\nabla} \cdot (\vec{v} p) \end{array}$

si17_e (Eq. 17.8)

where we have used the continuity equation (see Eq. 10.17) to set

∂vx∂x+∂vy∂y+∂vz∂z=∇→⋅v→=0 $\frac{\partial v_{x}}{\partial x} + \frac{\partial v_{y}}{\partial y} + \frac{\partial v_{z}}{\partial z} = \vec{\nabla} \cdot \vec{v} = 0$

We have now obtained a gradient notation for the pressure term. Next, we will discuss the convection term. Converting the Convection Term Into a Gradient. The second integral of Eq. 17.3 is

∫Ω(v→ηΔv→−v→(v→⋅∇→)v→)dV $\int_{Ω} (\vec{v} η Δ \vec{v} - \vec{v} (\vec{v} \cdot \vec{\nabla}) \vec{v}) d V$

si19_e (Eq. 17.9)

We will start by looking at the second term in this integral, which is the term due to convection. Reconverting from the vector notation we obtain

==v→(v→⋅∇→)v→v→⎛⎝⎜⎜⎜⎜vx∂vx∂x+vy∂vx∂y+vz∂vx∂zvx∂vy∂x+vy∂vy∂y+vz∂vy∂zvx∂vz∂x+vy∂vz∂y+vz∂vz∂z⎞⎠⎟⎟⎟⎟vxvx∂vx∂x+vyvx∂vx∂y+vzvx∂vx∂z+vxvy∂vy∂x+vyvy∂vy∂y+vzvy∂vy∂z+vxvz∂vz∂x+vyvz∂vz∂y+vzvz∂vz∂z $\begin{array}{l} \vec{v} (\vec{v} \cdot \vec{\nabla}) \vec{v} \\ = & \vec{v} (\begin{matrix} v_{x} \frac{\partial v_{x}}{\partial x} + v_{y} \frac{\partial v_{x}}{\partial y} + v_{z} \frac{\partial v_{x}}{\partial z} \\ v_{x} \frac{\partial v_{y}}{\partial x} + v_{y} \frac{\partial v_{y}}{\partial y} + v_{z} \frac{\partial v_{y}}{\partial z} \\ v_{x} \frac{\partial v_{z}}{\partial x} + v_{y} \frac{\partial v_{z}}{\partial y} + v_{z} \frac{\partial v_{z}}{\partial z} \end{matrix}) \\ = & v_{x} v_{x} \frac{\partial v_{x}}{\partial x} + v_{y} v_{x} \frac{\partial v_{x}}{\partial y} + v_{z} v_{x} \frac{\partial v_{x}}{\partial z} + v_{x} v_{y} \frac{\partial v_{y}}{\partial x} + v_{y} v_{y} \frac{\partial v_{y}}{\partial y} + v_{z} v_{y} \frac{\partial v_{y}}{\partial z} + v_{x} v_{z} \frac{\partial v_{z}}{\partial x} + v_{y} v_{z} \frac{\partial v_{z}}{\partial y} + v_{z} v_{z} \frac{\partial v_{z}}{\partial z} \end{array}$

si20_e (Eq. 17.10)

Applying the chain rule of differentiation (see Eq. 3.14) we note the following equivalences for v_x:

12∂v2x∂x12∂v2x∂y12∂v2x∂z===vx∂vx∂xvx∂vx∂yvx∂vx∂z $\begin{array}{l} \frac{1}{2} \frac{\partial v_{x}^{2}}{\partial x} & = & v_{x} \frac{\partial v_{x}}{\partial x} \\ \frac{1}{2} \frac{\partial v_{x}^{2}}{\partial y} & = & v_{x} \frac{\partial v_{x}}{\partial y} \\ \frac{1}{2} \frac{\partial v_{x}^{2}}{\partial z} & = & v_{x} \frac{\partial v_{x}}{\partial z} \end{array}$

si21_e

for v_y:

12∂v2y∂x12∂v2y∂y12∂v2y∂z===vy∂vy∂xvy∂vy∂yvy∂vy∂z $\begin{array}{l} \frac{1}{2} \frac{\partial v_{y}^{2}}{\partial x} & = & v_{y} \frac{\partial v_{y}}{\partial x} \\ \frac{1}{2} \frac{\partial v_{y}^{2}}{\partial y} & = & v_{y} \frac{\partial v_{y}}{\partial y} \\ \frac{1}{2} \frac{\partial v_{y}^{2}}{\partial z} & = & v_{y} \frac{\partial v_{y}}{\partial z} \end{array}$

si22_e

for v_z:

12∂v2z∂x12∂v2z∂y12∂v2z∂z===vz∂vz∂xvz∂vz∂yvz∂vz∂z $\begin{array}{l} \frac{1}{2} \frac{\partial v_{z}^{2}}{\partial x} & = & v_{z} \frac{\partial v_{z}}{\partial x} \\ \frac{1}{2} \frac{\partial v_{z}^{2}}{\partial y} & = & v_{z} \frac{\partial v_{z}}{\partial y} \\ \frac{1}{2} \frac{\partial v_{z}^{2}}{\partial z} & = & v_{z} \frac{\partial v_{z}}{\partial z} \end{array}$

si23_e

which we reinsert into Eq. 17.10 to find

v→(v→⋅∇→)v→==12(vx∂v2x∂x+vy∂v2x∂yvz∂v2x∂z+vx∂v2y∂x+vy∂v2y∂yvz∂v2y∂z+vx∂v2z∂x+vy∂v2z∂yvz∂v2z∂z)12(vx∂∂x(v2x+v2y+v2z)+vy∂∂y(v2x+v2y+v2z)+vz∂∂z(v2x+v2y+v2z)) $\begin{array}{l} \vec{v} (\vec{v} \cdot \vec{\nabla}) \vec{v} & = & \frac{1}{2} (v_{x} \frac{\partial v_{x}^{2}}{\partial x} + v_{y} \frac{\partial v_{x}^{2}}{\partial y} v_{z} \frac{\partial v_{x}^{2}}{\partial z} + v_{x} \frac{\partial v_{y}^{2}}{\partial x} + v_{y} \frac{\partial v_{y}^{2}}{\partial y} v_{z} \frac{\partial v_{y}^{2}}{\partial z} + v_{x} \frac{\partial v_{z}^{2}}{\partial x} + v_{y} \frac{\partial v_{z}^{2}}{\partial y} v_{z} \frac{\partial v_{z}^{2}}{\partial z}) \\ = & \frac{1}{2} (v_{x} \frac{\partial}{\partial x} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) + v_{y} \frac{\partial}{\partial y} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) + v_{z} \frac{\partial}{\partial z} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2})) \end{array}$

si24_e (Eq. 17.11)

We now do something seemingly strange, adding the following term to Eq. 17.11:

(v2x+v2y+v2z)(∂vx∂x+∂vy∂y+∂vz∂z)=(v2x+v2y+v2z)(∇→⋅v→)=0 $(v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) (\frac{\partial v_{x}}{\partial x} + \frac{\partial v_{y}}{\partial y} + \frac{\partial v_{z}}{\partial z}) = (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) (\vec{\nabla} \cdot \vec{v}) = 0$

si25_e

Obviously, the term ∇→⋅v→ $\vec{\nabla} \cdot \vec{v}$ is equal to 0 due to the continuity equation (see Eq. 10.17). But why would be want to add a term that is effectively 0 to Eq. 17.11? We do this because of the following equalities:

for v_x

∂∂x(vx(v2x+v2y+v2z))vx∂∂x(v2x+v2y+v2z)==∂vx∂x(v2x+v2y+v2z)+vx∂∂x(v2x+v2y+v2z)∂∂x(vx(v2x+v2y+v2z))−∂vx∂x(v2x+v2y+v2z) $\begin{array}{l} \frac{\partial}{\partial x} (v_{x} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2})) & = & \frac{\partial v_{x}}{\partial x} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) + v_{x} \frac{\partial}{\partial x} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) \\ v_{x} \frac{\partial}{\partial x} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) & = & \frac{\partial}{\partial x} (v_{x} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2})) - \frac{\partial v_{x}}{\partial x} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) \end{array}$

si27_e (Eq. 17.12)

for v_y

∂∂y(vy(v2x+v2y+v2z))vy∂∂y(v2x+v2y+v2z)==∂vy∂y(v2x+v2y+v2z)+vy∂∂y(v2x+v2y+v2z)∂∂y(vx(v2x+v2y+v2z))−∂vy∂y(v2x+v2y+v2z) $\begin{array}{l} \frac{\partial}{\partial y} (v_{y} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2})) & = & \frac{\partial v_{y}}{\partial y} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) + v_{y} \frac{\partial}{\partial y} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) \\ v_{y} \frac{\partial}{\partial y} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) & = & \frac{\partial}{\partial y} (v_{x} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2})) - \frac{\partial v_{y}}{\partial y} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) \end{array}$

si28_e (Eq. 17.13)

for v_z

∂∂z(vx(v2x+v2y+v2z))vz∂∂z(v2x+v2y+v2z)==∂vz∂z(v2x+v2y+v2z)+vz∂∂z(v2x+v2y+v2z)∂∂z(vx(v2x+v2y+v2z))−∂vz∂z(v2x+v2y+v2z) $\begin{array}{l} \frac{\partial}{\partial z} (v_{x} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2})) & = & \frac{\partial v_{z}}{\partial z} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) + v_{z} \frac{\partial}{\partial z} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) \\ v_{z} \frac{\partial}{\partial z} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) & = & \frac{\partial}{\partial z} (v_{x} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2})) - \frac{\partial v_{z}}{\partial z} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) \end{array}$

si29_e (Eq. 17.14)

Inserting Eq. 17.12, Eq. 17.13, and Eq. 17.14 into Eq. 17.11 yields

v→(v→⋅∇→)v→===12(∂∂x(vx(v2x+v2y+v2z))−∂vx∂x(v2x+v2y+v2z)+∂∂y(vy(v2x+v2y+v2z))−∂vy∂y(v2x+v2y+v2z)+∂∂z(vx(v2x+v2y+v2z))−∂vz∂z(v2x+v2y+v2z))12(∇→⋅v→(v2x+v2y+v2z)−(∂vx∂x+∂vy∂y+∂vz∂z)(v2x+v2y+v2z))12∇→⋅v→(v2x+v2y+v2z) $\begin{array}{l} \vec{v} (\vec{v} \cdot \vec{\nabla}) \vec{v} & = & \frac{1}{2} (\frac{\partial}{\partial x} (v_{x} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2})) - \frac{\partial v_{x}}{\partial x} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) + \frac{\partial}{\partial y} (v_{y} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2})) - \frac{\partial v_{y}}{\partial y} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) + \frac{\partial}{\partial z} (v_{x} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2})) - \frac{\partial v_{z}}{\partial z} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2})) \\ = & \frac{1}{2} (\vec{\nabla} \cdot \vec{v} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) - (\frac{\partial v_{x}}{\partial x} + \frac{\partial v_{y}}{\partial y} + \frac{\partial v_{z}}{\partial z}) (v_{x}^{2} + v_{y}^{2} + v_{z}^{2})) \\ = & \frac{1}{2} \vec{\nabla} \cdot \vec{v} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) \end{array}$

si30_e (Eq. 17.15)

where we have again used the continuity equation (see Eq. 10.17) to set

∂vx∂x+∂vy∂y+∂vz∂z=∇→⋅v→=0 $\frac{\partial v_{x}}{\partial x} + \frac{\partial v_{y}}{\partial y} + \frac{\partial v_{z}}{\partial z} = \vec{\nabla} \cdot \vec{v} = 0$

(Eq. 17.16)

We now also have a gradient notation for the convection term. It only remains to find a gradient notation for the viscosity term.

Converting the Viscosity Term Into a Gradient. The first term in the integral of Eq. 17.9 is due to the viscous momentum transport. As it turns out, this is the term we are most interested in. We first need to convert this term from vector notation resulting in

v→ηΔv→==η(vx(∂2vx∂x2+∂2vx∂y2+∂2vx∂z2)+vy(∂2vy∂x2+∂2vy∂y2+∂2vy∂z2)+vz(∂2vz∂x2+∂2vz∂y2+∂2vz∂z2))η(vx∂2vx∂x2+vx∂2vx∂y2+vx∂2vx∂z2+vy∂2vy∂x2+vy∂2vy∂y2+vy∂2vy∂z2+vz∂2vz∂x2+vz∂2vz∂y2+vz∂2vz∂z2) $\begin{array}{l} \vec{v} η Δ \vec{v} & = & η (v_{x} (\frac{\partial^{2} v_{x}}{\partial x^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}}) + v_{y} (\frac{\partial^{2} v_{y}}{\partial x^{2}} + \frac{\partial^{2} v_{y}}{\partial y^{2}} + \frac{\partial^{2} v_{y}}{\partial z^{2}}) + v_{z} (\frac{\partial^{2} v_{z}}{\partial x^{2}} + \frac{\partial^{2} v_{z}}{\partial y^{2}} + \frac{\partial^{2} v_{z}}{\partial z^{2}})) \\ = & η (v_{x} \frac{\partial^{2} v_{x}}{\partial x^{2}} + v_{x} \frac{\partial^{2} v_{x}}{\partial y^{2}} + v_{x} \frac{\partial^{2} v_{x}}{\partial z^{2}} + v_{y} \frac{\partial^{2} v_{y}}{\partial x^{2}} + v_{y} \frac{\partial^{2} v_{y}}{\partial y^{2}} + v_{y} \frac{\partial^{2} v_{y}}{\partial z^{2}} + v_{z} \frac{\partial^{2} v_{z}}{\partial x^{2}} + v_{z} \frac{\partial^{2} v_{z}}{\partial y^{2}} + v_{z} \frac{\partial^{2} v_{z}}{\partial z^{2}}) \end{array}$

si32_e (Eq. 17.16)

Applying the chain rule of differentiation (see Eq. 3.14) we note the following equalities

for v_x

∂∂x(vx∂vx∂x)∂∂y(vx∂vx∂y)∂∂z(vx∂vx∂z)===(∂vx∂x)2+vx∂2vx∂x2→vx∂2vx∂x2=∂∂x(vx∂vx∂x)−(∂vx∂x)2(∂vx∂y)2+vx∂2vx∂y2→vx∂2vx∂y2=∂∂y(vx∂vx∂y)−(∂vx∂y)2(∂vx∂z)2+vx∂2vx∂z2→vx∂2vx∂z2=∂∂z(vx∂vx∂z)−(∂vx∂z)2 $\begin{array}{l} \frac{\partial}{\partial x} (v_{x} \frac{\partial v_{x}}{\partial x}) & = & {(\frac{\partial v_{x}}{\partial x})}^{2} + v_{x} \frac{\partial^{2} v_{x}}{\partial x^{2}} \to v_{x} \frac{\partial^{2} v_{x}}{\partial x^{2}} = \frac{\partial}{\partial x} (v_{x} \frac{\partial v_{x}}{\partial x}) - {(\frac{\partial v_{x}}{\partial x})}^{2} \\ \frac{\partial}{\partial y} (v_{x} \frac{\partial v_{x}}{\partial y}) & = & {(\frac{\partial v_{x}}{\partial y})}^{2} + v_{x} \frac{\partial^{2} v_{x}}{\partial y^{2}} \to v_{x} \frac{\partial^{2} v_{x}}{\partial y^{2}} = \frac{\partial}{\partial y} (v_{x} \frac{\partial v_{x}}{\partial y}) - {(\frac{\partial v_{x}}{\partial y})}^{2} \\ \frac{\partial}{\partial z} (v_{x} \frac{\partial v_{x}}{\partial z}) & = & {(\frac{\partial v_{x}}{\partial z})}^{2} + v_{x} \frac{\partial^{2} v_{x}}{\partial z^{2}} \to v_{x} \frac{\partial^{2} v_{x}}{\partial z^{2}} = \frac{\partial}{\partial z} (v_{x} \frac{\partial v_{x}}{\partial z}) - {(\frac{\partial v_{x}}{\partial z})}^{2} \end{array}$

si33_e

for v_y

∂∂x(vy∂vy∂x)∂∂y(vy∂vy∂y)∂∂z(vy∂vy∂z)===(∂vy∂x)2+vy∂2vy∂x2→vy∂2vy∂x2=∂∂x(vy∂vy∂x)−(∂vy∂x)2(∂vy∂y)2+vy∂2vy∂y2→vy∂2vy∂y2=∂∂y(vy∂vy∂y)−(∂vy∂y)2(∂vy∂z)2+vy∂2vy∂z2→vy∂2vy∂z2=∂∂z(vy∂vy∂z)−(∂vy∂z)2 $\begin{array}{l} \frac{\partial}{\partial x} (v_{y} \frac{\partial v_{y}}{\partial x}) & = & {(\frac{\partial v_{y}}{\partial x})}^{2} + v_{y} \frac{\partial^{2} v_{y}}{\partial x^{2}} \to v_{y} \frac{\partial^{2} v_{y}}{\partial x^{2}} = \frac{\partial}{\partial x} (v_{y} \frac{\partial v_{y}}{\partial x}) - {(\frac{\partial v_{y}}{\partial x})}^{2} \\ \frac{\partial}{\partial y} (v_{y} \frac{\partial v_{y}}{\partial y}) & = & {(\frac{\partial v_{y}}{\partial y})}^{2} + v_{y} \frac{\partial^{2} v_{y}}{\partial y^{2}} \to v_{y} \frac{\partial^{2} v_{y}}{\partial y^{2}} = \frac{\partial}{\partial y} (v_{y} \frac{\partial v_{y}}{\partial y}) - {(\frac{\partial v_{y}}{\partial y})}^{2} \\ \frac{\partial}{\partial z} (v_{y} \frac{\partial v_{y}}{\partial z}) & = & {(\frac{\partial v_{y}}{\partial z})}^{2} + v_{y} \frac{\partial^{2} v_{y}}{\partial z^{2}} \to v_{y} \frac{\partial^{2} v_{y}}{\partial z^{2}} = \frac{\partial}{\partial z} (v_{y} \frac{\partial v_{y}}{\partial z}) - {(\frac{\partial v_{y}}{\partial z})}^{2} \end{array}$

si34_e

for v_z

∂∂x(vz∂vz∂x)∂∂y(vz∂vz∂y)∂∂z(vz∂vz∂z)===(∂vz∂x)2+vz∂2vz∂x2→vz∂2vz∂x2=∂∂x(vz∂vz∂x)−(∂vz∂x)2(∂vz∂y)2+vz∂2vz∂y2→vz∂2vz∂y2=∂∂y(vz∂vz∂y)−(∂vz∂y)2(∂vz∂z)2+vz∂2vz∂z2→vz∂2vz∂z2=∂∂z(vz∂vz∂z)−(∂vz∂z)2 $\begin{array}{l} \frac{\partial}{\partial x} (v_{z} \frac{\partial v_{z}}{\partial x}) & = & {(\frac{\partial v_{z}}{\partial x})}^{2} + v_{z} \frac{\partial^{2} v_{z}}{\partial x^{2}} \to v_{z} \frac{\partial^{2} v_{z}}{\partial x^{2}} = \frac{\partial}{\partial x} (v_{z} \frac{\partial v_{z}}{\partial x}) - {(\frac{\partial v_{z}}{\partial x})}^{2} \\ \frac{\partial}{\partial y} (v_{z} \frac{\partial v_{z}}{\partial y}) & = & {(\frac{\partial v_{z}}{\partial y})}^{2} + v_{z} \frac{\partial^{2} v_{z}}{\partial y^{2}} \to v_{z} \frac{\partial^{2} v_{z}}{\partial y^{2}} = \frac{\partial}{\partial y} (v_{z} \frac{\partial v_{z}}{\partial y}) - {(\frac{\partial v_{z}}{\partial y})}^{2} \\ \frac{\partial}{\partial z} (v_{z} \frac{\partial v_{z}}{\partial z}) & = & {(\frac{\partial v_{z}}{\partial z})}^{2} + v_{z} \frac{\partial^{2} v_{z}}{\partial z^{2}} \to v_{z} \frac{\partial^{2} v_{z}}{\partial z^{2}} = \frac{\partial}{\partial z} (v_{z} \frac{\partial v_{z}}{\partial z}) - {(\frac{\partial v_{z}}{\partial z})}^{2} \end{array}$

si35_e

which we reinsert into Eq. 17.16 to find

==η(∂∂x(vx∂vx∂x)−(∂vx∂x)2+∂∂y(vx∂vx∂y)−(∂vx∂y)2+∂∂z(vx∂vx∂z)−(∂vx∂z)2+∂∂x(vy∂vy∂x)−(∂vy∂x)2+∂∂y(vy∂vy∂y)−(∂vy∂y)2+∂∂z(vy∂vy∂z)−(∂vy∂z)2+∂∂x(vz∂vz∂x)−(∂vz∂x)2+∂∂x(vz∂vz∂y)−(∂vz∂y)2+∂∂z(vz∂vz∂z)−(∂vz∂z)2)η(∂∂x(vx∂vx∂x+vy∂vy∂x+vz∂vz∂x)+∂∂y(vx∂vx∂y+vy∂vy∂y+vz∂vz∂y)+∂∂z(vx∂vx∂x+vy∂vy∂z+vz∂vz∂z)−((∂vx∂x)2+(∂vx∂z)2+(∂vx∂z)2+(∂vy∂x)2+(∂vy∂y)2+(∂vy∂z)2+(∂vz∂x)2+(∂vz∂y)2+(∂vz∂z)2))η⎛⎝⎜⎜⎜⎜∇→⋅⎛⎝⎜⎜⎜⎜vx∂vx∂x+vy∂vy∂x+vz∂vz∂xvx∂vx∂y+vy∂vy∂y+vz∂vz∂yvx∂vx∂z+vy∂vy∂z+vz∂vz∂z⎞⎠⎟⎟⎟⎟−((∂vx∂x)2+(∂vx∂y)2+(∂vx∂z)2+(∂vy∂x)2+(∂vy∂y)2+(∂vy∂z)2+(∂vz∂x)2+(∂vz∂y)2+(∂vz∂z)2)⎞⎠⎟⎟⎟⎟ $\begin{array}{l} η (\frac{\partial}{\partial x} (v_{x} \frac{\partial v_{x}}{\partial x}) - {(\frac{\partial v_{x}}{\partial x})}^{2} + \frac{\partial}{\partial y} (v_{x} \frac{\partial v_{x}}{\partial y}) - {(\frac{\partial v_{x}}{\partial y})}^{2} + \frac{\partial}{\partial z} (v_{x} \frac{\partial v_{x}}{\partial z}) - {(\frac{\partial v_{x}}{\partial z})}^{2} + \frac{\partial}{\partial x} (v_{y} \frac{\partial v_{y}}{\partial x}) - {(\frac{\partial v_{y}}{\partial x})}^{2} + \frac{\partial}{\partial y} (v_{y} \frac{\partial v_{y}}{\partial y}) - {(\frac{\partial v_{y}}{\partial y})}^{2} + \frac{\partial}{\partial z} (v_{y} \frac{\partial v_{y}}{\partial z}) - {(\frac{\partial v_{y}}{\partial z})}^{2} + \frac{\partial}{\partial x} (v_{z} \frac{\partial v_{z}}{\partial x}) - {(\frac{\partial v_{z}}{\partial x})}^{2} + \frac{\partial}{\partial x} (v_{z} \frac{\partial v_{z}}{\partial y}) - {(\frac{\partial v_{z}}{\partial y})}^{2} + \frac{\partial}{\partial z} (v_{z} \frac{\partial v_{z}}{\partial z}) - {(\frac{\partial v_{z}}{\partial z})}^{2}) \\ = & η (\frac{\partial}{\partial x} (v_{x} \frac{\partial v_{x}}{\partial x} + v_{y} \frac{\partial v_{y}}{\partial x} + v_{z} \frac{\partial v_{z}}{\partial x}) + \frac{\partial}{\partial y} (v_{x} \frac{\partial v_{x}}{\partial y} + v_{y} \frac{\partial v_{y}}{\partial y} + v_{z} \frac{\partial v_{z}}{\partial y}) + \frac{\partial}{\partial z} (v_{x} \frac{\partial v_{x}}{\partial x} + v_{y} \frac{\partial v_{y}}{\partial z} + v_{z} \frac{\partial v_{z}}{\partial z}) - ({(\frac{\partial v_{x}}{\partial x})}^{2} + {(\frac{\partial v_{x}}{\partial z})}^{2} + {(\frac{\partial v_{x}}{\partial z})}^{2} + {(\frac{\partial v_{y}}{\partial x})}^{2} + {(\frac{\partial v_{y}}{\partial y})}^{2} + {(\frac{\partial v_{y}}{\partial z})}^{2} + {(\frac{\partial v_{z}}{\partial x})}^{2} + {(\frac{\partial v_{z}}{\partial y})}^{2} + {(\frac{\partial v_{z}}{\partial z})}^{2})) \\ = & η (\vec{\nabla} \cdot (\begin{array}{l} v_{x} \frac{\partial v_{x}}{\partial x} + v_{y} \frac{\partial v_{y}}{\partial x} + v_{z} \frac{\partial v_{z}}{\partial x} \\ v_{x} \frac{\partial v_{x}}{\partial y} + v_{y} \frac{\partial v_{y}}{\partial y} + v_{z} \frac{\partial v_{z}}{\partial y} \\ v_{x} \frac{\partial v_{x}}{\partial z} + v_{y} \frac{\partial v_{y}}{\partial z} + v_{z} \frac{\partial v_{z}}{\partial z} \end{array}) - ({(\frac{\partial v_{x}}{\partial x})}^{2} + {(\frac{\partial v_{x}}{\partial y})}^{2} + {(\frac{\partial v_{x}}{\partial z})}^{2} + {(\frac{\partial v_{y}}{\partial x})}^{2} + {(\frac{\partial v_{y}}{\partial y})}^{2} + {(\frac{\partial v_{y}}{\partial z})}^{2} + {(\frac{\partial v_{z}}{\partial x})}^{2} + {(\frac{\partial v_{z}}{\partial y})}^{2} + {(\frac{\partial v_{z}}{\partial z})}^{2})) \end{array}$

si36_e (Eq. 17.17)

where we succeeded in introducing one gradient expression. In addition to that we also obtained a total of nine terms containing the gradients of the velocity terms.

Motivation for Applying Gauss’s Theorem. In the following we will apply Gauss’s theorem (see Eq. 7.10) to the individual terms we found. Applying this theorem will allow us to convert the volume integrals of Eq. 17.3 into surface integrals. As an example, consider the rectangular cross-section microfluidic channel shown in Fig. 17.1. To the fluid of this channel, we apply a pressure drop along the x-axis, which leads to the formation of a stable velocity profile for t < t₀. At t₀ we remove the pressure drop. We observe a fluid control volume Ω as it travels along the x-axis. Over time the flow velocity will become smaller and smaller. The control volume loses kinetic energy due to viscous dissipation. The velocity decreases as a function of time, not as a function of x. In a more general case, this would also be true for the velocity along the y-direction and the z-direction. As we will see, applying Gauss’s theorem will significantly reduce the complexity of our equations. We will start with the pressure term, Eq. 17.8.

f17-01-9781455731411 — Fig. 17.1 Viscous dissipation in a microfluidic channel with rectangular cross-section. a) A pressure gradient along the x-axis is removed at t = t₀ and a fluid plug containing the control volume Ω with the surface ∂Ω travels along the x-axis at a given speed. The velocity profile v_x (y, z) is fully-developed at this stage due to the action of the pressure drop. b) After removing the pressure drop, the fluid plug slows down over time. c/d) After a given time interval the fluid plug will come to a standstill. The reduction of the velocity profile is due to viscous dissipation.

Applying Gauss’s Theorem to the Pressure Term. Applying Gauss’s theorem to Eq. 17.8 yields

∫Ωv→⋅(∇→p)dV==∫Ω∇→⋅(v→p)dV−∮∂Ωv→p⋅n→dA $\begin{array}{l} \int_{Ω} \vec{v} \cdot (\vec{\nabla} p) d V & = & \int_{Ω} \vec{\nabla} \cdot (\vec{v} p) d V \\ = & - \oint_{\partial Ω} \vec{v} p \cdot \vec{n} d A \end{array}$

si37_e (Eq. 17.18)

It is important to note that for the boundaries along the axis the surface normals n→ $\vec{n}$ are in opposing directions. This means that the normal for the positive boundary ∂Ω_+x is pointing in opposite direction of the negative boundary ∂Ω_−x. This is true for all other directions as well. In addition we choose our coordinate system to be traveling with the control volume, thus the velocity is not a function of the position but only of time. Therefore, the contribution to the surface integrals on the (for each axis respectively) positive and (for each axis respectively) negative boundaries will cancel as long as there is no change in pressure along the respective axis. For the most general case, the contribution amounts to

∫Ωv→⋅(∇→p)dV==−(∫(vxp−vx(p+∂p∂x))dydz+∫(vyp−vy(p+∂p∂y))dxdz+∫(vzp−vz(p+∂p∂z))dxdy)∫vx∂p∂xdydz+∫vy∂p∂ydxdz+∫vz∂p∂zdxdy $\begin{array}{l} \int_{Ω} \vec{v} \cdot (\vec{\nabla} p) d V & = & - (\int (v_{x} p - v_{x} (p + \frac{\partial p}{\partial x})) d y d z + \int (v_{y} p - v_{y} (p + \frac{\partial p}{\partial y})) d x d z + \int (v_{z} p - v_{z} (p + \frac{\partial p}{\partial z})) d x d y) \\ = & \int v_{x} \frac{\partial p}{\partial x} d y d z + \int v_{y} \frac{\partial p}{\partial y} d x d z + \int v_{z} \frac{\partial p}{\partial z} d x d y \end{array}$

si38_e (Eq. 17.19)

Applying Gauss’s Theorem to the Convection Term. The next contribution we will discuss is due to the convection term, Eq. 17.15. Applying Gauss’s theorem yields

∫Ωv→(v→⋅∇→)v→dV===∫Ω12∇→⋅v→(v2x+v2y+v2z)dV12∮∂Ωv→(v2x+v2y+v2z)⋅n→dA0 $\begin{array}{l} \int_{Ω} \vec{v} (\vec{v} \cdot \vec{\nabla}) \vec{v} d V & = & \int_{Ω} \frac{1}{2} \vec{\nabla} \cdot \vec{v} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) d V \\ = & \frac{1}{2} \oint_{\partial Ω} \vec{v} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) \cdot \vec{n} d A \\ = & 0 \end{array}$

si39_e

As before, the negative boundary and positive boundary contributions will cancel due to the fact that velocity does not change along the respective axis. As the surface normals on the respective boundaries point in the opposite direction with the same magnitude, the overall contribution to the dissipation is zero. Here it becomes very obvious why applying Gauss’s theorem is advantageous.

Applying Gauss’s Theorem to the Viscosity Term. The last term we will look at is the viscosity term Eq. 17.17. Applying Gauss’s theorem yields

∫Ωv→ηΔv→dV===η∫Ω∇→⋅⎛⎝⎜⎜⎜⎜vx∂vx∂x+vy∂vy∂x+vz∂vz∂xvx∂vx∂y+vy∂vy∂y+vz∂vz∂yvx∂vx∂z+vy∂vy∂z+vz∂vz∂z⎞⎠⎟⎟⎟⎟dV−η∫Ω((∂vx∂x)2+(∂vx∂y)2+(∂vx∂z)2+(∂vy∂x)2+(∂vy∂y)2+(∂vy∂z)2+(∂vz∂x)2+(∂vz∂y)2+(∂vz∂z)2)dVη∮⎛⎝⎜⎜⎜⎜vx∂vx∂x+vy∂vy∂x+vz∂vz∂xvx∂vx∂y+vy∂vy∂y+vz∂vz∂yvx∂vx∂z+vy∂vy∂z+vz∂vz∂z⎞⎠⎟⎟⎟⎟⋅n→dA−η∫Ω((∂vx∂x)2+(∂vx∂y)2+(∂vx∂z)2+(∂vy∂x)2+(∂vy∂y)2+(∂vy∂z)2+(∂vz∂x)2+(∂vz∂y)2+(∂vz∂z)2)dV−η∫Ω((∂vx∂x)2+(∂vx∂y)2+(∂vx∂z)2+(∂vy∂x)2+(∂vy∂y)2+(∂vy∂z)2+(∂vz∂x)2+(∂vz∂y)2+(∂vz∂z)2)dV $\begin{array}{l} \int_{Ω} \vec{v} η Δ \vec{v} d V & = & η \int_{Ω} \vec{\nabla} \cdot (\begin{matrix} v_{x} \frac{\partial v_{x}}{\partial x} + v_{y} \frac{\partial v_{y}}{\partial x} + v_{z} \frac{\partial v_{z}}{\partial x} \\ v_{x} \frac{\partial v_{x}}{\partial y} + v_{y} \frac{\partial v_{y}}{\partial y} + v_{z} \frac{\partial v_{z}}{\partial y} \\ v_{x} \frac{\partial v_{x}}{\partial z} + v_{y} \frac{\partial v_{y}}{\partial z} + v_{z} \frac{\partial v_{z}}{\partial z} \end{matrix}) d V - η \int_{Ω} ({(\frac{\partial v_{x}}{\partial x})}^{2} + {(\frac{\partial v_{x}}{\partial y})}^{2} + {(\frac{\partial v_{x}}{\partial z})}^{2} + {(\frac{\partial v_{y}}{\partial x})}^{2} + {(\frac{\partial v_{y}}{\partial y})}^{2} + {(\frac{\partial v_{y}}{\partial z})}^{2} + {(\frac{\partial v_{z}}{\partial x})}^{2} + {(\frac{\partial v_{z}}{\partial y})}^{2} + {(\frac{\partial v_{z}}{\partial z})}^{2}) d V \\ = & η \oint (\begin{matrix} v_{x} \frac{\partial v_{x}}{\partial x} + v_{y} \frac{\partial v_{y}}{\partial x} + v_{z} \frac{\partial v_{z}}{\partial x} \\ v_{x} \frac{\partial v_{x}}{\partial y} + v_{y} \frac{\partial v_{y}}{\partial y} + v_{z} \frac{\partial v_{z}}{\partial y} \\ v_{x} \frac{\partial v_{x}}{\partial z} + v_{y} \frac{\partial v_{y}}{\partial z} + v_{z} \frac{\partial v_{z}}{\partial z} \end{matrix}) \cdot \vec{n} d A - η \int_{Ω} ({(\frac{\partial v_{x}}{\partial x})}^{2} + {(\frac{\partial v_{x}}{\partial y})}^{2} + {(\frac{\partial v_{x}}{\partial z})}^{2} + {(\frac{\partial v_{y}}{\partial x})}^{2} + {(\frac{\partial v_{y}}{\partial y})}^{2} + {(\frac{\partial v_{y}}{\partial z})}^{2} + {(\frac{\partial v_{z}}{\partial x})}^{2} + {(\frac{\partial v_{z}}{\partial y})}^{2} + {(\frac{\partial v_{z}}{\partial z})}^{2}) d V \\ = & - η \int_{Ω} ({(\frac{\partial v_{x}}{\partial x})}^{2} + {(\frac{\partial v_{x}}{\partial y})}^{2} + {(\frac{\partial v_{x}}{\partial z})}^{2} + {(\frac{\partial v_{y}}{\partial x})}^{2} + {(\frac{\partial v_{y}}{\partial y})}^{2} + {(\frac{\partial v_{y}}{\partial z})}^{2} + {(\frac{\partial v_{z}}{\partial x})}^{2} + {(\frac{\partial v_{z}}{\partial y})}^{2} + {(\frac{\partial v_{z}}{\partial z})}^{2}) d V \end{array}$

si40_e (Eq. 17.20)

Here again, we have dropped the surface integral due to the fact that the velocities v_x (x) = v_x (x + dx), v_y (y) = v_y (y + dy), and v_z (z) = v_z (z + dz), and therefore the individual contributions to the surface integral cancel along each axis, respectively. The only remaining contribution originates from the gradients of the velocity profiles.

Equation for the Viscous Dissipation. Having found Eq. 17.18, Eq. 17.19, and Eq. 17.20, which we can reinsert into Eq. 17.3, we can now formulate the equation for the viscous dissipation as

∂Ekin∂t=∫Ω(vxk→x+vyk→y+vzk→z)dV+∫vx∂p∂xdydz+∫vy∂p∂ydxdz+vz∂p∂zdxdy−η∫Ω((∂vx∂x)2+(∂vx∂y)2+(∂vx∂z)2+(∂vy∂x)2+(∂vy∂y)2+(∂vy∂z)2+(∂vz∂x)2+(∂vz∂y)2+(∂vz∂z)2)dV $\frac{\partial E_{kin}}{\partial t} = \int_{Ω} (v_{x} {\vec{k}}_{x} + v_{y} {\vec{k}}_{y} + v_{z} {\vec{k}}_{z}) d V + \int v_{x} \frac{\partial p}{\partial x} d y d z + \int v_{y} \frac{\partial p}{\partial y} d x d z + v_{z} \frac{\partial p}{\partial z} d x d y - η \int_{Ω} ({(\frac{\partial v_{x}}{\partial x})}^{2} + {(\frac{\partial v_{x}}{\partial y})}^{2} + {(\frac{\partial v_{x}}{\partial z})}^{2} + {(\frac{\partial v_{y}}{\partial x})}^{2} + {(\frac{\partial v_{y}}{\partial y})}^{2} + {(\frac{\partial v_{y}}{\partial z})}^{2} + {(\frac{\partial v_{z}}{\partial x})}^{2} + {(\frac{\partial v_{z}}{\partial y})}^{2} + {(\frac{\partial v_{z}}{\partial z})}^{2}) d V$

si41_e (Eq. 17.21)

Eq. 17.21 contains all the relevant terms. We will now simplify this equation to the case depicted in Fig. 17.1. In this case, we have no volume forces and (for t > t₀) also no pressure drop. The velocities v_y = v_z = 0, and only v_x is non-zero. As the velocity is only dependent on time but not on position, there is no gradient in the velocity along the x-axis, and therefore ∂vx∂x=0 $\frac{\partial v_{x}}{\partial x} = 0$ . We also assume that there is no actual pressure drop along any of the axes driving the flow. As we know, pressure drops can be used to balance viscous dissipation and therefore drive the flow. Here we want to study the effect of dissipation in a flow that is not driven. We therefore set ∂p∂x=∂p∂y=∂p∂z=0 $\frac{\partial p}{\partial x} = \frac{\partial p}{\partial y} = \frac{\partial p}{\partial z} = 0$ . Therefore, we can simplify Eq. 17.21 to

∂Ekin∂t=−η∫Ω((∂vx∂y)2+(∂vx∂z)2)dV $\frac{\partial E_{kin}}{\partial t} = - η \int_{Ω} ({(\frac{\partial v_{x}}{\partial y})}^{2} + {(\frac{\partial v_{x}}{\partial z})}^{2}) d V$

si44_e

As stated in the introduction to this section, the decline in kinetic energy, i.e., the viscous dissipation, is dependent only on the gradient of the velocity profile and not (as one may expect) on the absolute velocities. Returning to our original example of the squirt gun this can be verified easily. If a higher velocity on the exit of the nozzle of the gun is desired, the central velocity in the (presumably circular) channel must be higher. Assuming no-slip boundary conditions the gradients of the velocity profile will be higher if the center velocity is higher. Please note that the increase in steepness of the velocity profile contributes as a square. Therefore small increases in the gradient will quickly lead to significant pressure drops in the nozzle.

Balancing Viscous Dissipation. Obviously, in a technical system, the flow velocity should usually be maintained. In this case, a driving external or volume force must be applied in order to compensate for the viscous dissipation. In case of the rectangular microchannel shown in Fig. 17.1 we can maintain the pressure drop along the x-axis in order to compensate the loss of kinetic energy. In this case we find

∂Ekin∂tη∫Ω((∂vx∂y)2+(∂vx∂z)2)dV====∫vx∂p∂xdydz−η∫Ω((∂vx∂y)2+(∂vx∂z)2)dV!−−−−0∫vx∂p∂xdydz∫vxΔpdydzΔp∫vxdydz $\begin{array}{l} \frac{\partial E_{kin}}{\partial t} & = & \int v_{x} \frac{\partial p}{\partial x} d y d z - η \int_{Ω} ({(\frac{\partial v_{x}}{\partial y})}^{2} + {(\frac{\partial v_{x}}{\partial z})}^{2}) d V \underline{\underline{!}} 0 \\ η \int_{Ω} ({(\frac{\partial v_{x}}{\partial y})}^{2} + {(\frac{\partial v_{x}}{\partial z})}^{2}) d V & = & \int v_{x} \frac{\partial p}{\partial x} d y d z \\ = & \int v_{x} Δ p d y d z \\ = & Δ p \int v_{x} d y d z \end{array}$

si45_e (Eq. 17.22)

=ΔpQ $= Δ p Q$

(Eq. 17.23)

where we have used Eq. 10.9 replacing the integral with the flowrate. Eq. 17.23 is a very important equation. It sets the pressure drop and the flowrate into relation.

Physical Unit. Looking at the left-hand side of Eq. 17.23 we can see that the unit of viscous dissipation is

Js−1=Nms−1=Nm3m−2s−1=Pam3s−1 ${Js}^{- 1} = N m s^{- 1} = N m^{3} m^{- 2} s^{- 1} = Pa m^{3} s^{- 1}$

which corresponds exactly to the units of the right-hand side of Eq. 17.23.

17.3 Hydraulic Resistance of Important flow Channel Geometries

17.3.1 Overview

We have now derived the correlation between the loss of kinetic energy due to viscous dissipation and the product of the flowrate and the pressure drop (Eq. 17.23). We already know from the various flow channel geometries we have studied so far that a constant pressure drop results in a constant flowrate:

Q ∝ ∆p

The most commonly encountered form of this correlation is the Hagen-Poiseuille law (see section 16.3.2) for which we found Eq. 16.22:

QΔp==−πR48η−8ηπR4QΔl $\begin{array}{l} Q & = & - \frac{π R^{4}}{8 η} \\ Δ p & = & - \frac{8 η}{π R^{4}} Q Δ l \end{array}$

si48_e

ΔpΔl=−Rn,hydQ $\begin{array}{l} \frac{Δ p}{Δ l} & = & - R_{n,hyd} Q \end{array}$

(Eq. 17.24)

Rn,hyd=−ΔpΔl1Q $\begin{array}{l} R_{n,hyd} & = & - \frac{Δ p}{Δ l} \frac{1}{Q} \end{array}$

(Eq. 17.25)

where we have introduced the normalized hydraulic resistance R_n,hyd as the proportionality factor. The unit of this factor is kg m⁻⁵ s⁻¹. It gives the pressure loss per unit of channel length l and per unit of the flowrate. The total hydraulic resistance can be calculated as

Rn,hydRhyd==RhydΔlRn,hydΔl $\begin{array}{l} R_{n,hyd} & = & \frac{R_{hyd}}{Δ l} \\ R_{hyd} & = & R_{n,hyd} Δ l \end{array}$

si51_e

by simply multiplying the normalized hydraulic resistance with the channel length. The unit of the hydraulic resistance is kg m⁻⁴ s⁻¹. It depends solely on the properties of the fluid and the geometry and can be derived for various geometries. Tab. 17.1 lists the most important channel cross-sections and the expressions for the respective hydraulic resistance.

Tab. 17.1

Hydraulic resistance values for some of the most important channel geometries

Cross-section	Normalized (per unit channel length) hydraulic resistance Rn,hyd=Rhydl $R_{n,hyd} = \frac{R_{hyd}}{l}$	Derived in
Elliptical	4ηπa2+b2a3b3 $\frac{4 η}{π} \frac{a^{2} + b^{2}}{a^{3} b^{3}}$	Eq. 16.5
Circular	8ηπr4 $\frac{8 η}{π r^{4}}$	Eq. 16.8
Planar infinitesimally extended	12ηwh3 $\frac{12 η}{w h^{3}}$	Eq. 16.14
Rectangular	12ηh3w⋅11−192hπ5w∑∞n=01(2n+1)5tanh((2n+1)πw2h) $\frac{12 η}{h^{3} w} \cdot \frac{1}{1 - \frac{192 h}{π^{5} w} \sum_{n = 0}^{\infty} \frac{1}{{(2 n + 1)}^{5}} \tanh ((2 n + 1) π \frac{w}{2 h})}$	Eq. 16.64
Square	12ηw4⋅11−192π5∑∞n=01(2n+1)5tanh((2n+1)π2) $\frac{12 η}{w^{4}} \cdot \frac{1}{1 - \frac{192}{π^{5}} \sum_{n = 0}^{\infty} \frac{1}{{(2 n + 1)}^{5}} \tanh ((2 n + 1) \frac{π}{2})}$	Eq. 16.65

Please note that the results of Tab. 17.1 are only valid in cases where the non-linear convection term (u→⋅∇→)u→ $(\vec{u} \cdot \vec{\nabla}) \vec{u}$ of the Navier-Stokes equation (see Eq. 11.40) can be neglected. This is due to the fact that all derivation of the flow profiles we made were based on the assumption of having Poiseuille flow, i.e., laminar flow conditions with the only velocity component being v_x and v_y = v_z = 0 (see Eq. 15.14). See section 15.4.1 for details on these assumptions.

In order to facilitate working with the equations of Tab. 17.1, Fig. 17.2, Fig. 17.3, and Fig. 17.4 show the pressure drops calculated for microchannels with the five different geometries when probed with water. Fig. 17.2 shows the pressure drop for elliptical and circular cross-sections, Fig. 17.3 shows the pressure drop for the parallel infinitesimally extended channel, and Fig. 17.4 shows the pressure drop for the rectangular and square channel cross-sections. The pressure drops are given as pressure drop per mm of channel length ΔpΔl $\frac{Δ p}{Δ l}$ per µl/ min of flowrate. In order to calculate the total pressure drop, the value read from the charts must therefore be multiplied by the flowrate (see section 17.3.2).

f17-02-9781455731411 — Fig. 17.2 Pressure drop in channels with elliptical (a) and circular (b) cross-sections calculated for water. The pressure drop is given as the pressure drop over channel length l per flowrate Q. In order to obtain the pressure drop per mm of channel length the value read must be multiplied with the flowrate (given in µl min⁻¹).

f17-03-9781455731411 — Fig. 17.3 Normalized hydraulic resistance in infinitesimally extended parallel channels with height h and width w calculated for water. The pressure drop is given as the pressure drop over channel length l per flowrate Q. In order to obtain the pressure drop per mm of channel length the value read must be multiplied with the flowrate (given in µl min⁻¹).

f17-04-9781455731411 — Fig. 17.4 Normalized hydraulic resistance in channels with rectangular (a) and square (b) cross-sections calculated for water. The pressure drop is given as the pressure drop over channel length l per flowrate Q. In order to obtain the pressure drop per mm of channel length the value read must be multiplied with the flowrate (given in µl min⁻¹).

17.3.2 Pressure Drop: Examples

As an example, for a microfluidic channel width a height of 50 µm and width of 100 µm we can derive a pressure drop per channel length and flowrate of 0.233 mbar min /mm µl from Tab. 17.1. This value can be roughly read also from Fig. 17.4a. If we want to use such a channel in an experiment where we intend to probe a flow of 50 µl/ min we expect a total pressure drop of 11.66 mbar mm⁻¹. This means that per mm of channel length, we would experience a pressure drop of 11.66 mbar.

As a second example, we consider a circular capillary with an inner diameter of 40 µm. For this geometry we derive a pressure drop per length and flowrate of 0.16 mbar min /mm µl from Tab. 17.1 and Fig. 17.2b. As you can see, this value is lower than the values we found for the rectangular channel cross-section that we choose to have approximately the same cross-section area. If we use this capillary in an experiment where we want to probe a total flowrate of 1 mL/ min through the system, we would expect a total pressure drop of 165.78 mbar mm⁻¹. Assuming we want to use a capillary of a total length of 1 m we expect a considerable pressure drop of about 165.78 bar. This value is pretty realistic for pressure drops expected, e.g., in high-performance liquid chromatography (HPLC) systems.

17.4 Simplification Approaches to Hydraulic Resistances

17.4.1 Normalized Hydraulic Resistance

The first simplification is one we already made implicitly. As the resistance depends linearly on the length of the channel it seems reasonable to define the normalized hydraulic resistance as

Rn,hyd=Rhydlwith the unit kg m−5s−1. $\begin{matrix} R_{n,hyd} = \frac{R_{hyd}}{l} \\ {with the unit kg m}^{- 5} s^{- 1} . \end{matrix}$

si54_e (Eq. 17.26)

17.4.2 Geometric Hydraulic Resistance

As we have seen from Tab. 17.1 the hydraulic resistance is dependent on the properties of the fluid and the geometry of the cross-section of the channel. As we can see from the equations derived, the only property of the fluid which must be taken into account is the viscosity η. No other fluid property is of importance. Therefore, it makes sense to define the hydraulic as a fluid-independent property that is referred to as geometric hydraulic resistance. It is defined as

RA,geom=Rn,hydη $R_{A,geom} = \frac{R_{n,hyd}}{η}$

which is again normalized to the channel length. As the name suggests, the geometric hydraulic resistance is only dependent on the geometry of the channel and not on the properties of the fluid. Fig. 17.2, Fig. 17.3, and Fig. 17.4 have been calculated for water. By dividing the values of the charts by the viscosity of water in mPa s (which is a division by 1 as the viscosity of water is 1 mPa s) we would obtain charts displaying the geometric hydraulic resistance. As the data would remain unchanged, we would simply need to adapt the vertical axis by dividing the units displayed by mPa s.

17.4.3 Cross-Section Hydraulic Resistance

Until now, we have derived hydraulic resistances that are independent of channel length and fluid. It only remains to find a version of the hydraulic resistance that is independent of the geometry. We have learned that the integral ∫_Ωv_xdA is used to define the viscous dissipation (see Eq. 17.22). Obviously, for the hydraulic resistance this integral ends up being moved into the denominator (see Tab. 17.1). Therefore it seems reasonable that as the differential of the cross-section of the channel dA is important, the total cross-section area A should play a major role. We therefore define the cross-section hydraulic resistance as

RA,hyd=ηlA2 $R_{A,hyd} = \frac{η l}{A^{2}}$

Please note that there is no derivation to this equation; it merely follows from intuition and heuristic analysis as well as experimental observation.

17.4.4 Correction Factor α

As we can see R_A,hyd is again linearly dependent on the viscosity and the length and could be easily normalized. R_A,hyd is obviously significantly easier to derive than the integral ∫_Ωv_xdA. We can now put the cross-section hydraulic resistance in relation to the actual hydraulic resistance as

RhydRA,hyd==A2ηlRhydA2RA,geom=α $\begin{array}{l} \frac{R_{hyd}}{R_{A,hyd}} & = & \frac{A^{2}}{η l} R_{hyd} \\ = & A^{2} R_{A,geom} = α \end{array}$

si57_e (Eq. 17.27)

where we have introduced the correction factor α. It defines how well the square of the area matches the integral ∫_Ωv_xdA used to calculate the hydraulic resistance. For α = 1 the integral is matched perfectly. In the following we will quickly calculate the correction factors for the geometries of Tab. 17.1.

Correction Factor for Elliptical Cross-Sections. For elliptical cross-sections we find the correction factor to be

αelliptical==(πab)24πa2+b2a3b34πa2+b2ab=4π(ab+ba) $\begin{array}{l} α_{elliptical} & = & {(π a b)}^{2} \frac{4}{π} \frac{a^{2} + b^{2}}{a^{3} b^{3}} \\ = & 4 π \frac{a^{2} + b^{2}}{a b} = 4 π (\frac{a}{b} + \frac{b}{a}) \end{array}$

si58_e

which we can simplify using the aspect ratio (see section 9.9.15) r=ba $r = \frac{b}{a}$ to

αelliptical=4π(1r+r) $α_{elliptical} = 4 π (\frac{1}{r} + r)$

(Eq. 17.28)

As can be seen, the correction factor for elliptical cross-sections is a function of the aspect ratio only.

Correction Factor for Circular Cross-Sections. For circular cross-sections we find the correction factor to be

αcircular==4π(rr+rr)8π $\begin{array}{l} α_{circular} & = & 4 π (\frac{r}{r} + \frac{r}{r}) \\ = & 8 π \end{array}$

si61_e (Eq. 17.29)

As can be seen, the correction factor for circular cross-sections is a constant.

Correction Factor for Planar Infinitesimally Extended Cross-Sections. For planar infinitesimally extended channel cross-sections we find the correction factor to be

αinfint.extended===(wh)212wh312wh12r $\begin{array}{l} α_{infint . e x t e n d e d} & = & {(w h)}^{2} \frac{12}{w h^{3}} \\ = & 12 \frac{w}{h} \\ = & \frac{12}{r} \end{array}$

si62_e (Eq. 17.30)

for which we have again used the aspect ratio. As can be seen, the correction factor scales inversely proportional to the aspect ratio.

Correction Factor for Rectangular Cross-Sections. For rectangular cross-sections we find the correction factor to be

αrectangular===(wh)212h3w11−192hπ5w∑∞n=01(2n+1)5tanh((2n+1)πw2h)12wh11−192hπ5w∑∞n=01(2n+1)5tanh((2n+1)πw2h)12r11−192rπ5∑∞n=01(2n+1)5tanh((2n+1)π12r) $\begin{array}{l} α_{rectangular} & = & {(w h)}^{2} \frac{12}{h^{3} w} \frac{1}{1 - \frac{192 h}{π^{5} w} \sum_{n = 0}^{\infty} \frac{1}{{(2 n + 1)}^{5}} \tanh ((2 n + 1) π \frac{w}{2 h})} \\ = & 12 \frac{w}{h} \frac{1}{1 - \frac{192 h}{π^{5} w} \sum_{n = 0}^{\infty} \frac{1}{{(2 n + 1)}^{5}} \tanh ((2 n + 1) π \frac{w}{2 h})} \\ = & \frac{12}{r} \frac{1}{1 - \frac{192 r}{π^{5}} \sum_{n = 0}^{\infty} \frac{1}{{(2 n + 1)}^{5}} \tanh ((2 n + 1) π \frac{1}{2 r})} \end{array}$

si63_e (Eq. 17.31)

where we have again used the aspect ratio.

Correction Factor for Square Cross-Sections. For square cross-sections we find the correction factor to be

αsquare==w412w411−192π5∑∞n=01(2n+1)5tanh((2n+1)π2)121−192π5∑∞n=01(2n+1)5tanh((2n+1)π2) $\begin{array}{l} α_{square} & = & w^{4} \frac{12}{w^{4}} \frac{1}{1 - \frac{192}{π^{5}} \sum_{n = 0}^{\infty} \frac{1}{{(2 n + 1)}^{5}} \tanh ((2 n + 1) \frac{π}{2})} \\ = & \frac{12}{1 - \frac{192}{π^{5}} \sum_{n = 0}^{\infty} \frac{1}{{(2 n + 1)}^{5}} \tanh ((2 n + 1) \frac{π}{2})} \end{array}$

si64_e (Eq. 17.32)

17.4.5 Compactness Factor C

Having introduced the correction factor we now introduce a second factor commonly used to characterize mi-crofluidic channel cross-sections, i.e., the compactness factor C. This factor is purely geometrically defined and puts the two relevant areas of a microfluidic channel into correlation: the cross-section and the circumference. It is defined as

C=(circumference)2cross-sectionarea=(∫∂Ωdl)2∫ΩdA $C = \frac{{(circumference)}^{2}}{cross-section area} = \frac{{(\int_{\partial Ω} d l)}^{2}}{\int_{Ω} d A}$

si65_e (Eq. 17.33)

As you can see, the compactness factor uses the square of the circumference. This is due to the fact that many effects in microfluidics are relying on surface effects. Due to the high surface-to-volume ratios in microfluidics these surface effects play a major role. As we have learned viscous dissipation is due to fluid friction which is highly dependent on the boundary. As in most microfluidic systems we have no-slip boundary conditions it is reasonable to assume that these effects play a major role in most fluid mechanical considerations.

Rationale for the Compactness Factor. Obviously, the compactness C is reasonably easy to calculate as it simply requires analytical equations for cross-section and circumference of a geometry. On the other hand, we have found a reasonably convenient way to calculate the geometric hydraulic resistance R_A,geom using the correction factor α. Now let us assume we would find a simple correlation between α and C, i.e., a function α (C). In this case we would be able to calculate R_A,geom pretty much without any effort by simply dividing the correction factor by the cross-section area using Eq. 17.27.

As we will see in a moment, for many cross-section geometries such a function does exist. We will demonstrate this using the five geometries we have studied so far.

Correlation Between Compactness and Correction Factor for Elliptical Cross-Sections. Interestingly, the first cross-section we will study is the most complicated one. This is due to the fact that the circumference of an ellipse cannot be evaluated analytically. The integral is

lc,ellipse=4∫π40(acosϕ)2+(bsinϕ)2−−−−−−−−−−−−−−−−−−√dϕ $l_{c,ellipse} = 4 \int_{0}^{\frac{π}{4}} \sqrt{{(a \cos ϕ)}^{2} + {(b \sin ϕ)}^{2}} d ϕ$

si66_e

We will retain this integral notation for now. The compactness factor for the ellipse can now be calculated as

Cellipse===16πab(∫π40(acosϕ)2+(bsinϕ)2−−−−−−−−−−−−−−−−−−−√dϕ)216aπb(∫π40(cosϕ)2+(basinϕ)2−−−−−−−−−−−−−−−−−−−√dϕ)216π1r(∫π40(cosϕ)2+(rsinϕ)2−−−−−−−−−−−−−−−−−√dϕ)2 $\begin{array}{l} C_{ellipse} & = & \frac{16}{π a b} {(\int_{0}^{\frac{π}{4}} \sqrt{{(a \cos ϕ)}^{2} + {(b \sin ϕ)}^{2}} d ϕ)}^{2} \\ = & \frac{16 a}{π b} {(\int_{0}^{\frac{π}{4}} \sqrt{{(\cos ϕ)}^{2} + {(\frac{b}{a} \sin ϕ)}^{2}} d ϕ)}^{2} \\ = & \frac{16}{π} \frac{1}{r} {(\int_{0}^{\frac{π}{4}} \sqrt{{(\cos ϕ)}^{2} + {(r \sin ϕ)}^{2}} d ϕ)}^{2} \end{array}$

si67_e (Eq. 17.34)

As you can see Eq. 17.34 is an analytical equation for the compactness factor that depends solely on the aspect ratio r=ba $r = \frac{b}{a}$ as the only independent variable. As stated, the elliptical circumference cannot be evaluated analytically so we must plug in numbers for r. On the other hand, we have found the correction factor for the ellipse (see Eq. 17.28) also to be an equation solely dependent on the aspect ratio as independent variable. We will now plug in values for r and plot the results in a chart using the values of C (r) on the abscissa and the values of α (r) on the ordinate. This plot is shown in Fig. 17.5 for different values of r. As you can see, the values fall on a line. The equation can be obtained by linear fitting to be

f17-05-9781455731411 — Fig. 17.5 Correction factor α over compactness factor C for elliptical cross-sections for different values of the aspect ratio r. As you can see, the data points fall on a line (α = 0.40C + 2.16, R² = 0.99).

αellipse(Cellipse)=0.40⋅Cellipse+2.16 $α_{ellipse} (C_{ellipse}) = 0.40 \cdot C_{ellipse} + 2.16$

(Eq. 17.35)

Obviously, this is a very useful finding. For a given aspect ratio we would need to calculate C_ellipse using Eq. 17.34 then apply Eq. 17.35 to find α_ellipse. With this value we can find the hydraulic resistance using Eq. 17.27 simply by division with the squared cross-section area. See section 17.4.5 for an example of such a calculation. Tab. 17.2 shows some of the most commonly encountered aspect ratios and the corresponding compactness C and correction values α for elliptical cross-sections.

Tab. 17.2

Compactness C and correction factor α values for elliptical cross-sections with different aspect ratios r

Aspect ratio r [-]	Compactness factor C [-]	Correction factor α [-]
0.001	5 092.998	12 566.383
0.010	509.576	1 256.763
0.100	52.572	126.920
1.000	12.566	25.133
2.000	14.939	31.416
5.000	28.102	65.345
10.000	52.572	126.920
20.000	102.851	251.956
50.000	255.137	628.570
100.000	509.576	1 256.763

Correlation Between Compactness and Correction Factor for Circular Cross-Sections. Obviously, circular cross-sections are simply a special case of elliptical cross-sections. Nevertheless we will quickly derive the equations and ensure, using Tab. 17.2 (for r = 1), that we end up with the correct values. The compactness value for circular cross-sections is simply

Ccircular=≈(2πr)2πr2=4π12.566 $\begin{array}{l} C_{circular} & = & \frac{{(2 π r)}^{2}}{π r^{2}} = 4 π \\ \approx & 12.566 \end{array}$

si70_e (Eq. 17.36)

which corresponds with the value in Tab. 17.2.

Using the expression we found for the correction factor Eq. 17.29 we find that

αcircular(Ccircular)=≈2Ccircular25.132 $\begin{array}{l} α_{circular} (C_{circular}) & = & 2 C_{circular} \\ \approx & 25.132 \end{array}$

(Eq. 17.37)

which again corresponds with the value in Tab. 17.2.

Correlation Between Compactness and Correction Factor for Planar Infinitesimally Extended Cross-Sections. For planar infinitesimally extended channel cross-sections we find the compactness factor to be

Cinfinit, extended===(2(w+h))2wh=4wh1(+wh)24r(1+r)24r+8+4r $\begin{array}{l} C_{infinit, extended} & = & \frac{{(2 (w + h))}^{2}}{w h} = 4 \frac{w}{h} 1 {(+ \frac{w}{h})}^{2} \\ = & \frac{4}{r} (1 + r) 2 \\ = & 4 r + 8 + \frac{4}{r} \end{array}$

si72_e (Eq. 17.38)

We will not try to express this as an expression of the correction factor given by Eq. 17.30. Instead we will again evaluate it again using numbers for r. As you can see from Fig. 17.6 we again end up with a linear correlation between compactness C and correction factor α given by the following equation obtained by linear fitting to the data points:

f17-06-9781455731411 — Fig. 17.6 Correction factor α over compactness factor C for planar infinitesimally extended cross-sections for different values of the aspect ratio r. As you can see, the data points fall on a line (α = 3.00 · C − 24.73, R² = 1.00).

Cinfinit.extended(Cinfinit.extended)=3.00⋅Cellipse−24.73 $C_{infinit . extended} (C_{infinit . extended}) = 3.00 \cdot C_{ellipse} - 24.73$

(Eq. 17.39)

Please note that our assumption for this flow profile states that the aspect ratio should be small. Thus it makes no sense to use this formula in the range of r ≥ 1. Tab. 17.3 shows some of the most commonly encountered aspect ratios and the corresponding compactness C and correction values α for this cross-section.

Tab. 17.3

Compactness C and correction factor α values for planar infinitesimally extended cross-sections with different aspect ratios r

Aspect ratio r [-]	Compactness factor C [-]	Correction factor α [-]
0.000 1	120 000	40 008.000 4
0.001 0	12 000	4 008.004 0
0.010 0	1 200	408.040 0
0.100 0	120	48.400 0

Correlation Between Compactness and Correction Factor for Rectangular Cross-Sections. We have already derived the compactness factor for rectangular cross-sections as this equation is identical to the compactness factor found for the infinitesimally extended channel cross-section (see Eq. 17.39). The correction factor is given by Eq. 17.31. Again, we calculate both factors for different aspect ratio values r. The resulting data points are shown in Fig. 17.7. As you can see, they again all fall on a line that can be fitted to the following equation:

f17-07-9781455731411 — Fig. 17.7 Correction factor α over compactness factor C for rectangular cross-sections for different values of the aspect ratio r. As you can see, the data points fall on a line (α = 0.33 · C + 6.29, R² = 0.99).

αrectangular(Crectangular)=0.33⋅Crectangular+6.29 $α_{rectangular} (C_{rectangular}) = 0.33 \cdot C_{rectangular} + 6.29$

(Eq. 17.40)

Tab. 17.4 shows some of the most commonly encountered aspect ratios and the corresponding compactness C and correction values α for this cross-section.

Tab. 17.4

Compactness C and correction factor α values for rectangular cross-sections with different aspect ratios r

Aspect ratio r [-]	Compactness factor C [-]	Correction factor α [-]
0.001	4 008.00	12 007.57
0.010	408.04	1 207.61
0.100	48.40	128.07
0.200	28.80	68.65
0.500	18.00	34.98
1.000	16.00	28.45
2.000	18.00	34.98
5.000	28.80	68.65
10.000	48.40	128.07
20.000	88.20	247.81
50.000	208.08	607.60
100.000	408.04	1 206.77

Correlation Between Compactness and Correction Factor for Square Cross-Sections. Obviously, square cross-sections are special cases of the more general rectangular cross-sections and can be read from Tab. 17.4 for the case r = 1. In order to verify we will quickly derive the compactness factor for square cross-sections, which is

Csquare=(4w)2w2=16 $C_{s q u a r e} = \frac{{(4 w)}^{2}}{w^{2}} = 16$

si75_e (Eq. 17.41)

This value is correct as we can see from Fig. 17.7. The correction factor is given by Eq. 17.32. It can be evaluated numerically, e.g., using Maple, to a total value of 28.45, which is also correct according to Tab. 17.4.

Derivation of the Hydraulic Resistance of a Rectangular Cross-Section Microchannel Using the Compactness C and Correction α Approach. Now the following question may arise: In what way are these compactness and correction factors actually useful? We will quickly illustrate this at a simple example. Suppose we want to calculate the pressure drop in a microfluidic channel with rectangular cross-section with a width of 100 µm and a height of 50 µm, thus with an aspect ratio of r = 0.5. In section 17.3.2 we have already derived the pressure drop at this channel using Fig. 17.4a. Using Tab. 17.4 we quickly find a correction factor of 34.98 for this aspect ratio. Thus, using Eq. 17.27 we can calculate the geometric resistance value of this cross-section to be

RA,geom, rect, r=0.5==αA2=34.98(100×50μm2)21.399×1015μm4=1.399×1015kgm−5s−1mPa−1s−1 $\begin{array}{l} R_{A, geom, rect, r = 0.5} & = & \frac{α}{A^{2}} = \frac{34.98}{{(100 \times 50 {μm}^{2})}^{2}} \\ = & 1.399 \times 10^{15} {μm}^{4} = 1.399 \times 10^{15} kg m^{- 5} s^{- 1} {mPa}^{- 1} s^{- 1} \end{array}$

si76_e

We now multiply with the viscosity of water, which is 1 mPa s, to find R_n,hyd using Eq. 17.26 as

Rn,hyd,rect, r=0.5==RA,geom,rect.,r=0.5η=1.399×1015kgm−5s−10.233mbarmin/mmμl $\begin{array}{l} R_{n, hyd, rect, r = 0.5} & = & R_{A, g e o m, r e c t ., r = 0.5} η = 1.399 \times 10^{15} kg m^{- 5} s^{- 1} \\ = & 0.233 mbar \min / m m μl \end{array}$

si77_e

In section 17.3.2 we have determined a value of 0.16 mbar min /mm µl for R_n,hyd analytically, which is reasonably close.

If we want to redo the calculation for a microchannel with the same aspect ratio but with twice the width and height, i.e., a width 200 µm and a height of 100 µm, we end up with the same correction factor, α = 34.98. We would only need to update the area as required for Eq. 17.27. Within a few seconds, we end up with the new resistance value of around 0.015 mbar min /mmµl, which is simply the previous value divided by 4² = 16. Remember, we doubled the length in each direction thus increasing the area by a factor of 4 reducing the pressure drop by a factor of 116 $\frac{1}{16}$ . If we evaluate the resistance value analytically using the exact equation (see Tab. 17.1), for rectangular channel cross-sections we end up with the same value of 0.015 mbar min /mm µl.

17.4.6 Concluding Remarks

As you can see, once a correlation between compactness C and correction α factor has been found, it is very easy to calculate the pressure drop for new channel geometries. This is a clear simplification especially for rectangular and square channel cross-sections where the pressure drop required the evaluation of a series (see Tab. 17.1). There are more microfluidic channel geometries for which similar simple correlations between C and α can be found. In this section we detailed the most important ones.

In general, the more sharp edges a cross-section has, the lower the expected flow velocities. This is due to the fact that the no-slip boundary conditions force the velocity profiles to rise only very slowly if two edges are in close proximity. Both edges retain the fluid; thus the overall fluid friction is higher and the gradient of the velocity profiles is less steep. Therefore, the overall velocities are smaller and therefore the associated pressure drop is higher. As a rule of thumb, the higher the boundary-to-area coefficient of a cross-section is the higher the influence of these retaining boundaries and the higher the associated pressure drops. Obviously, the circle is the best-case scenario, which is why the pressure drops associated with circular cross-section channels are smaller than for all other profiles.

17.5 Equivalent Circuit Theory

17.5.1 Introduction

In Eq. 17.23 we have derived the change of kinetic energy over time in a flowing fluid due to viscous dissipation. As you can see, this equation also corresponds to an equation known from electric circuits, namely the equation for Joule heating due to ohmic resistance.

Ohm’s Law. Ohm’s law¹ may well be considered the most elemental equation in electrical engineering. It is given by

R=UI $R = \frac{U}{I}$

(Eq. 17.42)

If the ohmic resistance of a physical body, e.g., a wire, is to be calculated, it is usually given by

R=ρlA $R = ρ \frac{l}{A}$

(Eq. 17.43)

where ρ is the specific resistance measured in Ω m, l is the length (e.g., of the wire), and A is the area, i.e., the cross-section of the wire. The specific resistance is often also given in form of its reciprocal value; this is referred to as conductivity κ=1ρ $κ = \frac{1}{ρ}$ , which is measured in siemens (S). For an electrolyte system, Eq. 17.43 can be used directly, interpreting l as the distance between the electrodes, A as the area of the electrodes, and ρ as the specific conductivity of the electrolyte.

Joule-Lenz Heating. Having introduced Ohm’s law we will quickly introduce the concept of Joule-Lenz heating. This phenomenon is observed in an electrical resistor through which an electrical current is flowing. The resistor will heat up over time. This is due to the fact that the flowing current, i.e., the moving electrons, will experience resistance in their movement due to the static atoms of the wire. This in terms will cause friction losses resulting in conversion of the energy of the electrons into heat. The higher the density of this flow the higher the resistance experienced. This phenomenon is often also referred to as Ohmic heating or resistive heating due to the fact that an Ohmic resistor is required. It was studied first by Joule² and later by Lenz³. Joule noted that the heat developed by a resistor when immersed in a water bath increased with the square root of the current and was proportional to the resistance and time, from which we derived

ΔT=cRI2Δt $Δ T = c R I^{2} Δ t$

As Joule measured the increase in temperature of the water, the density ρ and the heat capacity c_p of the water had to be taken into account, which accounts for the missing proportionality constant. Applying the first law of thermodynamics Eq. 6.32 we find

dEdEdtP=====mdUmcpdTRI2dtRI2RI2 $\begin{array}{l} d E & = & m d U \\ = & m c_{p} d T \\ = & R I^{2} d t \\ \frac{d E}{d t} & = & R I^{2} \\ P & = & R I^{2} \end{array}$

si83_e (Eq. 17.44)

=UI $\begin{array}{l} = & U I \end{array}$

(Eq. 17.45)

Today, we know the Joule-Lenz law of heating either in the form of Eq. 17.44 or Eq. 17.45.

Hydraulic Resistance. The analogy between electronics and fluid mechanics goes pretty far. Looking at Eq. 17.25 we see the resemblance of the pressure drop to Ohm’s law (see section 17.5.1):

RRhyd==UIΔpQ $\begin{array}{l} R & = & \frac{U}{I} \\ R_{hyd} & = & \frac{Δ p}{Q} \end{array}$

si85_e

where we have simplified the notation slightly in order to emphasize the resemblance. If we continue along these lines we can find the following equivalences between fluid mechanics and electronics:

• Hydraulic resistance R_hyd corresponds to ohmic resistance R.

• Volume V corresponds to charge Q.

• Flowrate Q (volume per time) corresponds to electric current I (charge per time).

• Pressure drop ∆p corresponds to voltage U (potential drop).

• Compliance, i.e., the change of volume due to pressure as a consequence of channel material elasticity corresponds to electrical capacitance C.

In the following we will briefly discuss series and parallel arrangements of hydraulic resistances further illustrating the similarity.

17.5.2 Serial Connection of Hydraulic Resistances

The first case is the connection of fluid resistances in series (see Fig. 17.8a). Here we obviously find

f17-08-9781455731411 — Fig. 17.8 Analogy between hydraulic and electrical resistance shown at a series of resistances (a) and two resistances in parallel (b).

Δp1Δp2==Rhyd,1Q1Rhyd,2Q2 $\begin{array}{l} Δ p_{1} & = & R_{hyd, 1} Q_{1} \\ Δ p_{2} & = & R_{hyd, 2} Q_{2} \end{array}$

whereby Q₁ = Q₂ = Q and ∆p₁ ≠ ∆p₂. The total pressure drop therefore is

Δptot.==Δp1+Δp2=(Rhyd,1+Rhyd,2)QRhyd,tot.Q $\begin{array}{l} Δ p_{tot .} & = & Δ p_{1} + Δ p_{2} = (R_{hyd, 1} + R_{hyd, 2}) Q \\ = & R_{hyd, tot .} Q \end{array}$

si87_e

As you can see, the two resistances simply add up to form the total resistance. This is the same behavior as expected for a series of electrical resistances.

17.5.3 Parallel Connection of Hydraulic Resistances

The next case we discuss is the parallel connection of hydraulic resistances (see Fig. 17.8b). Here we again find

Δp1Δp2==Rhyd,1Q1Rhyd,2Q2 $\begin{array}{l} Δ p_{1} & = & R_{hyd, 1} Q_{1} \\ Δ p_{2} & = & R_{hyd,2} Q_{2} \end{array}$

si88_e

whereby Q₁ + Q₂ = Q and ∆p₁ = ∆p₂ = ∆p. In this case we find

Δp(1Rhyd,1+1Rhyd,2)=Q $Δ p (\frac{1}{R_{hyd, 1}} + \frac{1}{R_{hyd,2}}) = Q$

si89_e

Δp==(1Rhyd,1+1Rhyd,2)−1QRhyd,1Rhyd,2Rhyd,1+Rhyd,2Q $\begin{array}{l} Δ p & = & {(\frac{1}{R_{hyd,1}} + \frac{1}{R_{hyd,2}})}^{- 1} Q \\ = & \frac{R_{hyd, 1} R_{hyd,2}}{R_{hyd, 1} + R_{hyd, 2}} Q \end{array}$

si90_e

which is obviously the expected behavior for a parallel arrangements of ohmic resistors as well.

17.5.4 Concluding Remarks

One important thing to keep in mind is that the analogy of hydraulic and electrical resistance has a certain limit. Whenever a microfluidic channel significantly widens in diameter or width the Poiseuille flow profile will expand or contract to fit this new channel cross-section. However, depending on the Reynolds number this expansion may lead to a contribution of the non-linear convection terms in the Navier-Stokes equation (see Eq. 11.40). If this is the case, the flow profile will not be due to a purely Poiseuille flow and the hydraulic resistance will increase. Therefore it is necessary to verify that for both cross-sections that the Reynolds numbers are small enough to neglect the contribution of the convection terms. The equations for calculating the hydraulic resistance of multiple hydraulic resistance in series or in parallel are only valid if the Reynolds number is Re 1 and the convection terms can be neglected.

17.6 Summary

In this section we have discussed the concept of hydraulic resistance which is an effect encountered whenever fluids are probed through technical systems. Hydraulic resistance gives rise to a pressure drop, i.e., a reduction of fluid pressure. The nature of this loss in pressure is due to the viscous friction within the fluid. Naturally, this effect is dependent on the shape of the channel geometry through which the fluid is probed. We have derived the hydraulic resistance for various commonly encountered flow channel geometries, e.g., circular or rectangular cross-sections. We have also detailed several simplifications that are commonly used to simplify the somewhat more complex equations that described the pressure drop.

..................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 17: Hydraulic Resistance

Create new playlist

Sign In

Sign Up