*4.1. Derivation*

In order to verify the correctness of the numerical results obtained from the implemented solver in Microsoft Excel, we chose the case of initiating two-dimensional flow in a rectangular channel cross-section with the no-slip boundary condition as an example. In this scenario, the flow in the channel is originally at rest. At *t* = 0, a driving pressure gradient is applied along the length of the channel thus initiating the flow whereby the characteristic velocity profile in a rectangular channel cross-section Poiseuille flow will form. The dynamics of the evolution of this flow profile will be studied using the derived analytical solution which yields the exact results. These will then be compared to the numerical results obtained from the spreadsheet.

The relevant partial di fferential equation for this case is Equation (2) which is rewritten to

$$\frac{\rho}{\eta} \frac{\partial v\_x}{\partial t} - \left(\frac{\partial^2 v\_x}{\partial y^2} + \frac{\partial^2 v\_x}{\partial z^2}\right) = -\frac{\Delta p}{\eta \Delta L} \tag{10}$$

We assume the solution to consist of a steady-state component which is time-independent and a transient solution which is time-dependent. The latter will be dominating during the initiation of the flow and decline as the flow achieves steady-state. Thus, the general solution will be

$$v\_{\mathbf{x}}(t, y, z) = v\_{\mathbf{x}, \text{steady}-\text{state}}(y, z) + v\_{\mathbf{x}, \text{transient}}(t, y, z) \tag{11}$$

If inserted into Equation (10) we obtain

$$\begin{split} \frac{\partial}{\partial \eta} \Big( \frac{v\_{\text{x,tramion}}}{\partial t} \Big) - \quad & \left( \frac{\partial^2 v\_{\text{x,tramly}-\text{stat}}}{\partial y^2} + \frac{\partial^2 v\_{\text{x,tramly}-\text{stat}}}{\partial z^2} + \frac{\partial^2 v\_{\text{x,tramion}}}{\partial y^2} + \frac{\partial^2 v\_{\text{x,tramion}}}{\partial z^2} \right) \\ & = -\frac{\Delta p}{\eta \Delta L} \end{split} \tag{12}$$

Please note that <sup>∂</sup>*vx*,steady−state ∂*t* = 0. The steady-state solution *vx*,steady−state for this flow case is known and was derived in a similar case in Richter et al. [7] as Equation 13' using an *eigenvalue* expansion for mixed boundary cases. The solution for no-slip boundary conditions is derived elsewhere ([8], Equation 16.62)

$$\begin{split} w\_{\text{x,steady-state}}(y, z) &= -\frac{16}{\eta \pi^2} \frac{\Delta p}{\Delta L} \sum\_{n=0}^{\otimes 0} \sum\_{m=0}^{\otimes 0} \sin \left( (2n+1)\pi \frac{y}{W} \right) \sin \left( (2m+1)\pi \frac{z}{H} \right) \\ & \left. \left( (2n+1)(2m+1) \left( \left( \frac{(2n+1)\pi}{W} \right)^2 + \left( \frac{(2m+1)\pi}{H} \right)^2 \right) \right)^{-1} \end{split} \tag{13}$$

The second-order partial differentials required for Equation (12) are given by

$$\begin{split} \frac{\partial^2 v\_{z, \text{steady-state}}}{\partial y^2} &= \frac{16}{\eta \cdot \pi^2} \frac{\Delta p}{\Delta L} \sum\_{n=0}^{\infty} \sum\_{m=0}^{\infty} \left( \frac{(2n+1)\pi}{W} \right)^2 \sin \left( (2n+1)\pi \frac{y}{W} \right) \sin \left( (2m+1)\pi \frac{z}{H} \right) \\ & \left( (2n+1)(2m+1) \left( \left( \frac{(2n+1)\pi}{W} \right)^2 + \left( \frac{(2m+1)\pi}{H} \right)^2 \right) \right)^{-1} \end{split} \tag{14}$$

$$\frac{\partial^2 v\_{x,\text{steady-state}}}{\partial z^2} = \frac{16}{\eta} \frac{\Delta p}{\Delta t} \sum\_{n=0}^{\infty} \sum\_{m=0}^{\infty} \left( \frac{(2n+1)\pi}{H} \right)^2 \sin \left( (2n+1)\pi \frac{y}{W} \right) \sin \left( (2m+1)\pi \frac{z}{H} \right)$$

$$\left( (2n+1)(2m+1) \left( \left( \frac{(2n+1)\pi}{W} \right)^2 + \left( \frac{(2m+1)\pi}{H} \right)^2 \right) \right)^{-1} \tag{15}$$

$$\frac{\partial^2 v\_{x,\text{steady-state}}}{\partial y^2} + \frac{\partial^2 v\_{x,\text{steady-state}}}{\partial z^2} = \frac{16}{\eta} \frac{\Delta p}{\Delta L} $$

$$\sum\_{n=0}^{\infty} \sum\_{m=0}^{\infty} \sin \left( (2n+1)\pi \frac{y}{W} \right) \sin \left( (2m+1)\pi \frac{z}{H} \right) \left( (2n+1)(2m+1) \right)^{-1} \tag{16}$$

Equation (16) is the representation of a constant by a two-dimensional Fourier series. Details on this can be found elsewhere ([8], Equation 4.44). We can thus simplify Equation (16) to

$$\frac{\partial^2 \upsilon\_{\text{x,steady-state}}}{\partial y^2} + \frac{\partial^2 \upsilon\_{\text{x,steady-state}}}{\partial z^2} = \frac{\Delta p}{\eta \Delta L} \tag{17}$$

Using Equation (17), we can rewrite Equation (12) to

$$\frac{\partial}{\partial t}\frac{v\_{\text{x,transient}}}{\partial t} - \frac{\partial^2 v\_{\text{x,transient}}}{\partial y^2} - \frac{\partial^2 v\_{\text{x,transient}}}{\partial z^2} = 0 \tag{18}$$

Equation (18) is a homogeneous partial differential equation which can be solved by a separation of variables approach using

$$w\_{\text{x,transient}} = T(t)\chi(y)Z(z) \tag{19}$$

Inserting Equation (19) into Equation (18) yields

$$\frac{\rho}{\eta} \frac{1}{T} \frac{d\mathbf{T}}{dt} - \frac{1}{Y} \frac{d^2 \mathbf{Y}}{dy^2} - \frac{1}{Z} \frac{d^2 \mathbf{Z}}{dz^2} = 0 \tag{20}$$

from which we derive the two second-order ordinary differential equations

$$\frac{1}{Y}\frac{d^2Y}{dy^2} = -\lambda\_Y^2\tag{21}$$

$$\frac{1}{Z}\frac{d^2Z}{dz^2} = -\lambda\_Z^2\tag{22}$$

and one first-order differential equation

$$\frac{\rho}{\eta} \frac{1}{T} \frac{\text{dT}}{dt} = - \left(\lambda\_Y^2 + \lambda\_Z^2\right) \tag{23}$$

*Biosensors* **2019**, *9*, 67

> The solution to Equation (23) is straight-forward and given by integration as

$$T(t) = C\_0 \, \text{g}^{-(\lambda\_Y^2 + \lambda\_Z^2)\frac{\eta}{\rho}} \, t \tag{24}$$

The solutions to Equations (21) and (22) are given by the *eigenfunctions* and derived in a similar fashion as shown in [7] to be

$$\mathcal{Y}(y) = \sum\_{n=0}^{\infty} \mathbb{C}\_n \sin(n\pi \frac{y}{W}) \,\,\lambda\_Y = \frac{n\pi}{W} \tag{25}$$

$$Z(z) = \sum\_{m=0}^{\infty} \mathcal{C}\_m \sin \left( m\pi \frac{z}{H} \right), \lambda\_Z = \frac{m\pi}{H} \tag{26}$$

Inserting Equations (24), (25) and (26) into Equation (19) yields the transient solution as

$$v\_{\text{x,transient}}(t,y,z) = e^{-(\lambda\_Y^2 + \lambda\_Z^2)\frac{\pi}{\rho}t} \sum\_{n=0}^{\infty} \sum\_{m=0}^{\infty} \mathbb{C}\_{mn} \sin\left(n\pi\frac{y}{W}\right) \sin\left(m\pi\frac{z}{H}\right) \tag{27}$$

We still lack the constant *Cnm*, which we can derive from the initial condition of the accelerating flow which requires the flow to be constant, i.e., *vx* = 0. This is the case for *vx*,transient(*<sup>t</sup>* = 0, *y*, *z*) = <sup>−</sup>*vx*,steady−state(0, *y*, *<sup>z</sup>*). In this case we find

$$\begin{aligned} \sum\_{n=0}^{\infty} \sum\_{m=0}^{\infty} C\_{nm} \quad \sin\left((2n+1)\pi\frac{\mathsf{y}}{W}\right) \sin\left((2m+1)\pi\frac{\mathsf{z}}{H}\right) \\ = \frac{16}{\eta\cdot\pi^2} \frac{\Delta p}{\Delta L} \sum\_{n=0}^{\infty} \sum\_{m=0}^{\infty} \sin\left((2n+1)\pi\frac{\mathsf{y}}{W}\right) \sin\left((2m+1)\pi\frac{\mathsf{z}}{H}\right) \\ = \left((2n+1)(2m+1)\left(\left(\frac{(2n+1)\pi}{W}\right)^2 + \left(\frac{(2m+1)\pi}{H}\right)^2\right)\right)^{-1} \end{aligned} \tag{28}$$

from which we find

$$C\_{nm} = \frac{16}{\eta \cdot \pi^2} \frac{\Delta p}{\Delta L} \Big| (2n+1)(2m+1) \Big| \left( \frac{(2n+1)\pi}{W} \right)^2 + \left( \frac{(2m+1)\pi}{H} \right)^2 \Big| \right)^{-1} \tag{29}$$

Assembling Equation (11) from Equations (13), (27) and (29) we find

$$\begin{split} v\_{x}(t,y,z) &= -\frac{16}{\eta \pi^{2}} \frac{\Delta p}{\Delta L} \sum\_{n=0}^{\infty} \sum\_{m=0}^{\infty} \begin{pmatrix} 1 - e^{-\left(\left(\frac{(2n+1)\pi}{W}\right)^{2} + \left(\frac{(2m+1)\pi}{H}\right)^{2}\right)} \frac{\eta}{\rho}t \\ \sin((2n+1)\pi \frac{y}{W}) \sin\left((2m+1)\pi \frac{z}{H}\right) \\ \end{pmatrix} \\ & \left((2n+1)(2m+1) \left[\left(\frac{(2n+1)\pi}{W}\right)^{2} + \left(\frac{(2m+1)\pi}{H}\right)^{2}\right]\right)^{-1} \end{split} \tag{30}$$

Visualization

In the following, Equation (30) is visualized using a microfluidic channel with 100 μm width and 100 μm height and choosing water (η = 1 mPa·s, ρ = 1 g/cm3) as the fluid. The Fourier series in Equation (30) is expanded to *n*max = *m*max = 10. Figure 2 shows the calculated profiles for 10 μs, 100 μs and 1000 μs. As can be seen the fluid is originally at rest. Over time, the characteristic profile of the Poiseuille flow in a rectangular channel cross-section evolves. By inspection of Equations (27) and (30), this behavior is expected. As the solution consists of a transient and a steady-state component, the first of which decays exponentially over time, the steady-state solution is effectively modulated by the exponential decay of the time-dependent contribution of the transient solution.

**Figure 2.** Velocity profile calculated from the analytical solution of Equation (30) for water as fluid and a channel with 100 μm width and 100 μm height. The Fourier series is expanded to *n*max = *m*max = 10. The expansion order of the Fourier series was 10 both along *y* and *z*. Velocity profiles for different points in time are plotted: 10 μs (dark blue), 100 μs (medium green) and 1000 μs (light yellow).

#### *4.2. Application of the Derived Spreadsheet*

#### 4.2.1. Initiating Two-Dimensional Flow in Rectangular Channel Cross-Sections

Given the analytical solution Equation (30) we now proceed to solving the same case using the numerical solver implemented in the spreadsheet. The step width in space was chosen as *hyz* = 2.5 μm and the step width in time as *ht* = 1 μs. Again, water as the fluid was assumed and thus η = 1 mPa·s, ρ = 1 g/cm<sup>3</sup> were set in the "Variables" section of the spreadsheet. For direct comparison with Figure 2, the spreadsheet was iterated 100 times to yield the velocity profile at *t* = 100 μs and *t* = 1000 μs. The resulting velocity profile was then compared with the analytical Page: 10 solution given by Equation (30). Figure 3 shows the relative error of the numerical output compared to the analytical solution. The maximum error in the whole computational domain is below 3% relative error in both cases. The errors are strongest in regions of steep gradients in the velocity profiles, i.e., at the edges. In order to decrease the error further, the computational domain can be increased by increasing the number of cells. However, even on this rather coarse and compact domain, the errors fall within acceptable limits.

**Figure 3.** Comparison of the results obtained from the numerical scheme implemented in Microsoft Excel against the analytical solution of Equation (30) for time point 100 μs (**a**) and 1000 μs (**b**). The graphs show the relative error. Both cases show very good agreements with maximum errors below 3%.

4.2.2. Complex Flow Cases: Di fferent Channel Cross-Sections

As shown, the analytical solution can be derived for the (rather simplistic) case of rectangular channel cross-sections. However, the availability of these solutions is limited if the channel geometries become more complex. However, these cross-sections are straight-forward to implement in the spreadsheet. Figure 4a shows the example of a rectangular channel with two fins constricting the flow.

**Figure 4.** Complex flow cases involving intricate domains such as a double-constricted channel cross-section. (**a**) These are implemented by adding additional boundary values, i.e., by copying the grey cell values which have static values. The image shows the flow after around 100 iterations. (**b**) By simple overwriting portions of the boundary values with domain values, i.e., reimplementing the conditional copy operation of the center values, one of the obstacles is removed at a given point in time. The scheme can then be iterated to see how the flow adapts to the changes in boundary conditions.

These can be implemented by copying additional boundary values into the center panel. By overwriting the formulae in these cells, the numerical scheme cannot iterate these values. However, they influence neighboring cells and thus show up as regions with no flow in the iterated panel on the right. Again, starting with a static flow, the scheme can be used to observe the development of the velocity profile over time. Figure 4a shows the profile at iteration #113, i.e., at *t* = 113 μs. Boundary conditions can be changed dynamically. As an example, the lower constriction in the channel of Figure 4a was removed at *t* = 113 μs. Figure 4b shows the adjusting flow 74 iterations after the object was removed. As can be seen, the flow re-expands into the regions blocked previously. Removing a boundary can be accomplished by simply recopying the formulae from the domain into the center panel thus reactivating the iteration.

#### 4.2.3. Boundary Conditions and Initial Conditions

As discussed, boundary conditions with fixed values (Dirichlet-type) can be implemented by overwriting cells in the center panel with fixed values. No-slip boundary conditions are implemented by setting this value to 0. Surfaces with fixed velocity as in the case of, e.g., Couette flow, can be implemented by setting the velocity of these boundary cells to a fixed non-zero value. As demonstrated in [7], it is possible to also implement Neumann-type boundary conditions by simply setting the cell values of the boundary to the value of the adjacent cell within the computational domain thus effectively generating a velocity gradient of 0 as required for, e.g., slip boundary conditions. Obviously, the gradient can also be set to a fixed value by simply setting the value of the boundary cell to the value of the neighboring cell within the computational domain plus a fixed o ffset value.

It is also possible to use non-zero initial values for computation. Until now, we assumed the flow to be static at the beginning of the experiment and to increase due to the application of the flow as a consequence of the application of a pressure gradient. This spreadsheet analyses the declining velocity profile of an initially moving fluid as a consequence of the removal of the driving pressure. The initial values are copied from a spreadsheet that iterated until the velocity profile was fully developed. The spreadsheet of Figure 5 sets the driving pressure gradient to 0. As you can see, after 1000 iterations the velocity profile has smeared slightly but the overall distribution is still similar to the case of the fully-developed profile. This is to be expected from the analytical solution to Equation (27). For a decelerating flow, the steady-state fully-developed flow profile would be modulated by an exponential term, thus generating an exponential decay. As a consequence, the flow profiles would simply be scaled but the shape of the profile would remain largely intact. As can be seen from Figure 5, after 1000 iterations, the maximum velocity in the channel has decreased to a value of below 0.5 mm/s.

**Figure 5.** Example of a flow case using boundary conditions. The demonstrated case is a decelerating flow which was initially fully-defined as indicated by the initial conditions. At *t* = 0 the driving pressure gradient is removed and the flow begins to stagnate. After more than 1500 iterations, the velocity profile is smeared and the velocity profile drops to maximum values below 0.5 mm/s.
