**1. Introduction**

The fundamental physics of flows in microchannels are pivotal for the precise control of dynamic effects underlying transport phenomena such as momentum, mass or heat transfer and ultimately define the behavior or the system at hand for a practical application. As an example, establishing the correct flow physics is crucial for modeling and understanding the transport kinetics of an analyte from the bulk of the flow to the surface of, e.g., a biosensor [1]. In particular, the response dynamics of a biosensor strongly depends on the di ffusion kinetics of the analyte within the bulk of the flow of the biosensor. The formation of the Nernst di ffusion layer is the limiting factor that defines the transport dynamics and thus the dynamic response of the biosensor [2]. This is particularly important during the initiation of the flow and in transitions of the flow such as, e.g., at the early stage of an experiment. These are tightly interlinked with the fluid mechanics of the system and thus a detailed understanding, specifically of transient e ffects in the bulk fluid flow, is pivotal for a correct system assessment [3]. In many respects, the fluid physics of microfluidics are straight-forward. Given the low Reynolds number flows commonly encounter in microfluidics, e ffects such as turbulence are rarely relevant in microfluidic flows [4].

However, given the fact that the underlying equation of fluid dynamics, i.e., the Navier–Stokes equation, is very di fficult to solve even in (seemingly) simple channel geometries and flow conditions, researchers commonly refer to numerical methods to model the flow physics [5]. However, correctly applying a complex numerical modeling or solver software package is usually beyond the scope of application-driven research in microfluidics and adjacent fields.

We have recently demonstrated, that the simplified Navier–Stokes equation for laminar flow can be conveniently solved in a spreadsheet analysis software such as Microsoft Excel [6,7]. Spreadsheet tools are convenient platforms for implementing schemes in an intuitive and documentable manner. Software packages such as Microsoft Excel are widely available, and researchers are used to working with these packages from their daily work routine. Thus, applying a scheme which is implemented on this platform comes with a significantly lower barrier than adopting a complex and costly numerical solver package.

In our two most recent contributions, we demonstrated that the Poiseuille equation, a simplified version of the Navier–Stokes equation for stationary laminar flows, can be implemented in a wide range of channel geometries. For geometries for which analytical solutions can be obtained, the numerical results from the spreadsheet deviate only by a few percent from the analytical solutions [6]. We also demonstrated that the spreadsheet can be used to implement di fferent boundary conditions besides the commonly employed no-slip (Dirichlet) boundary condition such as, e.g., Neumann-type boundary conditions which occur on slip surfaces as well as on open channel geometries [7]. These solutions are close-to-exact to the analytical solutions, which again, can only be derived for very simple channel geometries, usually with a high degree of symmetry such as, e.g., in circular channel (Hagen–Poiseuille) flows.

However, all of these solutions were derived for stationary flow scenarios, i.e., flows for which the time-dependency of the Poiseuille equation is ignored. This case is applicable for applications where the flow is assumed to be in steady-state and flow acceleration (i.e., during initiation of the flow by applying a pressure drop), as well as flow retardation (i.e., during stopping of the flow by removing the driving pressure drop) are not taken into account. However, there are many applications where time-dependence needs to be considered, especially in applications where the momentum di ffusion is overlaid by, e.g., mass di ffusion. In this paper we will show that the spreadsheets can be designed to reflect time dependency, allowing the study of transient e ffects during flow initiation and retardation, as well as intermediate changes in the driving pressure drop which modifies the flow conditions.

## **2. Numerical Scheme**

#### *2.1. Navier–Stokes Equation for Time-Dependent Flow*

In comparison to the stationary version of the Navier–Stokes equation used in [6], we now have to consider time-dependency. Neglecting volume forces, the Navier–Stokes equation simplifies to

$$
\rho \frac{\partial \overrightarrow{v}}{\partial t} = -\overrightarrow{\nabla}p + \eta \Delta \overrightarrow{v} \tag{1}
$$

which includes pressure contribution ( <sup>∇</sup>*p*), momentum di ffusion (ηΔ → *v* ), as well as time-dependency ( ∂ → *v* ∂*t* ). Again, we consider parallel flow which reduces the dependent variable → *v* (a vector) to only the contribution *vx* along the *x*-axis. Equation (1) can therefore be written as

→

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

where we also note that the only driving pressure drop is the drop along the x-axis of the channel given by Δ*p* Δ*L*.

#### *2.2. Numerical Scheme for the Second-Order Partial Di*ff*erential Equations*

In [6] we derived numerical schemes for second-order partial differentials from the Taylor series given by

$$\frac{d^2f}{d\mathbf{x}^2}(\mathbf{x}\_0) = \frac{f(\mathbf{x}\_0 + \Delta \mathbf{x}) + f(\mathbf{x}\_0 - \Delta \mathbf{x}) - 2f(\mathbf{x}\_0)}{\left(\Delta \mathbf{x}\right)^2} \tag{3}$$

which allowed us to rewrite the right-hand side of Equation (2) to

$$\begin{aligned} \rho \frac{\partial v\_x}{\partial t} &= -\frac{\Delta p}{\Delta L} + \\ \eta \frac{v\_x(y\_0 + \Lambda y\_x z\_0) + v\_x(y\_0 - \Lambda y\_x z\_0) - 2v\_x(y\_0, z\_0)}{\left(\Lambda y\right)^2} &+ \frac{v\_x(y\_0, z\_0 + \Lambda z) + v\_x(y\_0, z\_0 - \Lambda z) - 2v\_x(y\_0, z\_0)}{\left(\Lambda z\right)^2} \end{aligned} \tag{4}$$

which we can further simplify by using a common step width in space *hyz* = Δ*y* = Δ*z* to yield

$$\begin{aligned} \rho \frac{\partial v\_x}{\partial t} &= -\frac{\Delta p}{\Delta L} + \\ \eta \left( \frac{v\_x \left( y\_0 + h\_{yx} z\_0 \right) + v\_x \left( y\_0 - h\_{yx} z\_0 \right) + v\_x \left( y\_0 z\_0 + h\_{yx} \right) + v\_x \left( y\_0 z\_0 - h\_{yx} \right) - 4v\_x \left( y\_0 z\_0 \right)}{h\_{yx}^2} \right) \end{aligned} \tag{5}$$

This scheme uses a so-called finite difference scheme to approximate the second order partial differential by considering the changes in the function over a finite difference in the independent variables *y* and *z*. This scheme is a second-order scheme and thus numerically very stable.

Introducing the general notation *F*(*<sup>t</sup>*,*y*,*<sup>z</sup>*) for the value of *vx* at the position (*y*0, *<sup>z</sup>*0) at time point *t*, *F*(*<sup>t</sup>*,*y*+1,*z*) for the value of *vx* at the position (*y*0 + *hyz*, *<sup>z</sup>*0) at time point *t*, *F*(*<sup>t</sup>*,*y*−1,*z*) for the value of *vx* at the position (*y*0 − *hyz*, *<sup>z</sup>*0) at time point *t*, *F*(*<sup>t</sup>*,*y*,*z*+<sup>1</sup>) for the value of *vx* at the position (*y*0, *z*0 + *hyz*) at time point *t* and *F*(*<sup>t</sup>*,*y*,*z*−<sup>1</sup>) for the value of *vx* at the position (*y*0, *z*0 − *hyz*) at time point *t* we can rewrite Equation (5) to

$$\frac{\partial v\_{\chi}}{\partial t} = -\frac{\Delta p}{\Delta L} + \eta \left( \frac{F^{(t, y+1, z)} + F^{(t, y-1, z)} + F^{(t, y, z+1)} + F^{(t, y, z-1)} - 4F^{(t, y, z)}}{h\_{yz}^2} \right)$$

$$\frac{\partial v\_{\chi}}{\partial t} = -\frac{1}{\rho} \frac{\Delta p}{\Delta L} + \frac{\eta}{\rho} \left( \frac{F^{(t, y+1, z)} + F^{(t, y-1, z)} + F^{(t, y, z+1)} + F^{(t, y, z-1)} - 4F^{(t, y, z)}}{h\_{yz}^2} \right) \tag{6}$$

#### *2.3. Numerical Scheme for the First-Order Partial Di*ff*erential with Respect to Time*

The second partial differential required for setting up the numerical scheme is the first-order partial differential with respect to time. For this we expand the Fourier series in the positive direction given by

$$f(\mathbf{x}\_0 + \Delta \mathbf{x}) = \sum\_{n=0}^{n\_{\text{max}}} \frac{1}{n!} \left. \frac{d^n f}{d \mathbf{x}^n} \right|\_{\mathbf{x}\_0} \Delta \mathbf{x}^n + O\_{n\_{\text{max}} + 1} \tag{7}$$

to *n*max = 1, where *n*max is the "expansion order" and *On*max is the error function of order *n*max. We subsequently obtain

$$f(\mathbf{x}\_0 + \Delta \mathbf{x}) = f(\mathbf{x}\_0) + \frac{df}{d\mathbf{x}} \Delta \mathbf{x} + O\_2$$

$$\frac{df}{d\mathbf{x}} = \frac{f(\mathbf{x}\_0 + \Delta \mathbf{x}) - f(\mathbf{x}\_0)}{\Delta \mathbf{x}}\tag{8}$$

which is a first-order approximation for the first derivative in time (a so-called forward Euler scheme). Combining Equations (8) and (6) results in the numerical scheme

$$\frac{F^{(t+1,y,z)} - F^{(t,y,z)}}{h\_t} = -\frac{1}{\rho} \frac{\Delta p}{\Delta L} + \frac{\eta}{\rho} \left( \frac{F^{(t,y+1,z)} + F^{(t,y-1,z)} + F^{(t,y,z+1)} + F^{(t,y,z-1)} - 4F^{(t,y,z)}}{h\_{yz}^2} \right)$$

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

$$F^{(t+1,y,z)} = F^{(t,y,z)} - \frac{h\_t}{\rho} \frac{\Delta p}{\Delta L} + \frac{h\_{t'}\eta}{\rho} \left( \frac{F^{(t,y+1,z)} + F^{(t,y-1,z)} + F^{(t,y,z+1)} + F^{(t,y,z-1)} + F^{(t,y,z-1)}}{h\_{yz}^2} \right)$$

$$F^{(t+1,y,z)} = F^{(t,y,z)} \left( 1 - 4\frac{h\_{t'}\eta}{\rho h\_{yz}^2} \right) + \frac{h\_{t'}\eta}{\rho} \left( \frac{F^{(t,y+1,z)} + F^{(t,y-1,z)} + F^{(t,y,z+1)} + F^{(t,y,z-1)}}{h\_{yz}^2} \right) - \frac{h\_t}{\rho} \frac{\Delta p}{\Delta L}$$

$$F^{(t+1,y,z)} = F^{(t,y,z)} \left( 1 - 4\Omega \right) + \Omega \left( F^{(t,y+1,z)} + F^{(t,y-1,z)} + F^{(t,y,z+1)} + F^{(t,y,z+1)} \right) - \Gamma \tag{9}$$

with Ω = *ht*Δη <sup>ρ</sup>*h*2*yz* and Γ = *ht* ρ Δ*p* Δ*L* , *ht* being the step width in time. Equation (9) is a second-order scheme with respect to space and a first-order scheme with respect to time. It allows stepping forward in time from a known value *F*(*<sup>t</sup>*,*y*,*<sup>z</sup>*) at the location (*x*, *y*) at time *t* to the unknown value *F*(*<sup>t</sup>*+1,*y*,*<sup>z</sup>*) at the location (*x*, *y*) at time *t* + 1. Compared to the schemes used in [6] this scheme is not iterative as the solution at a given time point t does not have to fulfill the Navier–Stokes equation but only Equation (8). One point of note is the fact that the numerical scheme in time is only first-order. In order to not risk numerical instability the step width *ht* in time must be chosen su fficiently small. In general, first-order approximations are numerically significantly less stable than higher-order implementations but they are computationally very cost-e ffective and thus very simple to implement. The numerical stability of the overall scheme hinges mostly on the step width in time.
