34.3 Implementation of a Stationary Flow Numerical Solver

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

34.2.3.3 Instationary Flows

Momentum Equations. The second scenario we will consider involves instationary flows. For instationary flows we have to solve Eq. 34.1, Eq. 34.2, and Eq. 34.3 using a suitable numerical scheme for the timesteps $\frac{\partial v_{x}}{\partial t}, \frac{\partial v_{y}}{\partial t}$ $\frac{\partial v_{x}}{\partial t}, \frac{\partial v_{y}}{\partial t}$ and $\frac{\partial v_{z}}{\partial t}$ $\frac{\partial v_{z}}{\partial t}$ ,. Usually, we would choose a forward Euler method (see section 27.2.2) for the time step. We could then rewrite Eq. 34.1, Eq. 34.2, and Eq. 34.3 as

$\begin{array}{l} \frac{\partial v_{x}}{\partial t} + v_{x} \frac{\partial v_{x}}{\partial x} + v_{y} \frac{\partial v_{x}}{\partial y} + v_{z} \frac{\partial v_{x}}{\partial z} = - \frac{1}{ρ} \frac{\partial p}{\partial x} + \frac{η}{ρ} (\frac{\partial^{2} v_{x}}{\partial x^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}}) + \frac{k_{x}}{ρ} \\ \frac{\partial v_{y}}{\partial t} + v_{x} \frac{\partial v_{y}}{\partial x} + v_{y} \frac{\partial v_{y}}{\partial y} + v_{z} \frac{\partial v_{y}}{\partial z} = - \frac{1}{ρ} \frac{\partial p}{\partial y} + \frac{η}{ρ} (\frac{\partial^{2} v_{y}}{\partial x^{2}} + \frac{\partial^{2} v_{y}}{\partial y^{2}} + \frac{\partial^{2} v_{y}}{\partial z^{2}}) + \frac{k_{y}}{ρ} \\ \frac{\partial v_{x}}{\partial t} + v_{x} \frac{\partial v_{z}}{\partial x} + v_{y} \frac{\partial v_{z}}{\partial y} + v_{z} \frac{\partial v_{z}}{\partial z} = - \frac{1}{ρ} \frac{\partial p}{\partial x} + \frac{η}{ρ} (\frac{\partial^{2} v_{z}}{\partial x^{2}} + \frac{\partial^{2} v_{z}}{\partial y^{2}} + \frac{\partial^{2} v_{z}}{\partial z^{2}}) + k_{z} \end{array}$ $\begin{array}{l} \frac{\partial v_{x}}{\partial t} + v_{x} \frac{\partial v_{x}}{\partial x} + v_{y} \frac{\partial v_{x}}{\partial y} + v_{z} \frac{\partial v_{x}}{\partial z} = - \frac{1}{ρ} \frac{\partial p}{\partial x} + \frac{η}{ρ} (\frac{\partial^{2} v_{x}}{\partial x^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}}) + \frac{k_{x}}{ρ} \\ \frac{\partial v_{y}}{\partial t} + v_{x} \frac{\partial v_{y}}{\partial x} + v_{y} \frac{\partial v_{y}}{\partial y} + v_{z} \frac{\partial v_{y}}{\partial z} = - \frac{1}{ρ} \frac{\partial p}{\partial y} + \frac{η}{ρ} (\frac{\partial^{2} v_{y}}{\partial x^{2}} + \frac{\partial^{2} v_{y}}{\partial y^{2}} + \frac{\partial^{2} v_{y}}{\partial z^{2}}) + \frac{k_{y}}{ρ} \\ \frac{\partial v_{x}}{\partial t} + v_{x} \frac{\partial v_{z}}{\partial x} + v_{y} \frac{\partial v_{z}}{\partial y} + v_{z} \frac{\partial v_{z}}{\partial z} = - \frac{1}{ρ} \frac{\partial p}{\partial x} + \frac{η}{ρ} (\frac{\partial^{2} v_{z}}{\partial x^{2}} + \frac{\partial^{2} v_{z}}{\partial y^{2}} + \frac{\partial^{2} v_{z}}{\partial z^{2}}) + k_{z} \end{array}$

si123_e

from which we find

$\frac{\partial v_{x}}{\partial t} = - \frac{1}{ρ} \frac{\partial p}{\partial x} + \frac{η}{ρ} (\frac{\partial^{2} v_{x}}{\partial x^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}}) + \frac{k_{x}}{ρ} - (v_{x} \frac{\partial v_{x}}{\partial x} + v_{y} \frac{\partial v_{x}}{\partial y} + v_{z} \frac{\partial v_{x}}{\partial z})$ $\frac{\partial v_{x}}{\partial t} = - \frac{1}{ρ} \frac{\partial p}{\partial x} + \frac{η}{ρ} (\frac{\partial^{2} v_{x}}{\partial x^{2}} + \frac{\partial^{2} v_{x}}{\partial y^{2}} + \frac{\partial^{2} v_{x}}{\partial z^{2}}) + \frac{k_{x}}{ρ} - (v_{x} \frac{\partial v_{x}}{\partial x} + v_{y} \frac{\partial v_{x}}{\partial y} + v_{z} \frac{\partial v_{x}}{\partial z})$

si124_e (Eq. 34.39)

$\frac{\partial v_{y}}{\partial t} = - \frac{1}{ρ} \frac{\partial p}{\partial y} + \frac{η}{ρ} (\frac{\partial^{2} v_{y}}{\partial x^{2}} + \frac{\partial^{2} v_{y}}{\partial y^{2}} + \frac{\partial^{2} v_{y}}{\partial z^{2}}) + \frac{k_{y}}{ρ} - (v_{x} \frac{\partial v_{y}}{\partial x} + v_{y} \frac{\partial v_{y}}{\partial y} + v_{z} \frac{\partial v_{y}}{\partial z})$ $\frac{\partial v_{y}}{\partial t} = - \frac{1}{ρ} \frac{\partial p}{\partial y} + \frac{η}{ρ} (\frac{\partial^{2} v_{y}}{\partial x^{2}} + \frac{\partial^{2} v_{y}}{\partial y^{2}} + \frac{\partial^{2} v_{y}}{\partial z^{2}}) + \frac{k_{y}}{ρ} - (v_{x} \frac{\partial v_{y}}{\partial x} + v_{y} \frac{\partial v_{y}}{\partial y} + v_{z} \frac{\partial v_{y}}{\partial z})$

si125_e (Eq. 34.40)

$\frac{\partial v_{z}}{\partial t} = - \frac{1}{ρ} \frac{\partial p}{\partial z} + \frac{η}{ρ} (\frac{\partial^{2} v_{z}}{\partial x^{2}} + \frac{\partial^{2} v_{z}}{\partial y^{2}} + \frac{\partial^{2} v_{z}}{\partial z^{2}}) + \frac{k_{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})$ $\frac{\partial v_{z}}{\partial t} = - \frac{1}{ρ} \frac{\partial p}{\partial z} + \frac{η}{ρ} (\frac{\partial^{2} v_{z}}{\partial x^{2}} + \frac{\partial^{2} v_{z}}{\partial y^{2}} + \frac{\partial^{2} v_{z}}{\partial z^{2}}) + \frac{k_{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})$

si126_e (Eq. 34.41)

Here, we use the same numerical schemes as for the stationary flows, adding the forward Euler schemes for the timesteps:

$\begin{array}{l} v_{x}^{(t + 1, x, y, z)} = v_{x}^{(t, x, y, z)} + h_{t} {\frac{\partial v_{x}}{\partial t}}^{(t, x, y, z)} \\ v_{y}^{(t + 1, x, y, z)} = v_{y}^{(t, x, y, z)} + h_{t} {\frac{\partial v_{x}}{\partial t}}^{(t, x, y, z)} \\ v_{z}^{(t + 1, x, y, z)} = v_{z}^{(t, x, y, z)} + h_{t} {\frac{\partial v_{x}}{\partial t}}^{(t, x, y, z)} \end{array}$ $\begin{array}{l} v_{x}^{(t + 1, x, y, z)} = v_{x}^{(t, x, y, z)} + h_{t} {\frac{\partial v_{x}}{\partial t}}^{(t, x, y, z)} \\ v_{y}^{(t + 1, x, y, z)} = v_{y}^{(t, x, y, z)} + h_{t} {\frac{\partial v_{x}}{\partial t}}^{(t, x, y, z)} \\ v_{z}^{(t + 1, x, y, z)} = v_{z}^{(t, x, y, z)} + h_{t} {\frac{\partial v_{x}}{\partial t}}^{(t, x, y, z)} \end{array}$

si127_e

In contrast to the stationary flow cases, where we need to iterate the equations until we find no more changes in the dependent variables, we do not iterate the stationary solutions. However, we need to include the time index t since instationary solutions obviously need to step in time. Here, the index t refers to the current time step, whereas t + 1 refers to the following time step. We then find the following from Eq. 34.39, Eq. 34.40, and Eq. 34.41:

$\begin{array}{l} v_{x}^{(t + 1, x, y, z)} = v_{x}^{(t, x, y, z)} - \frac{h_{t}}{ρ h_{x y z}} \frac{p^{(t, x + 1, y, z)} - p^{(t, x - 1, y, z)}}{2} + \frac{h_{t} η}{h_{x y z}^{2} ρ} (\begin{array}{l} v_{x}^{(t, x + 1, y, z)} + v_{x}^{(t, x - 1, y, z)} + v_{x}^{(t, x, y + 1, z)} + \\ v_{x}^{(t, x, y - 1, z)} + v_{x}^{(t, x, y, z + 1)} + v_{x}^{(t, x, y, z - 1)} - 6 v_{x}^{(t, x, y, z)} + \frac{h_{t}}{ρ} k_{x} \end{array}) \\ - \frac{h_{t}}{h_{x y z}} (\begin{array}{l} v_{x}^{(t, x, y, z)} \frac{v_{x}^{(t, x + 1, y, z)} - v_{x}^{(t, x - 1, y, z)}}{2} + v_{y}^{(t, x, y, z)} \frac{v_{x}^{(t, x, y + 1, z)} - v_{x}^{(t, x, y - 1, z)}}{2} + \\ v_{z}^{(t, x, y, z)} \frac{v_{x}^{(t, x, y, z + 1)} - v_{x}^{(t, x, y, z - 1)}}{2} \end{array}) \end{array}$ $\begin{array}{l} v_{x}^{(t + 1, x, y, z)} = v_{x}^{(t, x, y, z)} - \frac{h_{t}}{ρ h_{x y z}} \frac{p^{(t, x + 1, y, z)} - p^{(t, x - 1, y, z)}}{2} + \frac{h_{t} η}{h_{x y z}^{2} ρ} (\begin{array}{l} v_{x}^{(t, x + 1, y, z)} + v_{x}^{(t, x - 1, y, z)} + v_{x}^{(t, x, y + 1, z)} + \\ v_{x}^{(t, x, y - 1, z)} + v_{x}^{(t, x, y, z + 1)} + v_{x}^{(t, x, y, z - 1)} - 6 v_{x}^{(t, x, y, z)} + \frac{h_{t}}{ρ} k_{x} \end{array}) \\ - \frac{h_{t}}{h_{x y z}} (\begin{array}{l} v_{x}^{(t, x, y, z)} \frac{v_{x}^{(t, x + 1, y, z)} - v_{x}^{(t, x - 1, y, z)}}{2} + v_{y}^{(t, x, y, z)} \frac{v_{x}^{(t, x, y + 1, z)} - v_{x}^{(t, x, y - 1, z)}}{2} + \\ v_{z}^{(t, x, y, z)} \frac{v_{x}^{(t, x, y, z + 1)} - v_{x}^{(t, x, y, z - 1)}}{2} \end{array}) \end{array}$

si128_e (Eq. 34.42)

$\begin{array}{l} v_{y}^{(t + 1, x, y, z)} = v_{y}^{(t, x, y, z)} - \frac{h_{t}}{ρ h_{x y z}} \frac{p^{(t, x + 1, y, z)} - p^{(t, x - 1, y, z)}}{2} + \frac{h_{t} η}{h_{x y z}^{2} ρ} (\begin{array}{l} v_{y}^{(t, x + 1, y, z)} + v_{y}^{(t, x - 1, y, z)} + v_{y}^{(t, x, y + 1, z)} + \\ v_{y}^{(t, x, y - 1, z)} + v_{y}^{(t, x, y, z + 1)} + v_{y}^{(t, x, y, z - 1)} - 6 v_{y}^{(t, x, y, z)} + \frac{h_{t}}{ρ} k_{y} \end{array}) \\ - \frac{h_{t}}{h_{x y z}} (\begin{array}{l} v_{x}^{(t, x, y, z)} \frac{v_{y}^{(t, x + 1, y, z)} - v_{y}^{(t, x - 1, y, z)}}{2} + v_{y}^{(t, x, y, z)} \frac{v_{y}^{(t, x, y + 1, z)} - v_{y}^{(t, x, y - 1, z)}}{2} + \\ v_{z}^{(t, x, y, z)} \frac{v_{y}^{(t, x, y, z + 1)} - v_{y}^{(t, x, y, z - 1)}}{2} \end{array}) \end{array}$ $\begin{array}{l} v_{y}^{(t + 1, x, y, z)} = v_{y}^{(t, x, y, z)} - \frac{h_{t}}{ρ h_{x y z}} \frac{p^{(t, x + 1, y, z)} - p^{(t, x - 1, y, z)}}{2} + \frac{h_{t} η}{h_{x y z}^{2} ρ} (\begin{array}{l} v_{y}^{(t, x + 1, y, z)} + v_{y}^{(t, x - 1, y, z)} + v_{y}^{(t, x, y + 1, z)} + \\ v_{y}^{(t, x, y - 1, z)} + v_{y}^{(t, x, y, z + 1)} + v_{y}^{(t, x, y, z - 1)} - 6 v_{y}^{(t, x, y, z)} + \frac{h_{t}}{ρ} k_{y} \end{array}) \\ - \frac{h_{t}}{h_{x y z}} (\begin{array}{l} v_{x}^{(t, x, y, z)} \frac{v_{y}^{(t, x + 1, y, z)} - v_{y}^{(t, x - 1, y, z)}}{2} + v_{y}^{(t, x, y, z)} \frac{v_{y}^{(t, x, y + 1, z)} - v_{y}^{(t, x, y - 1, z)}}{2} + \\ v_{z}^{(t, x, y, z)} \frac{v_{y}^{(t, x, y, z + 1)} - v_{y}^{(t, x, y, z - 1)}}{2} \end{array}) \end{array}$

si129_e (Eq. 34.43)

$\begin{array}{l} v_{z}^{(t + 1, x, y, z)} = v_{z}^{(t, x, y, z)} - \frac{h_{t}}{ρ h_{x y z}} \frac{p^{(t, x + 1, y, z)} - p^{(t, x - 1, y, z)}}{2} + \frac{h_{t} η}{h_{x y z}^{2} ρ} (\begin{array}{l} v_{z}^{(t, x + 1, y, z)} + v_{z}^{(t, x - 1, y, z)} + v_{z}^{(t, x, y + 1, z)} + \\ v_{z}^{(t, x, y - 1, z)} + v_{z}^{(t, x, y, z + 1)} + v_{z}^{(t, x, y, z - 1)} - 6 v_{z}^{(t, x, y, z)} + \frac{h_{t}}{ρ} k_{z} \end{array}) \\ - \frac{h_{t}}{h_{x y z}} (\begin{array}{l} v_{z}^{(t, x, y, z)} \frac{v_{z}^{(t, x + 1, y, z)} - v_{z}^{(t, x - 1, y, z)}}{2} + v_{z}^{(t, x, y, z)} \frac{v_{z}^{(t, x, y + 1, z)} - v_{z}^{(t, x, y - 1, z)}}{2} + \\ v_{z}^{(t, x, y, z)} \frac{v_{z}^{(t, x, y, z + 1)} - v_{z}^{(t, x, y, z - 1)}}{2} \end{array}) \end{array}$ $\begin{array}{l} v_{z}^{(t + 1, x, y, z)} = v_{z}^{(t, x, y, z)} - \frac{h_{t}}{ρ h_{x y z}} \frac{p^{(t, x + 1, y, z)} - p^{(t, x - 1, y, z)}}{2} + \frac{h_{t} η}{h_{x y z}^{2} ρ} (\begin{array}{l} v_{z}^{(t, x + 1, y, z)} + v_{z}^{(t, x - 1, y, z)} + v_{z}^{(t, x, y + 1, z)} + \\ v_{z}^{(t, x, y - 1, z)} + v_{z}^{(t, x, y, z + 1)} + v_{z}^{(t, x, y, z - 1)} - 6 v_{z}^{(t, x, y, z)} + \frac{h_{t}}{ρ} k_{z} \end{array}) \\ - \frac{h_{t}}{h_{x y z}} (\begin{array}{l} v_{z}^{(t, x, y, z)} \frac{v_{z}^{(t, x + 1, y, z)} - v_{z}^{(t, x - 1, y, z)}}{2} + v_{z}^{(t, x, y, z)} \frac{v_{z}^{(t, x, y + 1, z)} - v_{z}^{(t, x, y - 1, z)}}{2} + \\ v_{z}^{(t, x, y, z)} \frac{v_{z}^{(t, x, y, z + 1)} - v_{z}^{(t, x, y, z - 1)}}{2} \end{array}) \end{array}$

si130_e (Eq. 34.44)

As in the case of stationary flows, we solve Eq. 34.42, Eq. 34.43, and Eq. 34.44 one by one, updating the values after each equation.

Pressure-Poisson Equation. Having updated these values, we then need to calculate the corrected pressure. However, as we have done in case of the stationary flows, we need to be careful to use the correct velocity values when updating the pressure. If you recall, we developed the pressure-Poisson equation using the divergence of the Navier–Stokes equation. In stationary flows, there is only one velocity field. Although we use indices during the iteration (namely i and i + 1) we really only consider one velocity field for which we want to find a converging solution. However, in the case of instationary flows, there are effectively two velocity fields: the field at time step t and the field at time step t + 1. As stated in the case of the stationary flow, we want to update the pressure for the updated velocity field, i.e., for t + 1. For this, we use the pressure-Poisson equation for the instationary case, i.e.,Eq. 34.30. Compared to the pressure-Poisson equation for the stationary case, there are a couple of important things to note.

In the instationary case, we cannot invoke the continuity equation directly when stepping in time. If you look at Eq. 34.30 you will see that the equation is currently written such that it gives us the Laplacian of the pressure field at the time point t. If we use the equation in this form, we adapt the pressure field such that continuity is fulfilled at point t. Therefore, all terms involving $\vec{\nabla} {\vec{v}}^{(t)}$ $\vec{\nabla} {\vec{v}}^{(t)}$ could then be set to zero due to the continuity. However, this would not help us with our time step. We have written Eq. 34.42, Eq. 34.43, and Eq. 34.44 such that they return the velocities at time point t + 1. We should therefore rewrite Eq. 34.30 such that it returns the Laplacian of the pressure field at time point t + 1. This is actually very simple; we just exchange the index on the left-hand side, being left with

$Δ_{p}^{(t + 1)} = \frac{ρ}{h_{t}} (\vec{\nabla} {\vec{v}}^{(t)} - \vec{\nabla} {\vec{v}}^{(t + 1)}) - ρ \vec{\nabla} (({\vec{v}}^{(t)} \cdot \vec{\nabla}) {\vec{v}}^{(t)}) + \vec{\nabla} {\vec{k}}^{(t)} + η Δ (\vec{\nabla} {\vec{v}}^{(t)})$ $Δ_{p}^{(t + 1)} = \frac{ρ}{h_{t}} (\vec{\nabla} {\vec{v}}^{(t)} - \vec{\nabla} {\vec{v}}^{(t + 1)}) - ρ \vec{\nabla} (({\vec{v}}^{(t)} \cdot \vec{\nabla}) {\vec{v}}^{(t)}) + \vec{\nabla} {\vec{k}}^{(t)} + η Δ (\vec{\nabla} {\vec{v}}^{(t)})$

si132_e (Eq. 34.45)

However, if we do this, we cannot use the continuity equation for the time point t since all values we calculated (all three velocities and the pressure) are for time point t + 1. However, we can then use the continuity equation to set all terms $\vec{\nabla} {\vec{v}}^{(t)}$ $\vec{\nabla} {\vec{v}}^{(t)}$ to zero. This is in fact only one term, the second term in the first bracket of the right-hand side. We then simplify Eq. 34.45 to

$Δ_{p}^{(t + 1)} = \frac{ρ}{h_{t}} \vec{\nabla} {\vec{v}}^{(t)} - ρ \vec{\nabla} (({\vec{v}}^{(t)} \cdot \vec{\nabla}) {\vec{v}}^{(t)}) + \vec{\nabla} {\vec{k}}^{(t)} + η Δ (\vec{\nabla} {\vec{v}}^{(t)})$ $Δ_{p}^{(t + 1)} = \frac{ρ}{h_{t}} \vec{\nabla} {\vec{v}}^{(t)} - ρ \vec{\nabla} (({\vec{v}}^{(t)} \cdot \vec{\nabla}) {\vec{v}}^{(t)}) + \vec{\nabla} {\vec{k}}^{(t)} + η Δ (\vec{\nabla} {\vec{v}}^{(t)})$

si134_e (Eq. 34.46)

We now introduce a second simplification by using the updated values $v_{x}^{(t + 1)}, v_{y}^{(t + 1)}$ $v_{x}^{(t + 1)}, v_{y}^{(t + 1)}$ , and $v_{z}^{(t + 1)}$ $v_{z}^{(t + 1)}$ for calculating the diffusion term of the right-hand side. Doing so we obtain

$Δ_{p}^{(t + 1)} = \frac{ρ}{h_{t}} \vec{\nabla} {\vec{v}}^{(t)} - ρ \vec{\nabla} (({\vec{v}}^{(t)} \cdot \vec{\nabla}) {\vec{v}}^{(t)}) + \vec{\nabla} {\vec{k}}^{(t)} + η Δ (\vec{\nabla} {\vec{v}}^{(t + 1)})$ $Δ_{p}^{(t + 1)} = \frac{ρ}{h_{t}} \vec{\nabla} {\vec{v}}^{(t)} - ρ \vec{\nabla} (({\vec{v}}^{(t)} \cdot \vec{\nabla}) {\vec{v}}^{(t)}) + \vec{\nabla} {\vec{k}}^{(t)} + η Δ (\vec{\nabla} {\vec{v}}^{(t + 1)})$

si137_e

which allows us to drop the diffusion term due to the continuity equation $(\vec{\nabla} {\vec{v}}^{(t + 1)} = 0)$ $(\vec{\nabla} {\vec{v}}^{(t + 1)} = 0)$ in which case we are left with

$Δ_{p}^{(t + 1)} = \frac{ρ}{h_{t}} \vec{\nabla} {\vec{v}}^{(t)} + \vec{\nabla} {\vec{k}}^{(t)} - ρ \vec{\nabla} ({\vec{v}}^{(t)} \cdot \vec{\nabla}) {\vec{v}}^{(t)}$ $Δ_{p}^{(t + 1)} = \frac{ρ}{h_{t}} \vec{\nabla} {\vec{v}}^{(t)} + \vec{\nabla} {\vec{k}}^{(t)} - ρ \vec{\nabla} ({\vec{v}}^{(t)} \cdot \vec{\nabla}) {\vec{v}}^{(t)}$

si139_e (Eq. 34.47)

This is as far as we can simplify the pressure-Poisson equation. If you compare Eq. 34.47 with the pressure-Poisson equation found for the stationary case (see Eq. 34.24) you will see that they are identical, apart from the time-dependent term $\frac{ρ}{h_{t}} \vec{\nabla} {\vec{v}}^{(t)}$ $\frac{ρ}{h_{t}} \vec{\nabla} {\vec{v}}^{(t)}$ .

We can therefore simply copy the numerical scheme from the stationary pressure-Poisson equation (see Eq. 34.38) and add a suitable numerical scheme for the instationary term. For this contribution we use the following scheme:

$\vec{\nabla} {\vec{v}}^{(t)} = \frac{1}{h_{x y z}} \frac{v_{x}^{(x + 1, y, z)} - v_{x}^{(x - 1, y, z)} + v_{y}^{(x, y + 1, z)} - v_{y}^{(x, y - 1, z)} + v_{z}^{(x, y, z + 1)} - v_{z}^{(x, y, z - 1)}}{2}$ $\vec{\nabla} {\vec{v}}^{(t)} = \frac{1}{h_{x y z}} \frac{v_{x}^{(x + 1, y, z)} - v_{x}^{(x - 1, y, z)} + v_{y}^{(x, y + 1, z)} - v_{y}^{(x, y - 1, z)} + v_{z}^{(x, y, z + 1)} - v_{z}^{(x, y, z - 1)}}{2}$

si141_e

in which case we find the following numerical scheme:

$\begin{array}{l} p^{(t + 1, x, y, z)} = \frac{h_{x y z}}{6 h_{t}} \frac{v_{x}^{(t, x + 1, y, z)} - v_{x}^{(t, x - 1, y, z)} + v_{y}^{(t, x, y + 1, z)} - v_{y}^{(t, x, y - 1, z)} + v_{z}^{(t, x, y, z + 1)} - v_{z}^{(t, x, y, z - 1)}}{2} - \\ \frac{h_{x y z}}{6} \frac{k_{x}^{(x + 1, y, z)} - k_{x}^{(x - 1, y, z)} + k_{y}^{(x, y + 1, z)} - k_{y}^{(x, y - 1, z)} + k_{z}^{(x, y, z + 1)} - k_{z}^{(x, y, z - 1)}}{2} + \\ \frac{ρ}{6} (\begin{array}{l} {(\frac{v_{x}^{(i, x + 1, y, z)} - v_{x}^{(i, x - 1, y, z)}}{2})}^{2} + {(\frac{v_{y}^{(i, x, y + 1, z)} - v_{y}^{(i, x, y - 1, z)}}{2})}^{2} + \\ {(\frac{v_{z}^{(i, x, y, z + 1)} - v_{z}^{(i, x, y, z - 1)}}{2})}^{2} + 2 (\begin{array}{l} (\frac{v_{y}^{(i, x + 1, y, z)} - v_{y}^{(i, x - 1, y, z)}}{2}) (\frac{v_{x}^{(i, x, y + 1, z)} - v_{x}^{(i, x, y - 1, z)}}{2}) + \\ (\frac{v_{x}^{(i, x, y, z + 1)} - v_{x}^{(i, x, y, z - 1)}}{2}) (\frac{v_{z}^{(i, x + 1, y, z)} - v_{z}^{(i, x - 1, y, z)}}{2}) + \\ (\frac{v_{y}^{(i, x, y, z + 1)} - v_{y}^{(i, x, y, z - 1)}}{2}) (\frac{v_{z}^{(i, x, y + 1, z)} - v_{z}^{(i, x, y - 1, z)}}{2}) \end{array}) \end{array}) \\ + \frac{1}{6} (p^{(i, x + 1, y, z)} + p^{(i, x - 1, y, z)} + p^{(i, x, y + 1, z)} + p^{(i, x, y - 1, z)} + p^{(i, x, y, z + 1)} + p^{(i, x, y, z - 1)}) \end{array}$ $\begin{array}{l} p^{(t + 1, x, y, z)} = \frac{h_{x y z}}{6 h_{t}} \frac{v_{x}^{(t, x + 1, y, z)} - v_{x}^{(t, x - 1, y, z)} + v_{y}^{(t, x, y + 1, z)} - v_{y}^{(t, x, y - 1, z)} + v_{z}^{(t, x, y, z + 1)} - v_{z}^{(t, x, y, z - 1)}}{2} - \\ \frac{h_{x y z}}{6} \frac{k_{x}^{(x + 1, y, z)} - k_{x}^{(x - 1, y, z)} + k_{y}^{(x, y + 1, z)} - k_{y}^{(x, y - 1, z)} + k_{z}^{(x, y, z + 1)} - k_{z}^{(x, y, z - 1)}}{2} + \\ \frac{ρ}{6} (\begin{array}{l} {(\frac{v_{x}^{(i, x + 1, y, z)} - v_{x}^{(i, x - 1, y, z)}}{2})}^{2} + {(\frac{v_{y}^{(i, x, y + 1, z)} - v_{y}^{(i, x, y - 1, z)}}{2})}^{2} + \\ {(\frac{v_{z}^{(i, x, y, z + 1)} - v_{z}^{(i, x, y, z - 1)}}{2})}^{2} + 2 (\begin{array}{l} (\frac{v_{y}^{(i, x + 1, y, z)} - v_{y}^{(i, x - 1, y, z)}}{2}) (\frac{v_{x}^{(i, x, y + 1, z)} - v_{x}^{(i, x, y - 1, z)}}{2}) + \\ (\frac{v_{x}^{(i, x, y, z + 1)} - v_{x}^{(i, x, y, z - 1)}}{2}) (\frac{v_{z}^{(i, x + 1, y, z)} - v_{z}^{(i, x - 1, y, z)}}{2}) + \\ (\frac{v_{y}^{(i, x, y, z + 1)} - v_{y}^{(i, x, y, z - 1)}}{2}) (\frac{v_{z}^{(i, x, y + 1, z)} - v_{z}^{(i, x, y - 1, z)}}{2}) \end{array}) \end{array}) \\ + \frac{1}{6} (p^{(i, x + 1, y, z)} + p^{(i, x - 1, y, z)} + p^{(i, x, y + 1, z)} + p^{(i, x, y - 1, z)} + p^{(i, x, y, z + 1)} + p^{(i, x, y, z - 1)}) \end{array}$

si142_e (Eq. 34.48)

which completes the numerical scheme.

34.3 Implementation of a Stationary Flow Numerical Solver

34.3.1 Introduction

Having derived the numerical schemes for both stationary and instationary flows we will now proceed to implementing one of the schemes, thus creating a solver capable of solving fully three-dimensional flow problems. We will implement the stationary solver, which is the more complicated one due to the fact that it requires iteration. For this, we require Eq. 34.34, Eq. 34.35, Eq. 34.36, and Eq. 34.38.

Compared to the solvers we have been implementing so far, there are a couple of changes we need to account for.

34.3.2 Boundary Values

Until now, we have dealt with (mainly) Dirichlet boundary conditions. Wherever we encountered an open end on the domain, we assumed Neumann boundary conditions. As an example, refer to listing 33.8, where we simply copied the values of the second-to-last index of t to the last index. Doing so, we obtained Neumann-type boundary values on the end of the domain.

However, for three-dimensional flows, there is no “left” or “right” side of the domain and there is also no preferred calculation direction. We have to account for the boundary conditions in a more controlled fashion. Therefore, we will upgrade the boolean array BC, which we have previously used for the boundary values to a byte array. There are three values which we account for:

• 0 indicates that the cell has free values; these values must be updated during calculation.

• − 1 indicates a Neumann boundary condition; this value is updated during the calculation to the value of the neighboring cell.

• 0 indicates a Dirichlet boundary condition; this value is not updated since it is defined by the initial values.

The second modification we need to implement relates to the fact that we need to calculate velocities and pressure during the calculation. This also means that we require two sets of boundary conditions: one for the velocities (referred to as BCvel) and for the pressure (referred to as BCpress) as we can have boundary values for one but not for the other. Some typical scenarios are discussed in the following.

Open Boundaries. Open boundaries usually imply that the gradient in both pressure and velocity must be zero. Therefore, we would set Neumann-type boundary conditions for both the pressure and the velocity at these cells.

Applied Pressure Gradient. The classical case of pressure-driven flow, i.e., Poisson flow (see section 15.4), is generated by setting a given pressure at one end of the computation domain and a lower value at the other end. As we will see, our solver will find out on its own that the pressure drops linearly over a straight channel (see section 34.4.6). As we usually do not know the velocities in these cases, we simply apply Neumann-type boundary conditions for the cells for which we set the pressure.

Channel Wall: No-Slip Boundary Condition. A no-slip boundary condition is generated by applying Dirichlet boundary conditions for the velocity and setting the velocity to zero at these cells. For the pressure, we set Neumann-type boundary conditions.

Open Surface. An open surface can be generated by setting a Dirichlet boundary condition with a given value for the pressure. We do not have to set a boundary condition for the velocity, which can be defined by the solver.

34.3.3 Constants

Before doing the actual implementation we will modify Eq. 34.34, Eq. 34.35, Eq. 34.36, and Eq. 34.38, summing up the constants and correcting for the units.

When implementing our numerical scheme, we should make use of every simplification we can potentially come up with in order to make the code run faster. As you will see, several of the flow cases we calculate require thousands of iterations. If a calculation can be sped up by only a couple of instructions we will see a visible overall improvement in speed. Precalculating these factors is one of the more obvious ways of speeding up the calculation.

Desired Units. Before calculating and correcting the constants for the correct units, we obviously have to define the units first. We will use the following units:

• velocities v_x, v_y, v_z or velx, vely, velz: mms^− 1

• pressure p or press: mbar

• viscosity η or viscosity: mPa s

• density ρ or density: g cm^− 3

• volume forces k_x, k_y, k_z or kvolx, kvoly, kvolz: mN cm^− 3

• stepwidth in space h_xyz or hxyz: μm

The constants will take care of the correct units so that the result of Eq. 34.34, Eq. 34.35, and Eq. 34.36 will be given in mm s^− 1 and the result of Eq. 34.38 in mbar.

SOR Constants. The first constants we require are due to the fact that we require an iterative scheme. We will again implement a SOR scheme (see section 25.3.4) with a given relaxation parameter ω. In order to avoid recalculating these factors for every iteration we use the two constants omega1 and omega2, which are defined as

omega1 = 1 − omega [−]
omega2 = omega [−]

where omega is the user-supplied relaxation parameter ω.

Denominator Factor for the Velocities. If you refer to Eq. 34.34, Eq. 34.35, and Eq. 34.36 you can see that there is a constant in all of the denominators. This value is given by

$denominatorconstant = 1000 \cdot 6 \frac{η}{ρ h_{x y z}} [mm s^{- 1}]$ $denominatorconstant = 1000 \cdot 6 \frac{η}{ρ h_{x y z}} [mm s^{- 1}]$

where we have introduced the constant 1000 to adjust for the unit mms^− 1.

Constant for the Pressure Gradient in the Velocities. The terms originating from the pressure gradient in Eq. 34.34, Eq. 34.35, and Eq. 34.36 are preceded by a constant that is given by

$pressureconstantVel = \frac{100000}{ρ} [m m^{2} s^{- 2}]$ $pressureconstantVel = \frac{100000}{ρ} [m m^{2} s^{- 2}]$

which we have corrected for the desired unit mm² s^− 2.

Laplace Constant for the Velocities. The terms originating from the Laplace operator, i.e., due to momentum diffusion, are preceded by a constant in Eq. 34.34, Eq. 34.35, and Eq. 34.36. It is given by

$laplaceconstantVel = \frac{1000 η}{ρ h_{x y z}} [mm s^{- 1}]$ $laplaceconstantVel = \frac{1000 η}{ρ h_{x y z}} [mm s^{- 1}]$

which we have again corrected in order to obtain the correct unit, which is mm s^− 1.

Volume Force Constant for the Velocities. The velocity contribution due to the action of the volume forces is also preceded by a constant in Eq. 34.34, Eq. 34.35, and Eq. 34.36 which is given by

$volumeforceconstantVel = \frac{h_{x y z}}{ρ} [m m^{2} s^{- 2}]$ $volumeforceconstantVel = \frac{h_{x y z}}{ρ} [m m^{2} s^{- 2}]$

This factor requires no additional constant since the units add up correctly.

Volume Force Constants for the Pressure. Volume forces also give rise to a contribution in Eq. 34.38, where they are preceded by the constant

$volumeforceconstantVel = \frac{h_{x y z}}{100000} [mbar]$ $volumeforceconstantVel = \frac{h_{x y z}}{100000} [mbar]$

They have to be corrected by the constant 100 000 in order to end up with the correct unit mbar.

Convection Constant for the Pressure. The terms involving the velocity gradients, i.e., the contribution of the convection, are preceded by the following constant in Eq. 34.38:

$convectionconstantPress = \frac{ρ}{6 \cdot 100000} [mbar]$ $convectionconstantPress = \frac{ρ}{6 \cdot 100000} [mbar]$

which has to be corrected to obtain the desired unit mbar.

Pressure Constant for the Pressure. The last constant is simply given by

$pressureconstantPress = \frac{1}{6} [-]$ $pressureconstantPress = \frac{1}{6} [-]$

Obviously, this constant has no unit.

34.3.4 Structures

In the next step, we will define two structures in C, which will make the implementation of the solver slightly easier. The structures are given in listing 34.1.

Listing 34.1

[C] Helper structures used for the implementation of the three-dimensional solver. This listing is part of the NeptunLib and can be found in the file NeptunLib.h in the supplementary material.

1 typedef struct{
2  double p,m,d;
3 }currentflow3Ddirection;
4
5 typedef struct{
6  double value;
7  double* data;
8  currentflow3Ddirection x,y,z;
9 }currentflow3Dvariable;

The first structure is termed currentflow3Ddirection. It holds three variable p (for “plus”), m (for “minus”), and d (for “delta”). The structure currentflow3Ddirection is used in the second structure shown in listing 34.1, which is currentflow3Dvariable. Within our data structure a single currentflow3Dvariable holds information about the value of a dependent variable and its changes along each of the three axes. As an example, one of the dependent variables of our flow which we are interested in is the variable v_x, i.e., the x-component of the velocity $\vec{v}$ $\vec{v}$ . At a given location, i.e., for given values of x, y, and z, the cell v_x^(x,y,z) has a given numerical value. This value is set as the variable value in the structure currentflow3Dvariable. Besides the value at (x, y, z), we also require the changes of the dependent variable along the three axes x, y, and z. These are given as the variables x, y, and z in the structure currentflow3Dvariable. Each of these variables is a structure currentflow3Ddirection containing the value of the neighboring cells in the respective positive direction p, the negative direction m, and the change d, which is simply $\frac{p - m}{2}$ $\frac{p - m}{2}$ . Keep in mind that we are using central finite difference schemes, hence the method of calculating the change along the given axis.

Returning to our example, v_x^(x,y,z) would be stored in value, $v_{x}^{(x + 1, y, z)}$ $v_{x}^{(x + 1, y, z)}$ would be stored in x.p, v_x^{(x ‐ 1,y,z)} would be stored in x.m, $v_{x}^{(x, y + 1, z)}$ $v_{x}^{(x, y + 1, z)}$ would be stored in y.p, and so on. The structure currentflow3Dvariable also holds a pointer reference to the numerical data of the respective dependent variable. We need this reference to access the neighboring values.

34.3.5 Implementation

34.3.5.1 Main Listing

We are now ready to implement the three-dimensional solver. The listing is shown in listing 34.2. The function is called SolveFlow3D. It expects a number of arguments. The length L is the number of cells along the x-direction. The width W is the number of cells along the y-direction. The height H is the number of cells along the z-direction. The stepwidth h_xyz is given as the argument hxyz in μm. The density ρ is given as the argument density in gcm^− 3. The viscosity η is given as the argument η in mPas. The desired precision for the velocities is given as epsilonVel in mm s^− 1 and for the pressure as epsilonPress in mbar. The relaxation parameter ω is given as the parameter omega. The data arrays for the volume forces in x-direction kVolx, in y-direction kVoly, and z-direction kVolz are expected in mN cm^− 3. The arrays for the velocity data v_x in Velx, v_y in Vely, and v_z in Velz are expected in mm s^− 1. The pressure data Press is expected in mbar. The arrays for the boundary conditions for velocity BCvel and pressure BCpress are expected as byte arrays. The string log is used to pass information back to the user. The maximum number of iterations maxiterations is a user-supplied value, which indicates the maximum number of iterations.

Listing 34.2

[C] Body of the three-dimensional solver implementation. This listing is part of the NeptunLib and can be found in the file NeptunLib.h in the supplementary material.

 1 int DLLEXPORTIMPORT SolveFlow3D(int L, int W, int H, double hxyz, double density, double viscosity, double epsilonVel, double epsilonPress, double omega, double* kVolx, double* kVoly, double* kVolz, double* Velx, double* Vely, double* Velz, double* Press, char* BCvel, char* BCpress, char* log, int maxiterations){
 2   /* Block 1: initialize variables and set constants */
 3   do{
 4    index = 0;
 5    nonconverged = 0;
 6    isconverged = 1;
 7    for(int x = 0;x < L;x++){
 8     for(int y = 0;y < W;y++){
 9      for(int z = 0;z < H;z++){
10       if(BCvel[index]! = 0 && BCpress[index]! = 0) {
11        index++;
12        continue;
13       }
14       islocalconverged = 1;
15       /* Block 2: update the structures */
16       /* Block 3: update the velocities */
17       /* Block 4: update the pressure */
18       UpdateFlowVariableDoNeumann();
19       if(isconverged && !islocalconverged) isconverged = 0;
20       if(!islocalconverged)nonconverged++;
21       index++;
22      }
23     }
24    }
25    count++;
26    if (count==maxiterations){
27     sprintf(log, "Reached maximum number of iterations = %i, converged cells = %.2f percent",count,100.0*(totalcells-nonconverged)/totalcells);
28     return 0;
29    }
30   }while(!isconverged);
31  sprintf(log, "Converged! Number of iterations: %i, converged cells = %.2f percent",count,100.0*(totalcells-nonconverged)/totalcells);
32  return 1;
33 }

As you can see, we cut out several blocks of the listing, which we will discuss in the following. At this point, it is enough to quickly summarize what these blocks do. Block 1 sets up all the variables and calculates the prefactors we require during calculation. Please refer to section 3.5.2 and listing 34.4 for details. We then step into the iteration loop which will proceed until the convergence flag isconverged is true (see line 30) or the maximum number of iterations is reached (see line 26). In both cases, we inform the user with a suitable message indicating the total number of iterations count and the percentage of cells that have converged. This value is calculated using the total number of cells totalcells and the total number of cells that are not converged nonconverged.

Upon entering a new iteration, we first reset the index variable index and the number of cells that have not converged nonconverged (see line 5). We also reset the convergence flag isconverged (see line 14). We then iterate overall all cells in the domain. For each cell we check if it is a boundary value for both velocity and pressure. If this is the case, we simply skip the cell (see line 10). In case either the velocity or the pressure is not a boundary value, we process the cell. We start by resetting the local convergence flag islocalconverged, which indicates if the current cell has converged (see line 14). In the next step, we update the structures velx, vely, velz, press, kvolx, kvoly, and kvolz, all of which are required to calculate the updated values v_x^(x,y,z), v_y^(x,y,z), v_z^(x,y,z), and p^(x,y,z). Updating these structures is performed by the code in block 2. Please refer to section 3.5.3 and listing 34.5 for details. Block 3 and block 4 update the velocities and the pressure, respectively. They are detailed in section 3.5.4 and listing 34.10.

Having updated the values v_x^(x,y,z), v_y^(x,y,z), v_z^(x,y,z), and p^(x,y,z) we update the values of all neighboring cells with Neumann boundary conditions (see line 18). This macro is detailed in section 3.5.6. Once this update has been done we check if the current cell has converged (see line 19). If this is not the case, the convergence check for this iteration isconverged is reset and the solver is required to iterate further. In this case, the variable nonconverged is incremented.

At the end of each internal loop we increment the index variable index to setup the link into the arrays correctly for the next cell. As you can see, taken by itself, the implementation of the solver is rather simple and compact. However, we have moved the most important fractions of the code into several blocks which we will detail in the following.

Return Values. As you can see, the function SolveFlow3D returns an integer value depending on the condition that terminated the iteration. A return value of 0 indicates that the solver has reached the maximum number of iterations without converging (see line 28). A return value of 1 indicates the solution has converged (see line 32). The last potential return value is − 1, which indicates abortion due to numerical instability. This condition is checked by block 3 (see section 3.5.4) and block 4 (see section 3.5.5).

Usage of Macros. As you will see, the implementation of the solver makes heavy use of macros. Obviously, in terms of readability and portability it would be better to use functions instead of macros. However, for our solver, there are two important advantages when using macros.

First, as macros are expanded and the code simply “copied” into the listing before compilation they are faster to execute. Functions require arguments and return pointer to be put on the call stack, which takes a couple of processor cycles. In most applications, these differences are hardly noticeable. However, our numerical solver will have to perform thousands of iterations on potentially millions of cells. Each clock cycle saved will make a noticeable difference in the end. Even something small such as the implementation of an absolute value function will save clock cycles if not written as a function but as a macro. Listing 34.3 shows this implementation that we will use instead of the call to the standard function abs which we have been using so far.

Listing 34.3

[C] The macro ABS of the three-dimensional solver implementation that returns the absolute of a floating-point value. This listing is part of the NeptunLib and can be found in the file NeptunLib.h in the supplementary material.

l #define ABS(value) (((value) < 0) ? –(value) : (value))

Second, as macros are essentially “text replacements” using a macro to write a function is less error-prone. It is very easy to mix up a direction variable or a step variable when implementing more complex numerical schemes. Using macros allows us to make the overall listing shorter and avoid typing errors because the code generation is done by the macro. Therefore we do not end up with, e.g., a velx.y.p where we should have written vely.x.p.

34.3.5.2 Block 1 - Initializing Variables and Constants

The first block of code taken from listing 34.2 is block 1. The listing is shown in listing 34.4. This block initializes variables and sets constants. The first variable we require is the iteration counter count (see line 1). If a given maximum number of iterations is reached, the code execution is aborted and the user is informed that the solver did not converge. The second set of variables we require are the index variables, which help us with accessing the data in the array. The variable index holds the current absolute index, whereas correctedindexvel and correctedindexpress hold corrected indices (see line 2). These are required whenever we are at the edge of the domain and there is no neighboring cell for which we could query data. These variables are used in the function UpdateFlowVariableCorrectlndex shown in listing 34.6. We also require the two convergence flags isconverged, which indicates if the current iteration is converged, i.e., if no cells have been found that have not converged, and islocalconverged, which indicates if the current cell has converged (see line 30). As we will see in a moment, we have to check a total of four values for convergence: v_x^(x,y,z), v_y^(x,y,z), v_z^(x,y,z), and p^(x,y,z). We will have a converged solution only if the changes between the last and the current iterations of all four variables are smaller than a given threshold. In order to temporarily store the new values for v_x^(x,y,z), v_y^(x,y,z), v_z^(x,y,z), and p^(x,y,z) we use the variables storex, storey, storez, and storep (see line 4).

Listing 34.4

[C] Block 1 of the three-dimensional solver implementation that sets up variables and constants. This listing is part of the NeptunLib and can be found in the file NeptunLib.h in the supplementary material.

 1 int count = 0;
 2 int index = 0;int correctedindexvel = 0; int correctedindexpress = 0;
 3 char isconverged,islocalconverged;
 4 double storep = 0; double storex = 0; double storey = 0; double storez = 0;
 5 currentflow3Dvariable velx,vely,velz;
 6 currentflow3Dvariable kvolx,kvoly,kvolz;
 7 currentflow3Dvariable press;
 8 velx.data = Velx;
 9 vely.data = Vely;
10 velz.data = Velz;
11 press.data = Press;
12 kvolx.data = kVolx;
13 kvoly.data = kVoly;
14 kvolz.data = kVolz;
15 double omega1 = 1-omega;
16 double omega2 = omega;
17 double denominatorconstant = 1000.0*6*viscosity/(density*hxyz);
18 double pressureconstantVel = 100000.0/density;
19 double laplaceconstantVel = 1000.0*viscosity/(density*hxyz);
20 double volumeforceconstantVel = hxyz/density;
21 double volumeforceconstantPress = hxyz/(100000.0*6.0);
22 double convectionconstantPress = density/(100000.0*6.0);
23 double pressureconstantPress = 1.0/6.0;
24 int xmax = L-1; int ymax = W-1; int zmax = H-1;
25 int xstep = W*H;int ystep = H; int zstep = 1;
26 long long nonconverged;
27 long long totalcells = L*W*H;
28 char NeumannVel_x_p = 0;char NeumannVel_x_m = 0;char NeumannVel_y_p = 0;char NeumannVel_y_m = 0;char NeumannVel_z_p = 0;char NeumannVel_z_m = 0;
29 char NeumannPress_x_p = 0;char NeumannPress_x_m = 0;char NeumannPress_y_p = 0;char NeumannPress_y_m = 0; char NeumannPress_z_p = 0;char NeumannPress_z_m = 0;

Listing 34.5

[C] Block 2 of the three-dimensional solver implementation that updates the current cell values and the values of the neighboring cells as well as the differentials. This listing is part of the NeptunLib and can be found in the file NeptunLib.h in the supplementary material.

 1 velx.value = velx.data[index];
 2 vely.value = vely.data[index];
 3 velz.value = velz.data[index];
 4 press.value = press.data[index];
 5 kvolx.value = kvolx.data[index];
 6 kvoly.value = kvoly.data[index];
 7 kvolz.value = kvolz.data[index];
 8 //x-direction
 9 UpdateFlowVariableCorrectIndex(x,xmax,xstep,NeumannVel_x_p,NeumannPress_x_p);
10 UpdateFlowVariableValuesAll(x,p);
11 UpdateFlowVariableCorrectIndex(x,0,-xstep,NeumannVel_x_m,NeumannPress_x_m);
12 UpdateFlowVariableValuesAll(x,m);
13 UpdateFlowVariableValuesDifferenceAll(x);
14 //y-direction
15 UpdateFlowVariableCorrectIndex(y,ymax,ystep,NeumannVel_y_p,NeumannPress_y_p);
16 UpdateFlowVariableValuesAll(y,p);
17 UpdateFlowVariableCorrectIndex(y,0,-ystep,NeumannVel_y_m,NeumannPress_y_m);
18 UpdateFlowVariableValuesAll(y,m);
19 UpdateFlowVariableValuesDifferenceAll(y);
20 //z-direction
21 UpdateFlowVariableCorrectIndex(z,zmax,zstep,NeumannVel_z_p,NeumannPress_z_p);
22 UpdateFlowVariableValuesAll(z,p);
23 UpdateFlowVariableCorrectIndex(z,0,-zstep,NeumannVel_z_m,NeumannPress_z_m);
24 UpdateFlowVariableValuesAll(z,m);
25 UpdateFlowVariableValuesDifferenceAll(z);

Listing 34.6

[C] Macro UpdateFlowVariableCorrectIndex of the three-dimensional solver implementation that corrects the index. This listing is part of the NeptunLib and can be found in the file NeptunLib.h in the supplementary material.

 1 #defineUpdateFlowVariableCorrectIndex(direction,boundary,step,NeumannVel,NeumannPress)
 2  if(direction==boundary){ 
 3   correctedindexvel = index; 
 4   correctedindexpress = index; 
 5   NeumannVel = 0; 
 6   NeumannPress = 0; 
 7  }else{
 8   correctedindexvel = index + (step); 
 9   if(BCvel[correctedindexvel]== − 1){ 
10    correctedindexvel = index; 
11    NeumannVel = − 1; 
12   }else{ 
13    NeumannVel = 0; 
14   } 
15   correctedindexpress = index + (step); 
16   if(BCpress[correctedindexpress]== − 1){ 
17    correctedindexpress = index; 
18    NeumannPress = − 1; 
19   }else{ 
20    NeumannPress = 0; 
21   } 
22  }

We then proceed to defining our calculation variables velx, vely, velz, and press (see line 11), all of which use the structure currentflow3Dvariable introduced in listing 34.1. We also use this structure to hold the volume forces kvolx, kvoly, and kvolz (see line 12). We do not have to calculate the volume forces as they are supplied by the user. However, we require the gradients of the volume forces for updating the velocity and the pressure values. Using the structure currentflow3Dvariable makes accessing these values significantly easier. We link up the values for all variables by assigning the structure element data to the arrays given by the user.

In the next step, we calculate the constants that we have already introduced in section 3.3 (see line 15). As stated, precalculating these factors in a way that they account for the correct units will make the overall calculation significantly faster. We then define the maximum values along each axis given by xmax, ymax, and zmax as well as the steps along each axis given by xstep, ystep, and zstep (see line 24). We will require all of these values during the calculation, therefore precalculating them saves computational effort.

As stated, we will keep track of the number of converged cells in the current iteration. This allows us to inform the user which degree of precision our iteration has obtained. We therefore introduce the variables nonconverged, which counts the number of cells that have not converged during the current iteration, and the variable totalcells, which stores the total number of cells (see line 26).

The last set of variables we require account for Neumann boundary conditions in neighboring cells (see line 28). As already stated, Neumann boundary values are updated whenever the value of the neighboring cell has been calculated. We will see this implementation when discussing the macro UpdateFlowVariableDoNeumann (see listing 34.11). Essentially, whenever the value of the cell (x, y, z) is updated, we test if the neighboring cells in positive and negative directions along each axis has a Neumann boundary value. If this is the case, we set the respective variable to 1. As an example, if we find that the boundary values for the velocity BCvel indicate that the cell in positive x-direction is a Neumann boundary value we set NeumannVel_x_p to 1. If the cell in negative direction is a Neumann boundary value, we set the flag NeumannVel_x_m to 1 etc. Once we determined the new values v_x^(x,y,z), v_y^(x,y,z), v_z^(x,y,z), for (x, y, z) we can check the variables NeumannVel_x_p, NeumannVel_x_m, NeumannVel_y_p, NeumannVel_y_m, NeumannVel_z_p, and NeumannVel_z_m to see if we should copy the freshly calculated values into the neighboring cells in order to account for Neumann boundary conditions.

34.3.5.3 Block 2 - Update the Structures

The second block of code taken from listing 34.2 is block 2. The listing is shown in listing 34.5. In this block, we update the structures velx, vely, velz, press, kvolx, kvoly, and kvolz. For all of these structures we first set the value in the current cell (x, y, z) (see line 1). This value is simply given by accessing the element in the array stores in the field data of the respective structure. The value is then stored in the field value of the respective structure. We then proceed by updating the values along the positive and negative x-direction. For this we first call the macro UpdateFlowVariableCorrectlndex, which is given in listing 34.6. Before continuing with the implementation of block 2 we will briefly look at the implementation of this macro.

The Macro UpdateFlowVariableCorrectIndex. As you can see from listing 34.5, line 9 the macro is called by passing the relevant direction variable, the boundary value of this index variable, the stepwidth as well as the relevant flags for Neumann boundary values for velocity and pressure. As we intend to correct the index for positive x-values, we pass x as the direction variable and xmax as the boundary value. The stepwidth is xstep and the relevant Neumann boundary flags are NeumannVel_x_p for the velocity and NeumannPress_x_p for the pressure.

The aim of this macro is to determine whether or not we can actually query a neighboring value to the current cell in the given direction. There are two cases we have to account for. The first case is that the neighboring cell in the given direction does not exist because we are already at the boundary of the domain. This case will not happen unless the user fails to declare the boundary values correctly. In the simplest case, all cells that are on the boundary should be defined as being boundary values. However, if the user fails to do so, we need to make sure that the solver does not access data that is not actually there. Therefore this first case must be accounted for. Whenever we come across this scenario, we simply assume a Neumann boundary condition.

The second case is that the neighboring cell in the given direction is a boundary value. If the neighboring cell is a Dirichlet boundary, the value can be accessed. If the cell is a Neumann boundary the value of the cell should not be directly accessed. In this case, we correct the index to access the value of the current cell instead. In addition we then need to set the correct flag to indicate that once the current cell’s value has been updated, we should update the value of the boundary condition too.

The macro first checks for the first case, i.e., if the direction variable is already at the boundary (see line 2). If this is the case, the corrected index variables correctedindexvel and correctedindexpress are set to index, which is the index to access the value of the current cell. In this case, instead of accessing the value of the neighboring cell (x + 1, y, z), which does not exist, we access the value of the current cell (x, y, z) thereby emulating a Neumann boundary condition. Since there are no neighboring cells, we cannot update the cell’s, value. Therefore we reset the flags for Neumann boundary conditions on the velocity and the pressure.

The second case is the more common and it is accounted for by the else block starting in line 7. As we are not on the boundary, we can correct the index variables correctedindexvel and correctedindexpress by simply adding step to the current index variable index (see line 8). We then check the velocity for boundary values by inspecting the value in BCvel. A value of − 1 indicated a Neumann boundary value (see line 9). If this value is found, we reset correctedindexvel to index and set the flag for the Neumann boundary values to 1. This effectively means that instead of accessing the value of neighboring cell (x + 1, y, z) we access the value of the current cell (x, y, z). Once this value has been updated, we copy the value into the neighboring cell (see line 10). If we find a value different from − 1 in the boundary conditions we leave the index corrected and reset the Neumann boundary value flag to 0 (see line 13).

After checking the velocity, we perform the same procedure on the pressure field correcting the index variable correctedindexpress and setting the Neumann boundary value flag NeumannPress.

After calling the macro UpdateFlowVariableCorrectIndex we can safely access the value of the neighboring cell using correctedindexvel for the velocity and correctedindexpress for the pressure fields. Even if the current cell happens to lie on the edge of the domain without being a boundary we have accounted for that. Furthermore, if the neighboring cell happened to be a Neumann boundary we have correctly set the flags NeumannVel for the velocity and NeumannPress for the pressure field and “deviated” the access to the cell’s value to the value of the current cell instead.

The Macros UpdateFlowVariableValuesAll and UpdateFlowVariableValuesSingle. Before returning to the implementation of block 2 in listing 34.5 we will briefly introduce the macros UpdateFlowVariableValuesAll and UpdateFlowVariableValuesSingle (see listing 34.7). Having updated the values correctedindexvel and correctedindexpress of the neighboring values we can now update the values along the given axis direction in all of our structure variables. This is done by calling the macro UpdateFlowVariableValuesAll with argument direction (either x, y, or z) and the increment (either p for “plus” or m for “minus”). The macro simply calls the macro UpdateFlowVariableValuesSingle, which expects the variable to set variable, the direction direction (again either x, y, or z), the increment (again p or m), and the corrected index correctedindex, which should be correctedindexvel when accessing the velocity field and correctedindexpress when accessing the pressure and force fields.

Listing 34.7

[C] Macros UpdateFlowVariableValuesAll and UpdateFlowVariableValuesSingle of the three dimensional solver implementation that update the values of the structures. This listing is part of the NeptunLib and can be found in the file NeptunLib.h in the supplementary material.

 1 #define UpdateFlowVariableValuesSingle(variable,direction,increment,correctedindex) 
 2  variable.direction.increment = variable.data[correctedindex]
 3
 4 #define UpdateFlowVariableValuesAll(direction,increment) 
 5  UpdateFlowVariableValuesSingle(velx,direction,increment,correctedindexvel); 
 6  UpdateFlowVariableValuesSingle(vely,direction,increment,correctedindexvel); 
 7  UpdateFlowVariableValuesSingle(velz,direction,increment,correctedindexvel); 
 8  UpdateFlowVariableValuesSingle(press,direction,increment,correctedindexpress); 
 9  UpdateFlowVariableValuesSingle(kvolx,direction,increment,correctedindexpress); 
10  UpdateFlowVariableValuesSingle(kvoly,direction,increment,correctedindexpress); 
11  UpdateFlowVariableValuesSingle(kvolz,direction,increment,correctedindexpress)

When calling UpdateFlowVariableValuesAll with the arguments x and p the macro will update all values {velx, vely, velz, press, kvolx, kvoly, kvolz}.x.p using the corrected indices correctedindexvel and correctedindexpress. Please note that these values have just been updated because we have called UpdateFlowVariableCorrectIndex right before calling UpdateFlowVariableValuesAll in listing 34.5, line 9.

Returning to Block 2. Having corrected the indices for the velocity and pressure fields in positive direction in listing 34.5, line 9 and setting the correct values in our structures velx, vely, velz, press, kvolx, kvoly, and kvolz by calling UpdateFlowVariableValuesAll for the positive x-direction we then proceed by doing the same thing for the negative x-direction (see listing 34.5, line 11). Obviously, we have to pass the lower boundary 0 and the negative step as well as the Neumann boundary flags NeumannVel_x_m and NeumannPress_x_m. Having updated the values for the positive and negative changes in x-direction we proceed by updating the difference, i.e., the differential along the x-direction by calling the macro UpdateFlowVariableValuesDifferenceAll (see listing 34.5, line 13). We will discuss this macro that is shown in listing 34.8 in a moment.

Listing 34.8

[C] Macros UpdateFlowVariableValuesDifferenceAll and UpdateFlowVariableValuesDifference of the three-dimensional solver implementation that update the differential values of the structures. This listing is part of the NeptunLib and can be found in the file NeptunLib.h in the supplementary material.

 1 #define UpdateFlowVariableValuesDifference(variable,direction) 
 2  variable.direction.d = 0.5*(variable.direction.p-variable.direction.m)
 3
 4 #define UpdateFlowVariableValuesDifferenceAll(direction) 
 5  UpdateFlowVariableValuesDifference(velx,direction); 
 6  UpdateFlowVariableValuesDifference(vely,direction); 
 7  UpdateFlowVariableValuesDifference(velz,direction); 
 8  UpdateFlowVariableValuesDifference(press,direction); 
 9  UpdateFlowVariableValuesDifference(kvolx,direction); 
10  UpdateFlowVariableValuesDifference(kvoly,direction); 
11  UpdateFlowVariableValuesDifference(kvolz,direction)

Having updated the values for both positive and negative changes as well as the differential along the x-direction, we proceed likewise for the y-direction (see listing 34.5, line 15) and the z-direction (see listing 34.5, line 21).

At the end of block 2 we have updated all fields in the our structures velx, vely, velz, press, kvolx, kvoly, and kvolz. With these updated values, we can now proceed to calculating the updated values for the current cells $v_{x}^{(x, y, z)}, v_{y}^{(x, y, z)}, v_{z}^{(x, y, z)},$ $v_{x}^{(x, y, z)}, v_{y}^{(x, y, z)}, v_{z}^{(x, y, z)},$ and p^{(x, y, z)}.

The MacrosUpdateFlowVariableValuesDifferenceAllandUpdateFlowVariableValuesDifference. As stated, we are still missing two macros for the implementation of block 2. The first macro we require is UpdateFlowVariableValuesDifferenceAll (see listing 34.8). It is called with a single parameter direction that is either x, y, or z. It then calls the macro UpdateFlowVariableValuesDifference, which updates the field .d and a given direction direction (either x, y, or z) on a variable. As you can see line 2 is the central difference scheme (see Eq. 4.6), which simply uses the value of the neighboring cell in positive and negative directions, which we just updated.

34.3.5.4 Block 3 - Update the Velocities

The third block of code taken from listing 34.2 is block 3. In this block, we update the velocity values $v_{x}^{(x, y, z)}, v_{y}^{(x, y, z)}$ $v_{x}^{(x, y, z)}, v_{y}^{(x, y, z)}$ and $v_{z}^{(x, y, z)}$ $v_{z}^{(x, y, z)}$ using the equations we derived in section 2.3.2. The listing is shown in listing 34.9. We first check if the given cell is a boundary value by inspecting BCvel. If a value different from 0 is found, we assume a boundary value and skip the updating of the velocity.

Listing 34.9

[C] Block 3 of the three-dimensional solver implementation that updates the current cell’s velocities values. This listing is part of the NeptunLib and can be found in the file NeptunLib.h in the supplementary material.

 1 if(BCvel[index]==0){
 2  storex = omega1*velx.value + omega2*(
 3    − (vely.value*velx.y.d + velz.value*velx.z.d)
 4    − pressureconstantVel*press.x.d
 5    + laplaceconstantVel*(velx.x.p + velx.x.m + velx.y.p + velx.y.m + velx.z.p + velx.z.m)
 6    + volumeforceconstantVel*kvolx.value
 7   )/(velx.x.d + denominatorconstant);
 8  storey = omega1*vely.value + omega2*(
 9    − (velx.value*vely.x.d + velz.value*vely.z.d)
10    − pressureconstantVel*press.y.d
11    + laplaceconstantVel*(vely.x.p + vely.x.m + vely.y.p + vely.y.m + vely.z.p + vely.z.m)
12    + volumeforceconstantVel*kvoly.value
13   )/(vely.y.d + denominatorconstant);
14  storez = omega1*velz.value + omega2*(
15    − (velx.value*velz.x.d + vely.value*velz.y.d)
16    − pressureconstantVel*press.z.d
17    + laplaceconstantVel*(velz.x.p + velz.x.m + velz.y.p + velz.y.m + velz.z.p + velz.z.m)
18    + volumeforceconstantVel*kvolz.value
19   )/(velz.z.d + denominatorconstant);
20  if(islocalconverged && (ABS(storex-velx.value) > epsilonVel || ABS(storey-vely.value) > epsilonVel || ABS(storez-velz.value) > epsilonVel)) islocalconverged = 0;
21  if(ABS(storex) > STABILITYTHRESHOLD ||ABS(storey) > STABILITYTHRESHOLD || ABS(storez) > STABILITYTHRESHOLD){
22   sprintf(log, “Numerical instability detected! Aborting…”);
23   return − 1;
24  }
25  Velx[index] = storex;
26  Vely[index] = storey;
27  Velz[index] = storez;
28 }

Listing 34.10

[C] Block 4 of the three-dimensional solver implementation that updates the current cell’s pressure value. This listing is part of the NeptunLib and can be found in the file NeptunLib.h in the supplementary material.

1 if(BCpress[index]==0){
2  storep = omega1*press.value + omega2*(
3    − volumeforceconstantPress*(kvolx.x.d + kvoly.y.d + kvolz.z.d)
4    + convectionconstantPress*(velx.x.d*velx.x.d + vely.y.d*vely.y.d + velz.z.d*velz.y.d + 2*(vely.x.d* velx.y.d + velx.z.d*velz.x.d + vely.z.d*velz.y.d))
5    + pressureconstantPress*(press.x.p + press.x.m + press.y.p + press.y.m + press.z.p + press.z.m)
6  );
7  if(islocalconverged && ABS(storep-press.value) > epsilonPress) islocalconverged = 0;
8  Press[index] = storep;
9 }

If no boundary value was found, the velocity values must be updated. We begin with $v_{x}^{(x, y, z)}$ $v_{x}^{(x, y, z)}$ using Eq. 34.34 to calculate storex by means of an SOR scheme (see line 2). As you can see, we make use of the precalculated factors at this point, which slims down the implementation significantly. The updated value for $v_{x}^{(x, y, z)}$ $v_{x}^{(x, y, z)}$ which is temporarily stored in storex is not directly assigned because we want to check convergence and numerical stability before assigning the value. We proceed likewise with $v_{y}^{(x, y, z)}$ $v_{y}^{(x, y, z)}$ which implements Eq. 34.35 for updating storey (see line 8) and with $v_{z}^{(x, y, z)}$ $v_{z}^{(x, y, z)}$ which implements Eq. 34.36 for updating storez (see line 14).

Having assigned storex, storey, and storez we check convergence in line 20 by comparing the absolute difference between the respective values of the last and the current iterations. If any of these differences is larger than epsilonVel we reset the local convergence flag islocalconverged. In the next step, we check for numerical stability. This can be done conveniently by checking if any of the updated values, storex, storey, or storez, are larger than the given stability threshold value STABILITYTHRESHOLD. This constant is currently set to 10 × 10⁸. We can safely assume that the velocity in our domain should never achieve such high velocity values. If such high values are encountered, we can assume its due to the fact that the numerical scheme has become instable. If any of the values exceed the threshold, we inform the user by setting log and exit the function with the exit code − 1, which indicates abortion due to numerical instability.

In case the scheme was not found to be numerically instable, we copy the updated values for $v_{x}^{(x, y, z)}, v_{y}^{(x, y, z)}$ $v_{x}^{(x, y, z)}, v_{y}^{(x, y, z)}$ , and $v_{z}^{(x, y, z)}$ $v_{z}^{(x, y, z)}$ from storex, storey, and storez into the arrays Velx, Vely, and Velz at the position given by the index index.

34.3.5.5 Block 4 - Update the Pressure

The fourth block of code taken from listing 34.2 is block 4. In this block, we update the pressure value p^{(x, y, z)}.In analogy to the velocity field we first check whether or not the current cell is a boundary value for the pressure.

If this is the case, the cell is skipped, otherwise the updated pressure value is calculated using Eq. 34.38. The value is assigned to storep and checked for convergence. Please note that we support different threshold values for velocity and pressure. In case of the pressure we compare against the value epsilonPress. The value is then copied into the array Press.

34.3.5.6 The macrosUpdateFlowVariableDoNeumannandUpdateFlowVariableDoNeumannSingle

The last two macros that are missing for our implementations are UpdateFlowVariableDoNeumann and UpdateFlowVariableDoNeumannSingle, shown in listing 34.11. These macros account for potential Neumann boundary conditions in the vicinity of the current cell. As you can see, UpdateFlowVariableDoNeumann calls UpdateFlowVariableDoNeumannSingle using the respective Neumann boundary value flags of the velocities for each axis NeumannVel_x, NeumannVel_y, and NeumannVel_z for positive and negative directions _p and _m, respectively. It passes on the stepsize along the respective direction as well as the value to copy (storex, storey, storez, and storep).

Listing 34.11

[C] The macros UpdateFlowVariableDoNeumann and UpdateFlowVariableDoNeumannSingle of the three dimensional solver implementation that updates the neighboring cells that have Neumann boundary conditions. This listing is part of the NeptunLib and can be found in the file NeptunLib.h in the supplementary material.

 1 #define UpdateFlowVariableDoNeumannSingle(Neumann,data,step,value)
 2  if((Neumann)== − 1) data[index + (step)] = (value); 
 3
 4 #define UpdateFlowVariableDoNeumann( )
 5  UpdateFlowVariableDoNeumannSingle(NeumannVel_x_p,Velx,xstep,storex); 
 6  UpdateFlowVariableDoNeumannSingle(NeumannVel_x_m,Velx,-xstep,storex);
 7  UpdateFlowVariableDoNeumannSingle(NeumannVel_x_p,Vely,xstep,storey); 
 8  UpdateFlowVariableDoNeumannSingle(NeumannVel_x_m,Vely,-xstep,storey);
 9  UpdateFlowVariableDoNeumannSingle(NeumannVel_x_p,Velz,xstep,storez); 
10  UpdateFlowVariableDoNeumannSingle(NeumannVel_x_m,Velz,-xstep,storez); 
11  
12  UpdateFlowVariableDoNeumannSingle(NeumannVel_y_p,Velx,ystep,storex);
13  UpdateFlowVariableDoNeumannSingle(NeumannVel_y_m,Velx,-ystep,storex); 
14  UpdateFlowVariableDoNeumannSingle(NeumannVel_y_p,Vely,ystep,storey);
15  UpdateFlowVariableDoNeumannSingle(NeumannVel_y_m,Vely,-ystep,storey); 
16  UpdateFlowVariableDoNeumannSingle(NeumannVel_y_p,Velz,ystep,storez);
17  UpdateFlowVariableDoNeumannSingle(NeumannVel_y_m,Velz,-ystep,storez); 
18  
19  UpdateFlowVariableDoNeumannSingle(NeumannVel_z_p,Velx,zstep,storex);
20  UpdateFlowVariableDoNeumannSingle(NeumannVel_z_m,Velx,-zstep,storex); 
21  UpdateFlowVariableDoNeumannSingle(NeumannVel_z_p,Vely,zstep,storey);
22  UpdateFlowVariableDoNeumannSingle(NeumannVel_z_m,Vely,-zstep,storey); 
23  UpdateFlowVariableDoNeumannSingle(NeumannVel_z_p,Velz,zstep,storez);
24  UpdateFlowVariableDoNeumannSingle(NeumannVel_z_m,Velz,zstep,storez); 
25  
26  UpdateFlowVariableDoNeumannSingle(NeumannPress_x_p,Press,xstep,storep)
27  UpdateFlowVariableDoNeumannSingle(NeumannPress_x_m,Press,xstep,storep
28  UpdateFlowVariableDoNeumannSingle(NeumannPress_y_p,Press,ystep,storep)
29  UpdateFlowVariableDoNeumannSingle(NeumannPress_y_m,Press,-ystep,storep);
30  UpdateFlowVariableDoNeumannSingle(NeumannPress_z_p,Press,zstep,storep)
31  UpdateFlowVariableDoNeumannSingle(NeumannPress_z_m,Press,-zstep,storep)

UpdateFlowVariableDoNeumannSingle essentially checks if the given Neumann boundary value flag is − 1, in which case it simply copies value into the given array data at the position index + step. After calling UpdateFlowVariableDoNeumann the values of all neighboring cells of the current cell (x, y, z) that have been identified as Neumann boundary values will have the same value of the current cell.

34.4 Usage of the Numerical Solver

34.4.1 Introduction

We have now finished the implementation of the three-dimensional numerical solver and are ready to test the solver using various flow cases. We will study a total of 6 cases. For each case, we will write a custom function to set up the boundary and the initial values in Maple whereupon we call the function SolveFlow3D for solving the flow problem.

Before we start, we first introduce the first few lines of the Maple worksheet that call the DLL. They are shown in listing 34.12. We first restart the Maple server and include Core.txt and CoreFunctions.txt. We then define the Maple function C_SolveFlow3D, which calls the function SolveFlow3D written in C via the NeptunLib.dll (see line 3). The parameters of the function have already been introduced. All data arrays are passed by reference as arrays of size L · W · H with 8-byte floating point values. The boundary conditions BCvel and BCpress are not passed as boolean arrays but rather as byte arrays since we require the values − 1, 0, and 1. We pass the string log with a maximum length of 256 characters and expect the function to return an integer value.

Listing 34.12

[C] Maple worksheet for solving the flow cases using the three-dimensional solver. These are the first few lines of the Maple worksheet that uses the function SolveFlow3D implemented in C and called via a DLL. This listing is part of the listing ThreedimensionalFlowsNumerical.mw that can be found in the supplementary material.

 1 restart:read “Core.txt”: read “CoreFunctions.txt”:
 2
 3  C_SolveFlow3D:=define_external(‘SolveFlow3D’,LIB = “NeptunLib.dll”,L::(integer[4]),W::(integer[4]),H::(integer[4]),hxyz::(float[8]),density::(float[8]),viscosity::(float[8]),epsilonVel::(float[8]),epsilonPress::(float[8]),omega::(float[8]),kVolx::ARRAY(1..L*W*H,float[8]),kVoly::ARRAY(1..L*W*H,float[8]),kVolz::ARRAY(1..L*W*H,float[8]),dataVelx’::ARRAY(1..L*W*H,float[8]),‘dataVely’::ARRAY(1..L*W*H,float[8]),‘dataVelz’::ARRAY(1..L*W*H,float[8]),‘dataPress’::ARRAY(1..L*W*H,float[8]),‘BCvel’::ARRAY(1.L*W*H,integer[1]), ‘BCpress’::ARRAY(1..L*W*H,integer[1]), ‘ log’ ::(string[256]),maxiterations::(integer[4]),RETURN::(integer[4]));
 4
 5 PrepareFlow3D:=proc(L,W,H) local dataVelx, dataVely, dataVelz, dataPress, kVolx, kVoly, kVolz, BCvel, BCpress:
 6  dataVelx:=Array(1..L,1..W,1..H,datatype = float[8],order = C_order):
 7  dataVely:=Array(1..L,1..W,1..H,datatype = float[8],order = C_order):
 8  dataVelz:=Array(1..L,1..W,1..H,datatype = float[8],order = C_order):
 9  dataPress:=Array(1..L,1..W,1..H,datatype = float[8],order = C_order):
10  kVolx:=Array(1..L,1..W,1..H,datatype = float[8],order = C_order):
11  kVoly:=Array(1..L,1..W,1..H,datatype = float[8],order = C_order):
12  kVolz:=Array(1..L,1..W,1..H,datatype = float[8],order = C_order):
13  BCvel:=Array(1..L,1..W,1..H,datatype = integer[1],order = C_order):
14  BCpress:=Array(1..L,1..W,1..H,datatype = integer[1],order = C_order):
15  return L,W,H,dataVelx,dataVely,dataVelz,dataPress,kVolx,kVoly,kVolz,BCvel,BCpress:
16 end proc:
17
18  SolveFlow3D:=proc(hxyz,density,viscosity,L,W,H,dataVelx,dataVely,dataVelz,dataPress,kVolx,kVoly, kVolz,BCvel,BCpress,maxiterations) local log:="":
19   C_SolveFlow3D(L,W,H,hxyz,density,viscosity,1E-6,1E 6,1.8,kVolx,kVoly,kVolz,dataVelx,dataVely, dataVelz,dataPress,BCvel,BCpress,'log',maxiterations):
20  print(log);
21  return [(x,y,z) → dataVelx[trunc(x),trunc(y),trunc(z)],(x,y,z) > dataVely[trunc(x),trunc(y),trunc(z )],(x,y,z) → dataVelz[trunc(x),trunc(y),trunc(z)],(x,y,z) → dataPress[trunc(x),trunc(y),trunc(z) ],dataVelx,dataVely,dataVelz,dataPress,L,W,H]:
22 end proc:

In the next step, we implement the function PrepareFlow3D (see line 5). This function creates the required arrays in the correct size, data format, and layout and returns an array containing the length L, width W, height H, and the pointers to the arrays dataVelx, dataVely, dataVelz, dataPress, kVolx, kVoly, kVolz, BCvel, and BCpress. For all of the following flow cases, we will first call PrepareFlow3D to prepare the arrays. In the next step, we will modify the data to account for the boundary values and the initial values, which are different for each of the cases we study.

The last function we define is SolveFlow3D (see line 18). This function invokes the C function SolveFlow3D after being provided with the required input values. The first three arguments this function expects are the step width h_xyz, the density ρ, and the viscosity η of the fluid. Please note that the following 12 parameters are in the same order as the return value of PrepareFlow3D. We can therefore simply plug in the return value of this function here. The last parameter that we need to modify is the maximum number of iterations maxiterations, which may be different for each flow case. The function then proceeds by calling the C function and printing the log message stored in log. The return value of the Maple function SolveFlow3D is an array of four functions which simply returns the values for v_x, v_y, v_z, and p as functions of x, y, and z. The use of the function trunc allows these functions to also be supplied with noninteger values. Normally, Maple would not allow using floating-point values as array indices, which makes it more difficult to plot the velocity and pressure fields as continuous functions. By using this little trick we simply round floating-point values to integer values, which allows us to use the regular Maple plot commands for displaying the computation results.

Having written these functions we are now ready to test our solver using the first flow case.

34.4.2 Case 1: Chamber By-Flow

The first case we will discuss is a flow that flows by a chamber. The geometry is shown in Fig. 34.1a. As you can see the channel is smaller than the chamber and the flow is expected to extend into the chamber to a certain degree. Please note that this flow case is close to impossible to tackle analytically. However, as we will in a moment, our custom-written solver finds a solution conveniently.

f34-01-9781455731411 — Fig. 34.1 Usage of the three-dimensional solver for case 1: channel by-flow. a) Geometry of the flow problem. A microfluidic channel of width 20 μm flows along a chamber of length L = 400 μm and width W = 100 μm. The height of both the chamber and the channel is H = 100 μm. For this case, we define a pressure drop of − 0.01 mbar across the domain. The stepwidth was h_xyz = 2 μm. The relaxation parameter was ω = 1.8 and the precision for both velocity and pressure was set to 1 × 10^− 6. The fluid was chosen to be water with ρ = 1gcm^− 3and η = 1mPas. b) Numerical output obtained with the implemented three-dimensional solver. The output is the velocity field in the xy-plane at z = 50 μm (see a) overlaid on top of the pressure field. As you can see, the pressure is highest at the entrance and lowest at the exit. The flow extends into the chamber creating a by-flow.

As stated, we need a custom function to set up the initial values and the boundary values for this case. This function is termed SetupBCIVCase1 and shown in listing 34.13. It accepts an inlet pressure p0 and the entry width entrywidth of the inflowing channel. The remaining arguments are the 12 values returned by PrepareFlow3D. If you look at line 32 you can see that this makes it convenient to call the setup function by simply calling PrepareFlow3D in the asrguments.

Listing 34.13

[C] Maple worksheet for solving case 1 using the three-dimensional solver. This listing is part of the listing ThreedimensionalFlowsNumerical.mw, which can be found in the supplementary material.

 1 SetupBCIVCase1:=proc(p0,entrywidth,L,W,H,dataVelx,dataVely,dataVelz,dataPress,kVolx,kVoly,kVolz, BCvel, BCpress) local x, y, z:
 2  for y from 1 to W do
 3   for z from 1 to H do
 4    if y < entrywidth then
 5     BCvel[1, y, z]:= − 1:
 6     BCvel[L, y, z]:= − 1:
 7     BCpress[1,y,z]:=1:
 8     BCpress[L, y, z]:=1:
 9     dataPress[1,y,z]:=p0:
10     dataPress[L, y, z]:=0:
11    else
12     BCvel[1,y,z]: = 1:
13     BCvel[L,y,z]: = 1:
14     BCpress[1,y,z]:= − 1:
15     BCpress[L,y,z]:= − 1:
16    end if:
17   end do:
18  end do:
19  for x from 1 to L do
20   for y from 1 to W doss
21    for z from 1 to H do
22     if y = 1 or y = W or z = 1 or z = H then
23      BCvel[x, y, z]:=1:
24      BCpress[x, y, z]:= − 1:
25     end if
26    end do:
27   end do:
28  end do:
29  return L,W,H,dataVelx,dataVely,dataVelz,dataPress,kVolx,kVoly,kVolz,BCvel,BCpress:
30 end proc:
31
32 dataCase1:=SetupBCIVCase1(0.01,10,PrepareFlow3D(200,50,50)):
33 sol1:=SolveFlow3D(2,1,1,dataCase1,20000):
34
35 fieldplot([(x,y) → sol1[1](x,y,25),(x,y) → sol1[2](x,y,25)],1..200,1..50);
36 densityplot((x,y) → sol1[4](x,y,25),1..200,1..50,colorscheme = ["Blue","Green","Yellow","Red"]);
37  fieldplot3d([(x,y,z) → sol1[1](x,y,z),(x,y,z) → sol1[2](x,y,z),(x,y,z) → sol1[3](x,y,z) ],1..200,0..50,0..50,grid = [8,8,8],arrows = THICK);

At the entrance and the exit of the domain, i.e., for x = 1 and x = L, the function then sets boundary values depending on the value of y. For all values y <entrywidth the boundary values are Neumann type for velocity and Dirichlet type for the pressure (see line 5). At the entrance (x = 1) the pressure is set to the given value p0 whereas at the exit (x = L) the pressure is set to 0. For values y ≥ entrywidth the boundary values are Dirichlet type for the velocity and Neumann type for the pressure (see line 12). This corresponds to a no-slip boundary condition.

Having set the values for the entrance and the exit of the domain, we then proceed setting the no-slip boundary condition for all values of (x, y, z) wherever y = 1 or y = W or z = 1 or z = H (see line 22). Again, no-slip boundary conditions correspond to a Dirichlet boundary condition for the velocity and a Neumann boundary condition for the pressure. The function then simply returns the 12 variables it obtained from PrepareFlow3D and modified to suit the boundary conditions and the initial values.

Pressure Drop and Absolute Pressure. As you can see from line 32 we use an entry pressure of 0.01 mbar in our example. The outlet pressure is set to be 0 mbar, which amounts to a pressure drop of − 0.01 mbar along the x-direction in our domain. The exact values of the inlet and the outlet pressure can be offset by any given values. This is due to the fact that the Navier–Stokes (see Eq. 11.40) does not account for absolute pressure but only for gradients in the pressure, i.e., pressure changes. We can therefore give the inlet and the outlet pressure as boundary values. They must be regarded as being “offsets” to, e.g., ambient pressure. As long as the pressure drop is the same, the obtained velocity field will also be identical. The absolute pressure values do not matter in this regard.

A First Calculation. In line 32 we have set up a domain with 200 cells along the x-direction and 50 cells along the y-direction and the z-direction, respectively. The channel width was chosen to be 10 cells. In the next step, we call SolveFlow3D (see line 33) providing a stepwidth sh_xyz = 2 pm as well as the density and the viscosity of water. The step width will scale our domain to a length L = 400 pm, a width W = 100 pm, and a height H = 100 pm. The inflowing channel will then have a width of 20 pm. We request a maximum number of 20 000 iterations because we cannot be sure how many iterations we actually require. As we will see, we will require a total of 18 899 iterations. However, this calculation will take some time and it may be tedious to wait for the solver to complete. The way it is currently written, we have no feedback on the current number of iterations and the total percentage of converged cells. Therefore it makes sense to start the calculation with a smaller number of maximum iterations, e.g., with 1000. Depending on the system you are running the solver on, this should not take more than a minute. Once completed, the solver will output that it reached the maximum number of iterations but did not converge. It will also give a total percentage of converged cells, which will be the percentage of cells holding boundary values. These cells obviously do not have to converge.

Visualization. As stated, the solver will take 18 899 iterations to complete. Once it returns, it provides us with functions for plotting v_x (x, y, z), v_y (x, y, z), v_z (x, y, z), and p (x, y, z).Fig. 34.1b shows the velocity field at z = 25, i.e., z = 50 μm overlaid onto the pressure distribution. These plots are generated by line 35 and line 36. Fig. 34.2 shows a three-dimensional view of the velocity field created by line 37. As you can see, the flow extends into the chamber. In the three-dimensional view you can see that the flow velocity is smaller at the boundaries along the z-axis due to the no-slip boundary conditions. The pressure is highest at the inlet and lowest at the outlet. It is also not uniform across the chamber due to the way the velocity field expands in a nonuniform way. Overall this result agrees very well with what we would expect in such a scenario. This is already a very complex example of a three-dimensional flow, and, as stated, it is very simple to solve using our custom-implemented numerical solver.

f34-02-9781455731411 — Fig. 34.2 Three-dimensional flow field for case 1: channel by-flow. Fig. 34.1b shows a cut of the velocity field at constant z. This graph shows the entire three-dimensional flow field. You can see that the by-flow extends along the z-axis but becomes smaller as z approaches the no-slip boundaries at the top and the bottom of the channel.

34.4.3 Case 2: Chamber Through-Flow

The second case we will discuss is similar to the first; only this time, we arrange the channel such that it crosses the chamber in the middle. The geometry is shown in Fig. 34.3a. Again, we would expect to see an extending flow but this time symmetrically.

f34-03-9781455731411 — Fig. 34.3 Usage of the three-dimensional solver for case 2: channel through-flow. a) Geometry of the flow problem. A microfluidic channel of width w = 20 μm flows through a chamber of length L = 400 μm and width W = 100 μm. The height of both the chamber and the channel is H = 100 μm. For this case, we define a pressure drop of − 0.01 mbar across the domain. The stepwidth was h_xyz = 2 μm. The relaxation parameter was ρ = 1.8 and the precision for both velocity and pressure was set to 1 × 10^− 6. The fluid was chosen to be water with ρ = 1gcm^− 3and η = ImPas. b) Numerical output obtained with implemented three-dimensional solver. The output is the velocity field in the xy-plane at z = 50 μm (see a) overlaid on top of the pressure field. As you can see, the pressure is highest at the entrance where the microfluidic channel flows into the chamber. The lowest value is at the exit. The flow extends symmetrically into the chamber.

As you can see the channel is smaller than the chamber and the flow is expected to extend into the chamber to a certain degree. Please note that this flow case is close-to-impossible to tackle analytically. However, our custom-written solver finds a solution conveniently. The funsction for setting up this boundary condition is similar to SetupBCIVCase1 of listing 34.13 and shown in listing 34.14. It checks if the current value for y is sufficiently close the chamber center halfwidth. The rest of the function is identical to SetupBCIVCase1.

Listing 34.14

[Maple] Maple worksheet for solving case 2 using the three-dimensional solver. This listing is part of the listing ThreedimensionalFlowsNumerical.mw, which can be found in the supplementary material.

 1 SetupBCIVCase2:=proc(p0,entrywidth,L,W,H,dataVelx,dataVely,dataVelz,dataPress,kVolx,kVoly,kVolz, BCvel,BCpress) local x,y,z,halfwidth,halfchannelwidth:
 2  halfwidth:=W/2;
 3  halfchannelwidth:=entrywidth/2:
 4  for y from 1 to W do
 5  for z from 1 to H do
 6   if abs(halfwidth − y) < halfchannelwidth then
 7    BCvel[1,y,z]:= − 1:
 8    BCvel[L,y,z]:= − 1:
 9    BCpress[1,y,z]:=1:
10    BCpress[L,y,z]:=1:
11    dataPress[1,y,z]:=p0:
12   else
13    BCvel[1,y,z]:=1:
14    BCvel[L,y,z]:=1:
15    BCpress[1,y,z]:= − 1:
16    BCpress[L,y,z]:= − 1:
17   end if:
18  end do:
19 end do:
20 for x from 1 to L do
21  for y from 1 to W do
22   for z from 1 to H do
23    if y = 1 or y = W or z = 1 or z = H then
24     BCvel[x,y,z]:=1:
25     BCpress[x,y,z]:= − 1:
26    end if
27   end do:
28  end do:
29 end do:
30 return L,W,H,dataVelx,dataVely,dataVelz,dataPress,kVolx,kVoly,kVolz,BCvel,BCpress:
31 end proc:
32
33 dataCase2:=SetupBCIVCase2(0.01,10,PrepareFlow3D(200,50,50)):
34 sol2:=SolveFlow3D(2,1,1,dataCase2,20000):
35
36 fieldplot([(x,y) → sol2[1](x,y,25),(x,y) → sol2[2](x,y,25)],1..200,1..50);
37 densityplot((x,y) → sol2[4](x,y,25),1..200,1..50,colorscheme = ["Blue","Green","Yellow","Red"]);
38 fieldplot3d([(x,y,z) → sol2[1](x,y,z),(x,y,z) → sol2[2](x,y,z),(x,y,z) → sol2[3](x,y,z) ],1..200,0..50,0..50,grid = [8,8,8],arrows = THICK);

Visualization. Again we use a pressure drop of − 0.01 mbar across the x-direction of the domain and an entrywidth of 10 cells, which corresponds to a value of 20 μm. We call SolveFlow3D with the properties of water for the fluid and a maximum number of iterations of 20 000 (see line 34). The solver requires 15 481 iterations to complete. Again, we plot the velocity field in the xy-plane at z = 50 μm, which is shown with the overlay of the pressure field in Fig. 34.3b. As you can see, the flow expands into the chamber pretty rapidly requiring only very short amounts of time. If you recall, the speed at which momentum diffuses in water (see Fig. 9.6) is high: For diffusing 50 μm less than 10 ms are required, which is why we see effective momentum diffusion into the chamber in the output data. Again, if you consider the pressure field you can see that the highest pressure is at the inlet of the channel, whereas the lowest pressure is at the outlet. Fig. 34.4 shows the expansion of the flow at the entrance and the recompression of the flow at the outlet. This plot is created by line 38. If you rotate the three-dimensional view in Maple you can see that the flow velocity at the region close to the boundary in the z-direction is significantly smaller due to the no-slip boundary conditions.

f34-04-9781455731411 — Fig. 34.4 Three-dimensional flow field for case 2: channel through-flow. Fig. 34.3b shows a cut of the velocity field at constant z. This graph shows the entire three-dimensional flow field. You can see that the by-flow extends along the z-axis but becomes smaller as z approaches the no-slip boundaries at the top and the bottom of the channel.

A Comment on the Pressure Field. The resulting pressure field in Fig. 34.4b looks slightly unsymmetri-cal, i.e., there is a blue line showing on the left-hand, top, and bottom side whereas there is no such line at the right-hand edge of the domain. This is only due to the fact that Maple discretizes the data during plotting. If you expand the view or check the numerical data, you will find that there is a blue line on the right-hand side as well.

34.4.4 Case 3: Step Expansion

The third case we will discuss is similar to the second case. This time, the microfluidic channel simply expands to a larger channel, where we again expect to see a symmetrical expansion of the flow field. The geometry is shown in Fig. 34.5a.

f34-05-9781455731411 — Fig. 34.5 Usage of the three-dimensional solver for case 3: step expansion. a) Geometry of the flow problem. A mi-crofluidic channel of width ω = 20 flows into a chamber of length L = 400 and width W = 100 μm. The height of both the chamber and the channel is H = 100 μm. For this case, we define a pressure drop of − 0.01 mbar across the domain. The stepwidth was h_xyz = 2 μm. The relaxation parameter was ω = 1.8 and the precision for both velocssity and pressure was set to 1 × 10^− 6. The fluid was chosen to be water with ρ = 1gcm^− 3and η = 1mPas. b) Numerical output obtained with implemented three-dimensional solver. The output is the velocity field in the xy-plane at z = 50 μm (see a) overlaid on top of the pressure field. The flow extends symmetrically into the chamber.

The function for setting up the boundary and initial values SetupBCIVCase3 is shown in listing 34.15. The only difference to SetupBCIVCase2 shown in listing 34.14 is the fact that the boundary condition is only set for x = 1 and not for x = L(see line 6). For x = L we define Neumann boundary conditions for the pressure and Dirichlet boundary conditions for the velocity (see line 14).

Listing 34.15

[Maple] Maple worksheet for solving case 3 using the three-dimensional solver. This listing is part of the listing ThreedimensionalFlowsNumerical.mw, which can be found in the supplementary material.

 1  SetupBCIVCase3:=proc(p0,entrywidth,L,W,H,dataVelx,dataVely,dataVelz,dataPress,kVolx,kVoly,kVolz, BCvel,BCpress) local x,y,z,halfwidth,halfchannelwidth:
 2  halfwidth:=W/2;
 3  halfchannelwidth:=entrywidth/2:
 4  for y from 1 to W do
 5   for z from 1 to H do
 6    if abs(halfwidth-y) < halfchannelwidth then
 7     BCvel[1,y,z]:= − 1:
 8     BCpress[1,y,z]:=1:
 9     dataPress[1,y,z]:=p0:
10    else
11     BCvel[1,y,z]:=1:
12     BCpress[1,y,z]:= − 1:
13    end if:
14    BCvel[L,y,z]:= − 1:
15    BCpress[L,y,z]:=1:
16   end do:
17  end do:
18  for x from 1 to L do
19   for y from 1 to W do
20    for z from 1 to H do
21     if y = 1 or y = W or z = 1 or z = H then
22     BCvel[x,y,z]:=1:
23     BCpress[x,y,z]:= − 1:
24    end if
25   end do:
26   end do:
27  end do:
28  return L,W,H,dataVelx,dataVely,dataVelz,dataPress,kVolx,kVoly,kVolz,BCvel,BCpress:ss
29 end proc:
30.
31 dataCase3:=SetupBCIVCase3(0.01,10,PrepareFlow3D(200,50,50)):
32 sol3:=SolveFlow3D(2,1,1,dataCase3,20000):
33
34 fieldplot([(x,y) → sol3[1](x,y,25),(x,y) → sol3[2](x,y,25)],1..200,1..50);
35 densityplot((x,y) → sol3[4](x,y,25),1..200,1..50,colorscheme = [“Blue”,“Green”,“Yellow”,“Red”]);
36 fieldplot3d([(x,y,z) → sol3[1](x,y,z),(x,y,z) → sol3[2](x,y,z),(x,y,z) → sol3[3](x,y,z) ],1..200,0..50,0..50,grid = [8,8,8],arrows = THICK);

Visualization. As before, we use a pressure drop of − 0.01 mbar across the x-direction of the domain and an entrywidth of 10 cells, which corresponds to a value of 20 μm. SolveFlow3D will require a total of 13 589 iterations to complete (see line 32). Fig. 34.5b shows the velocity field in the xy-plane at z = 50 pm, where we can again see the quickly expanding velocity field due to the fast momentum diffusion. Fig. 34.6 shows the three-dimensional view of the velocity field, where you can see the flow expansion even better when rotating the graph in Maple.

f34-06-9781455731411 — Fig. 34.6 Three-dimensional flow field for case 3: step expansion. Fig. 34.5b shows a cut of the velocity field at constant z. This graph shows the entire three-dimensional flow field. You can see that the by-flow extends along the z-axis but becomes smaller as z approaches the no-slip boundaries at the top and the bottom of the channel.

34.4.5 Case 4: Flow Around Objects

The fourth case we will discuss uses a classical rectangular microfluidic channel of width and height 100 pm. In the center of the channel, we place a cylindrical obstacle such as the pillars commonly used for deterministic lateral displacement devices, which are often referred to as microfluidic bumper arrays [3]. Obviously, this pillar introduces no-slip boundary conditions in the middle of the flow. We will see how the solver aligns the flow field in order to cope with this obstacle.

We first define a new function SetupBCIVCase4 shown in listing 34.16, which sets up our domain and the boundary conditions. The function accepts the pressure at the inlet p0 as well as the diameter d of the cylinder. The function then sets up the inlet and the outlet of our domain, i.e., for x = 1 and x = L as Neumann boundary condition for the velocity and as Dirichlet boundary condition for the pressure (see line 3). For all values of x it then sets the edges of the domain (y = 1 or y = w or z = 1 or z = H) as no-slip boundary conditions (see line 15). Whenever ${(\frac{W}{2} - y)}^{2} + {(\frac{L}{2} - x)}^{2} < {(\frac{d}{2})}^{2}$ ${(\frac{W}{2} - y)}^{2} + {(\frac{L}{2} - x)}^{2} < {(\frac{d}{2})}^{2}$ we have a pair of coordinates (x, y) that are inside of the obstacle (see line 19). At these locations, we set no-slip boundary conditions, i.e., Dirichlet boundary conditions for the velocity and Neumann boundary conditions for the pressure.

Listing 34.16

[Maple] Maple worksheet for solving case 4 using the three-dimensional solver. This listing is part of the listing ThreedimensionalFlowsNumerical.mw, which can be found in the supplementary material.

 1 SetupBCIVCase4:=proc(p0,d,L,W,H,dataVelx,dataVely,dataVelz,dataPress,kVolx,kVoly,kVolz,BCvel, BCpress) local x,y,z,radiusterm:
 2   radiusterm:=evalf((d/2)^2):
 3   for y from 1 to W do
 4    for z from 1 to H do
 5     BCvel[1,y,z]:= − 1:
 6     BCvel[L,y,z]:= − 1:
 7     BCpress[1,y,z]: = 1:
 8     dataPress[1,y,z]:=p0:
 9     BCpress[L,y,z]: = 1:
10    end do:
11   end do:
12   for x from 1 to L do
13    for y from 1 to W do
14     for z from 1 to H do
15      if y = 1 or y = W or z = 1 or z = H then
16       BCvel[x,y,z]:=1:
17       BCpress[x,y,z]:= − 1:
18      end if:
19      if (W/2 − y)^2 + (L/2 − x)^2 < radiusterm then
20       BCvel[x,y,z]:=1:
21       BCpress[x,y,z]:= − 1:
22      end if:
23     end do:
24    end do:
25   end do:
26   return L,W,H,dataVelx,dataVely,dataVelz,dataPress,kVolx,kVoly,kVolz,BCvel,BCpress:
27 end proc:
28 dataCase4:=SetupBCIVCase4(0.01,20,PrepareFlow3D(200,50,50)):
29 sol4:=SolveFlow3D(2,1,1,dataCase4,20000):
30 fieldplot([(x,y) → sol4[1](x,y,25),(x,y) → sol4[2](x,y,25)],1..200,1..50);
31 densityplot((x,y) → sol4[4](x,y,25),1..200,1..50,colorscheme = [“Blue”,“Green”,"Yellow","Red"]);
32 fieldplot3d([(x,y,z) → sol4[1](x,y,z),(x,y,z) → sol4[2](x,y,z),(x,y,z) → sol4[3](x,y,z) ],1..200,0..50,0..50,grid = [8,8,8],arrows = THICK);

Visualization. In this example we use a pillar diameter of 40 pm and a pressure drop of − 0.01 mbar across the x-direction (see line 29). The solver will require a total of 11302 iterations. Fig. 34.7b shows the velocity field in the xy-plane at z = 50 μm overlaid by the pressure field. Due to the scaling constrains our cylinder looks elliptical. However, what we can see very clearly is how the flow extends around the obstacle. Please also note the distortion in the pressure field. Looking at the three-dimensional view of the velocity field in Fig. 34.8 you can see the flow moving around the obstacle in the center.

f34-07-9781455731411 — Fig. 34.7 Usage of the three-dimensional solver for case 4: flow around objects. a) Geometry of the flow problem. A cylindrical object with diameter 40 μm is placed in the middle of a microfluidic channel. Such structures occur regularly in deterministic lateral displacement arrays. The microfluidic channel is of length L = 400 μm, width W = 100 pm, and height H = 100 μm. We define a pressure drop of− 0.01 mbar across the channel and use a stepwidth h_xyz = 2 pm. The relaxation parameter was ω = 1.8 and the precision for both velocity and pressure was set to 1 × 10^− 6. The fluid was chosen to be water with ρ = 1 g cm^− 3and η = 1 mPas. b) Numerical output obtained with implemented three-dimensional solver. The output is the velocity field in the xy-plane at z = 50 μm (see a) overlaid on top of the pressure field. As you can see, the flow extends around the object as expected.

f34-08-9781455731411 — Fig. 34.8 Three-dimensional flow field for case 4: flow around objects. Fig. 7b shows a cut of the velocity field at constant z. This graph shows the entire threedimensional flow field. You can see that the flow moves around the cylindrical object. Toward the top and the bottom as well as the sides of the channel it is additionally slowed due to the no-slip boundary conditions.

A Comment on the Pressure Field. In Fig. 34.7b we can see the obstacle showing up as the blue region in the pressure field, i.e., having a pressure of 0 mbar. Obviously the question arises: How can the pressure be 0 at these locations? The reason for this is that this area is a collection of cells defined as boundary values. Our solver currently only sets boundary values next to the cells of the domains that are updated during the iteration. It does not update the boundary values any further than just the neighboring cell. This means that if several boundary value cells are located next to each other, as is the case for the obstacle, only the cells on the edge of the boundary domain will be updated. Therefore the original values in the pressure data array are not changed and remain at their default value 0. Although this is physically incorrect there is no other way for our solver to act in this case. The problem is that the obstacle is a solid object. The physics of pressure distribution in a solid object is different from a fluid. Unless we implement additional routines that account for the physics of pressure guidance in solid objects, it would be physically incorrect to simply copy down the pressure value of one boundary value cell into the next. After all, the pressure on the left-hand side of the obstacle is different than the pressure on the right-hand side. Merely copying down the values from left to right would change the physics of the flow significantly and yield wrong results.

Therefore, the outmost boundary value “layer” of the object holds the correct surface pressure values. How this pressure is guided through the solid object must be calculated using different routines.

34.4.6 Case 5: Rectangular Channel Flow

The fifth case we will discuss uses a classical rectangular microfluidic channel of width and height 100 pm. Compared to the flow cases we have been calculating so far, this is actually a very easy case. However, it is interesting to discuss since it allows us to compare the analytical results obtained in section 16.4 with the results obtained using our custom-written solver. The function SetupBCIVCase5 that sets up our boundary and initial values is very compact. We simply define Neumann boundary values for the velocity and Dirichlet boundary conditions for the pressure at x = 1 and x = L. Across the domain, we define no-slip boundary conditions at the edges, i.e., for y = 1, y = W, z = 1, or z = H. The pressure drop is again − 0.01 mbar across the x-direction.

Visualization. The solver will require a total of 11664 iterations to converge. Fig. 34.9b shows the velocity field in the xy-plane at z = 50 pm overlaid with the pressure field. A three dimensional view is given in Fig. 34.10. As you can see, even though we only defined the pressure at the entrance and the exit of the domain, the solver determines a constant pressure drop across the channel. This is expected from the analytical solution, but it is interesting to see that the solver finds these values as well. Furthermore, the velocity is constant across the domain, which is (again) expected from the analytical solution assuming a fully developed flow profile. You can see the profiles in Fig. 34.12.

f34-09-9781455731411 — Fig. 34.9 Usage of the three-dimensional solver for case 5: rectangular channel flow. a) Geometry of the flow problem. The microfluidic channel is of length L = 400 μm, width W = 100 μm, and height H = 100 μm. We define a pressure drop of − 0.01 mbar across the channel and use a stepwidth h_xyz = 2 μm. The relaxation parameter was ω = 1.8 and the precision for both velocity and pressure was set to 1 × 10^− 6. The fluid was chosen to be water with ρ = 1gcm^− 3and η = 1 mPas. b) Numerical output obtained with implemented three-dimensional solver. The output is the velocity field in the xy-plane at z = 50 pm (see a) overlaid on top of the pressure field. We obtain the same linear pressure drop and channel flow profiles as obtained from the analytical solution.

f34-10-9781455731411 — Fig. 34.10 Three-dimensional flow field for case 5: rectangular channel flow.Fig. 34.9b shows a cut of the velocity field at constant z. This graph shows the entire three-dimensional flow field. We obtain the fully developed flow profiles we expect from rectangular cross-sections.

Going one step further, we plot the velocity distribution v_x (y, z) at y = 50 pm in direct comparison with the expected flow profile in a rectangular channel cross-section given by Eq. 16.53. As you can see from Fig. 34.11 the results are in good agreement. The analytical solution predicts a center velocity of 1.84mms^− 1 (expansion order of n_max = 10) whereas the numerical solver finds a value of 1.78mms^− 1.

f34-11-9781455731411 — Fig. 34.11 Comparison of the flow profiles in a rectangular channel cross-section found using the analytical solution (see Eq. 16.53) (dark blue, expansion order n_max = 10)and the numerical results found using the custom-implemented three-dimensional solver (medium green) using a channel width L = 400 μm and a channel height and width of W = H = 100 μm, a pressure drop of $- \frac{0.01}{0.4}$ $- \frac{0.01}{0.4}$ mbar mm^− 1and water as fluid. As you can see, there is a very good agreement between the analytical and the numerical solutions.

f34-12-9781455731411 — Fig. 34.12 Derived velocity v_x (x) (a) and pressure profile p (x) (b) along the x-axis for case 5 at y = 50 μm. As you can see, the solver recovers the expected constant velocity v_x and the linear pressure drop $\frac{d p}{d x}$ $\frac{d p}{d x}$ , which we expect from the analytical solution of the Poisson flow in rectangular channel cross-sections discussed in section 16.4.

It is important to compare the results obtained from a numerical solver with analytically determined values. As a matter of fact, case 5 is the only scenario that is sufficiently simple to allow an analytical solution. However, for this case, we find very good agreement of the analytical and the numerical solutions.

34.4.7 Case 6: Double-Fin Channel Flow

The sixth case we will discuss uses a classical rectangular microfluidic channel of width and height 100 pm into which two fin structures are introduced. The fins are 80 μm long, 75 pm wide, and 66.6 pm high. They are located at x = 0.1 L and x = 0.3 L. The first fin is attached to the front/bottom corner, whereas the second fin is attached to the back/top corner of the channel. The geometry is shown in Fig. 34.13a as an overview and in Fig. 34.13b as a view along the channel axis in the negative x-direction. The function for setting up the initial values and the boundary values is SetupBCIVCase6 and is shown in listing 34.18. For x = 1 and x = L we defined Neumann boundary values for the velocity and Dirichlet boundary values for the pressure (see line 4). Throughout the domain we define no-slip boundary values wherever y = 1, y = W, z = 1, or z = H (see line 14). Line 18 creates the first fin, and line 22 creates the second fin, both of which constitute obstacles with no-slip boundary conditions, i.e., Neumann boundary condition for the pressure and Dirichlet boundary condition for the velocity.

f34-13-9781455731411 — Fig. 34.13 Usage of the three-dimensional solver for case 6: double-fin channel flow. Geometry of the flow problem. The microfluidic channel is of length L = 400 μm, width W = 100 μm, and height H = 100 μm. We introduce two fins as obstacles into the channel each of which is 0.1 L long, 0.75 W wide, and 0.66 H high. The first fin is attached to the front/bottom, whereas the second one to the back/top corner of the channel. We use a pressure drop of − 0.01 mbar across the channel and a stepwidth of h_xyz = 2 μm. The relaxation parameter was ω = 1.8 and the precision for both velocity and pressure was set to 1 × 10^− 6. The fluid was chosen to be water with ρ = 1 gcm^− 3and η = 1 mPas. a) Bird’s eye view. b) View along the channel axis is negative x-direction.

Listing 34.17

[Maple] Maple worksheet for solving case 5 using the three-dimensional solver. This listing is part of the listing ThreedimensionalFlowsNumerical.mw, which can be found in the supplementary material.

 1  SetupBCIVCase5:=proc(p0,L,W,H,dataVelx,dataVely,dataVelz,dataPress,kVolx,kVoly,kVolz,BCvel,BCpress) local x, y, z:
 2  for y from 1 to W do
 3   for z from 1 to H do
 4    BCvel[1,y,z]:= − 1:
 5    BCvel[L,y,z]:= − 1:
 6    BCpress[1,y,z]: = 1:
 7    BCpress[L,y,z]: = 1:
 8    dataPress[1,y,z]:=p0:
 9   end do:
10  end do:
11  for x from 1 to L do
12   for y from 1 to W do
13    for z from 1 to H do
14     if y = 1 or y = W or z = 1 or z = H then
15      BCvel[x,y,z]:=1:
16      BCpress[x,y,z]:= − 1:
17     end if
18    end do:
19   end do:
20  end do:
21  return L,W,H,dataVelx,dataVely,dataVelz,dataPress,kVolx,kVoly,kVolz,BCvel,BCpress:
22 end proc:
23
24  dataCase5:=SetupBCIVCase5(0.01,PrepareFlow3D(200,50,50)):
25  sol5:=SolveFlow3D(2,1,1,dataCase5,20000):
26
27 fieldplot([(x,y) → sol5[1](x,y,25),(x,y) → sol5[2](x,y,25)],1..200,1..50);
28 densityplot((x,y) → sol5[4](x,y,25),1..200,1..50,colorscheme = ["Blue","Green","Yellow","Red"]);
29 fieldplot3d([(x,y,z) → sol5[1](x,y,z),(x,y,z) → sol5[2](x,y,z),(x,y,z) → sol5[3](x,y,z) ],1..200,0..50,0..50,grid = [8,8,8],arrows = THICK);

Listing 34.18

[Maple] Maple worksheet for solving case 6 using the three-dimensional solver. This listing is part of the listing ThreedimensionalFlowsNumerical.mw, which can be found in the supplementary material.

 1 SetupBCIVCase6:=proc(p0,L,W,H,dataVelx,dataVely,dataVelz,dataPress,kVolx,kVoly,kVolz,BCvel,BCpress) local x,y,z:
 2  for y from 1 to W do
 3   for z from 1 to H do
 4    BCvel[1,y,z]:= − 1:
 5    BCvel[L,y,z]:= − 1:
 6    BCpress[1,y,z]:=1:
 7    BCpress[L,y,z]:=1:
 8    dataPress[1,y,z]:=p0:
 9   end do:
10  end do:
11  for x from 1 to L do
12   for y from 1 to W do
13    for z from 1 to H do
14     if y = 1 or y = W or z = 1 or z = H then
15      BCvel[x,y,z]: = 1:
16      BCpress[x,y,z]:= − 1:
17     end if:
18     if x > L/5 and x < 2/5*L and y < 3/4*W and z < 2/3*H then
19      BCvel[x,y,z]: = 1:
20      BCpress[x,y,z]:= − 1:
21     end if:
22     if x > 3/5*L and x < 4/5*L and y > W/4 and z > 1/3*H then
23      BCvel[x,y,z]: = 1:
24      BCpress[x,y,z]:= − 1:
25     end if:
26    end do:
27   end do:
28  end do:
29  return L,W,H,dataVelx,dataVely,dataVelz,dataPress,kVolx,kVoly,kVolz,BCvel,BCpress:
30 end proc:
31
32 dataCase6:=SetupBCIVCase6(0.01,PrepareFlow3D(200,50,50)):
33 sol6:=SolveFlow3D(2,1,1,dataCase6,20000):
34
35 fieldplot([(x,y) → sol6[1](x,y,25),(x,y) → sol6[2](x,y,25)],1..200,1..50);
36 densityplot((x,y) → sol6[4](x,y,25),1..200,1..50,colorscheme = [“Blue”,“Green”,“Yellow”,“Red”]);
37 fieldplot([(x,z) → sol6[1](x,5,z),(x,z) → sol6[3](x,5,z)],1..200,1..50);
38 densityplot((x,z) → sol6[4](x,5,z),1..200,1..50,colorscheme = [“Blue”,“Green”,“Yellow”,“Red”]);
39 fieldplot3d([(x,y,z) → sol6[1](x,y,z),(x,y,z) → sol6[2](x,y,z),(x,y,z) → sol6[3](x,y,z)],1..200,0..50,0..50,grid = [8,8,8],arrows = THICK);

Visualization. Again, we use a pressure drop of − 0.01 mbar along the x-direction of the domain. The solver requires a total of 12 503 of iterations to complete. Fig. 34.14 shows the velocity field in the xy-plane at z = 50 μm overlaid by the pressure field. You can see the two fins showing up as blue areas in the plot. The flow moves around the fins, which accounts for the rather complex flow field. Fig. 34.15 shows the velocity field in the xz-plane at y = 10 μm, where we can again see the fins as areas of zero pressure. However, in this view, we can see the flow passing in front of the second fin.

f34-14-9781455731411 — Fig. 34.14 Numerical output for case 6 in the xy-plane at z = 50 μm.The image shows the velocity field overlaid on top of the pressure field. As you can see, the flow is very complex, avoiding the obstacles by flowing around it. The flow over the obstacles can be seen in Fig. 34.15.

f34-15-9781455731411 — Fig. 34.15 Numerical output for case 6 in the xz-plane at y = 10 μm. The image shows the velocity field overlaid on top of the pressure field. Again, you can see that the flow avoids the obstacles behind and before it. The flow around the obstacles can be seen in Fig. 34.14.

In order to get the most accurate view of the flow field, it is useful to inspect the three-dimensional views. The flow field is shown under two different angles in Fig. 34.16 where you can see the complex flow patterns around the two fin structures.

f34-16-9781455731411 — Fig. 34.16 Three-dimensional flow field for case 6:double-fin channel flow. The image shows two views of the three-dimensional flow fields at slightly different angles showing how the flow passes over and under as well as around the two obstacles. Please note that it is near-impossible to derive an analytical solution for this scenario.

As you can see, we can obtain solutions to flow problems as complicated as this within minutes using our custom-written numerical solver. Analytical solutions to these types of flows are near-impossible to derive. This shows the potential of using numerical solvers for similar flow problems.

34.4.8 Modifications

34.4.8.1 Volume Forces

As you can see, none of our examples have actually used volume forces; we merely left the array containing volume data empty. Feel free to experiment with volume forces, e.g., gravity, and see how they change the output of the solver. You will see, that in our microfluidic flow fields gravitational force can be neglected. We can also derive this by inspecting the Froude numbers (see section 9.9.11), which are very small, indicating that gravity can in fact be neglected. However, you may also apply different volume forces, e.g., electrokinetic forces or the like.

34.4.8.2 Reuse of Data

One nice aspect of the way we implemented the solver is that it can “reuse” previous data. As stated, it is useful to start with a lower number of maximum iterations in order to assess the speed at which the solver progresses. On very large computational domains, it may take several seconds to perform a total of 1000 iterations. However, we do not lose any time by progressing “in small steps.” If the solver performs a couple of iterations and then returns stating that the solution has not converged yet, we can simply execute the calculation again using the same data set. The solver will pick up the calculation on the nonconverged solution and make it better. This step can be repeated until we finally reach convergence. Once convergence is reached, SolveFlow3D will return instantaneously even if called multiple times because the convergence check will be positive every time. Therefore we cannot “over-process” the data.

By first running a couple of iterations and checking the total time required, we can estimate if the overall calculation will require minutes or rather hours to complete. This is useful because the way the code of our solver is currently written, it does not produce intermediary output to inform the user about the state of the calculation. Obviously, it is easy to add such a logging capability.

34.4.8.3 Saving and Loading Data

The fact that the calculations can be taken up using previously processed data also allows for another feature that would be a good addition: The capability to save and reload data. We would then first call the function PrepareFlow3D to set up the arrays and run several iterations on the data, whereupon we save the content of the arrays dataVelx, dataVely, dataVelz, and dataPress to file. Upon reloading, we simply call PrepareFlow3D again and load the saved data into the arrays. Obviously, it does not matter if we implement this piece of code in Maple or in C.

34.4.8.4 Stability Considerations

Another aspect about the reuse of data is stability considerations. Given the fact that finite difference method solvers tend to become instable if the gradients become too high we may come across a flow scenario that our solver will not be able to solve directly. This may happen if we apply, e.g., a pressure gradient that is too high. Another potential source of instability is choosing a stepwidth in space h_xyz that is too big. In any of these cases, the solver may abort with the message that the solver became instable.

In such cases, it may be possible for the solver to find a solution by relaxing the boundary conditions and starting up with, e.g., a smaller pressure gradient. This problem may produce a stable solution. Once this solution has been found, we can simply modify the boundary conditions by increasing the pressure gradient. We can do this by simply modifying the data in the pressure data dataPress. Usually, increasing the boundary conditions in steps of 10% is a good first guess. We then call the solver again, which has to only “adapt” the solution that is already in the data to match the new boundary conditions. This is equivalent to providing the solver with an approximated solution. As we have seen previously (see the function buildestimation in section 24.4.2.4) in many cases convergent solutions can be found if a suitable estimated solution is provided. By slowly increasing the boundary values until the desired value is reached, stable solutions can be obtained iteratively.

Obviously, this is a process that can be automated. As our solver reports an occurring numerical instability by returning a value of − 1 we could write a routine that slowly increases a boundary value checking for a convergent solution each time. In case a convergent solution is found, this solution is saved. The routine would then increase the boundary value and iterate this process. Wherever an instability is detected, the routine would reload the last successfully converged solution and choose a smaller increment in the boundary value until a new, “improved” convergent solution is found. Many commercial solvers implement routines that follow this protocol.

34.5 Summary

With this rather complex flow problem we will finish the application demonstration of our custom-written three-dimensional solver. As you have seen, one can implement a very powerful solver that readily solves very complex flow patterns without much effort. We have also seen that the solver produces results that are in good agreement with the analytical solutions if we can provide one. As you have seen from the examples we just discussed the cases where we can come up with analytical solutions are very limited and there are only a small number of scenarios in microfluidics that can be reduced to such simple flow cases. We have discussed the most important ones in section 15. For all other applications, it is necessary to refer to numerical solutions.

As this section demonstrated, there really is nothing special about numerical solvers. It is true that we have used Maple as an interface for inputting the data into our solver and displaying the results adequately. In fact, these are two very important points for a numerical solver software. Providing the user with means of inputting the geometry and the properties of the fluid, e.g., directly from computer aided design (CAD) systems make a solver significantly more comfortable to use than inputting the data the way our solver currently requires it. However, “under the hood” commercial numerical packages do nothing else than derive boundary conditions and set up the computational domain using suitable discretization. Obviously, a good interface for displaying the output is equally important. Here, we have used Maple routines that allow creating graphs of flow fields and the like. However, the actual work of solving the flow case was done by our custom-implemented routines.

..................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 34.3 Implementation of a Stationary Flow Numerical Solver

Create new playlist

Sign In

Sign Up

34.2.3.3 Instationary Flows

34.3 Implementation of a Stationary Flow Numerical Solver

34.3.1 Introduction

34.3.2 Boundary Values

34.3.3 Constants

34.3.4 Structures

34.3.5 Implementation

34.3.5.1 Main Listing

34.3.5.2 Block 1 - Initializing Variables and Constants

34.3.5.3 Block 2 - Update the Structures

34.3.5.4 Block 3 - Update the Velocities

34.3.5.5 Block 4 - Update the Pressure

34.3.5.6 The macrosUpdateFlowVariableDoNeumannandUpdateFlowVariableDoNeumannSingle

34.4 Usage of the Numerical Solver

34.4.1 Introduction

34.4.2 Case 1: Chamber By-Flow

34.4.3 Case 2: Chamber Through-Flow

34.4.4 Case 3: Step Expansion

34.4.5 Case 4: Flow Around Objects

34.4.6 Case 5: Rectangular Channel Flow

34.4.7 Case 6: Double-Fin Channel Flow

34.4.8 Modifications

34.4.8.1 Volume Forces

34.4.8.2 Reuse of Data

34.4.8.3 Saving and Loading Data

34.4.8.4 Stability Considerations

34.5 Summary

Table of Contents for
34.3 Implementation of a Stationary Flow Numerical Solver