Finite Difference Method
30.1 Introduction
As we have seen, weighted residual methods form a class of methods that can be used to solve differential equations. They make use of approximation functions that are appropriately weighted in order to find a solution which approximates the solution to the differential equations as closely as possible. Weighted residual methods are used in several other commonly encountered methods for solving differential equations numerically, mostly FVM and FEM.
Before introducing these methods, we will quickly introduce the finite difference method (FDM). We have (indirectly) worked with FDM already, as it forms the basis for most of the numerical schemes we used in section 25.
The concept of FDM is focused on approximating differentials. In contrast to this, weighted residual methods evaluate the integral of a differential equation and optimize an approximation such that the integrals of the correct and the approximated solutions match on a given domain. Therefore, these equations use integral approximations. FDM uses an approximation of the differential of the differential equation. It therefore is a differential approximation. The mathematics of FDM is based on Taylor series approximations (see section 4.2). The most common equations are
• Central finite difference schemes (see Eq. 4.6) for approximating first derivatives
• Forward finite difference schemes (see Eq. 4.7) for approximating first derivatives
• Backward finite difference schemes (see Eq. 4.8) for approximating first derivatives
• Central finite difference schemes (see Eq. 4.10) for approximating second derivatives
Obviously, these schemes are used in many forms in numerical solvers and we have discussed many numerical methods which use these approximations in various forms. For details on these methods, please refer to section 27.2. The difference in the solution, i.e., the finite change of the solution is approximated on a very small finite interval using one of these equations. All of these equations are linear, i.e., the solution is linearly approximated. Obviously, this approximation is only correct if the interval on which the function is linearized is sufficiently small. Otherwise, the solution becomes inexact.
30.2 Advantages and Disadvantages
The main advantage of FDM is the fact that it is a very exact method. The solutions obtained from FDM are usually significantly closer to the solution, as the results obtained from, e.g., weighted residual methods. This additional precision comes with several disadvantages.
Ordered Grids. As we have already seen when using FDM in section 27.2, we must split the domain on which we perform the linear approximation of the solution. The increments, which are usually referred to as the stepwidth in space, must be sufficiently small in order for the solution to be exact. If the domain is large and geometrically complex, this may amount to a very high number of intervals on which we need to operate. If a domain can be described by a simple geometry, e.g., a simple rectangular channel, the discretization is straightforward as the domain is simply split into n intervals of stepwidth h. Usually the size of the stepwidth is defined by the smallest relevant geometrical feature. As an example, for a microfluidic channel of width 100 µ m and height 50 µ m we could chose the stepwidth to be h = 5 µ m both along the channel width and height and end up with a pretty good result. However, if the microfluidic channel features groove structures of around 1 µ m width and depth as common, e.g., mixer structures, we would need to refine the stepsize. If the flow in the groove should be resolved, we would require a stepwidth of around h = 0.1 µ m potentially even less.
For FDM the smallest feature size defines the stepwidth h locally. However, the method relies on a certain order in the grid. In order to calculate the central finite difference, we require a cell directly before and after the cell we are currently working on along both directions. This implies that the grid must be ordered, i.e., each grid has two neighbors in each direction. The cells are arranged in a parallel fashion and there is no local adaptation of the cells to suit, e.g., the curvature of the channel wall. This means that a curved surface must always be approximated, i.e., it becomes “pixelated.” You can see an example of this process in Fig. 30.6a. The only method of decreasing this pixelation is by using smaller stepwidths.
Another important implication of the mathematics underlying FDM is the fact that neighboring cells should usually the same stepwidth h. This means that the smallest geometrical feature size of the domain will determine the stepwidth and therefore, the number of cells in the domain. Returning to our example with the microfluidic channel without the groove structure, we determined that a grid with 200 cells would be sufficient to obtain a good solution. If the groove is added to the channel, the stepwidth decreases to 0.1 µ m and the total number of cells increases to 500 000 which is a significant boost in complexity.
Numerical Stability. Another point we need to consider is the fact that finite difference schemes have certain limits to the numerical stability. We have discussed problems with numerical stability of the forward finite difference scheme when discussing the forward Euler method (see section 27.2.2). We have also discussed numerous variations of these numerical methods with significantly better numerical stability. One of these examples is the backward Euler method (see section 27.2.3).
Comparison to Other Methods. In general, the most severe disadvantage is the requirement for ordered grids. FDM usually does not cope with complex or multiscale geometries very well. In such geometries, FDM is computationally very expensive. This is one of the main advantages of FEM and FVM that have significantly less strict requirements on grid structuring. Both methods allow locally adapting the grid to suit the local geometry. Returning to our example with the grooved microfluidic channel, the region of the groove could be resolved at significantly higher accuracy while relaxing the resolution throughout the rest of the microfluidic channel.
However, finite difference schemes are essential for both FEM and FVM. The reason for this is that many differential equations involve time-dependency, which gives rise to a time-dependent differential such as the one found on the left-hand side of the Navier-Stokes equation (see Eq. 11.40). Many numerical solvers will use FEM and FVM to solve the space-dependent terms of a differential equation and then use a central finite, forward finite, or backward finite scheme to step in time. Obviously, the main solution of the differential equation is then obtained using a combination of FEM or FVM and FDM.
30.3 FDM in Microsoft Excel
We have already used FDM indirectly in section 25 as all numerical schemes used for solving ODEs are based on the same mathematics. We will expand these concepts further in section 33 where we will use FDM to solve time- and space-transient solutions to the Navier-Stokes equation and in section 34 where we will implement a fully function three-dimensional numerical solver to calculate complex flows in microfluidics.
At this point we will derive solutions to the Poisson equation which we found for pressure-driven Poiseuille flow (see Eq. 15.14). In section 15.4 we derived analytical solutions to this PDE. Here, we will use FDM to solve this equation using a common spreadsheet analysis tool, e.g., Microsoft Excel[1].
30.3.1 Derivation
The PDE we want to solve is given by Eq. 15.14 as
∂2vx∂z2+∂2vx∂y2=1ηdpdx
where we neglect the volume forces and assume the pressure drop along the channel axis dpdx 
to be a constant. Recall that this PDE can be solved to yield the velocity profile vx (y, z) across the channel cross-section.
We now use a central finite difference scheme (see Eq. 4.10) for approximating the second derivatives as
∂2vx∂y2=F(y+1,z)+F(y−1,z)−2F(y,z)h2∂2vx∂z2=F(y,z+1)+F(y,z−1)−2F(y,z)h2
which when reinserted into Eq. 30.1 yields
F(y+1,z)+F(y−1,z)+F(y,z+1)+F(y,z−1)−4F(y,z)h2=1ηdpdxF(y,z)=14(F(y+1,z)+F(y−1,z)+F(y,z+1)+F(y,z−1)−h2ηdpdx)
We have already derived this equation (see Eq. 28.2) in section 28 where we used sequential dense LU-decomposition and SOR to implement this scheme. Here, we will use Microsoft Excel to do the same. If you recall, Eq. 4.10 is a second-order method which means we can expect a decent numerical stability.
30.3.2 Preparing the Spreadsheet
In the following we will design a spreadsheet which we will use to implement Eq. 30.2. You can find a prepared version of this spreadsheet as file FDM_in_MSExcel.xlsx in the supplementary material. In the following, we will recreate this sheet step by step. We start by creating a new Excel file which we name FDM_in_MSExcel.xlsx.
Grid. On the first sheet we reserve a grid of 40 × 40 cells. This is the domain we will be working on. We have this grid start in the cell C3. We surround this field on all sides with cells in which we enter the value 0. These are the boundary conditions which we assume to be no-slip conditions. In order to illustrate that these are boundary values, we color the cells gray. In addition we add values to indicate the lateral coordinates for y and z. You can see the domain grid in Fig. 30.1a.
Scale Bar. Next to the grid we make a scale bar. For this we first determine the maximum and the minimum values of the grid. For this we use the Excel functions max and min on the whole grid. We then design the scale bar by simply creating a total of 10 cells with increasing equidistant values from the found minimum to the maximum value according to
F(i)=F(i−1)+max−min9
where we set the first value
F(0)=min
In Fig. 30.1a you can see that the value 5 was inserted at an arbitrary location of the grid. The scale bar now shows that values 0.00, 0.56, 1.11, etc. as discreet intervals. You can also see that the values are copied in a second column. We will use this column for coloring.
Adding Color-Coding.Microsoft Excel has a very nice feature that allows representing data in the intuitive “color-coded” format which we are familiar with from numerical software packages. For this, highlight the cells for which the data is to be color-coded and select conditional formatting–>color scales. If you select the second entry red-yellow-green you will receive a color-code which will put the highest value to red and the lowest values to green. This gives the characteristic heat-map-look. If we do this for the second column of our color bar you can see that the scale bar starts to look like a real scale bar is supposed to look like.
Removing the Cells’ Values. In order to improve the look of the color-coded output you can remove the cells’ values by customizing the formatting. For this right-click and select format cells–>number–>custom and enter ;;; as format string. This will remove the cell values. If you look at Fig. 30.1b you can see the output.
Removing the Cells’ Values and Adding Color-Coding to the Grid. We also apply these modifications to the cells of the grid. As a consequence, the lone value of 5 will be converted in a single red block.
Variables. You can also see that the spreadsheet has a region where we can input variables. These variables are the pressure drop dpdx in mbar mm−1, the stepwidth in µ m, the viscosity in mPa s, and the values of the boundary at the top, bottom, left, and right edge of the grid. These values are expected to be given in mm s−1. In order to make these variables usable, we link them to the values of the boundary of our grid.
Working Variables. In this block, we precalculate the value Delta which simply is
Delta=0.1h2ηdpdx
which accounts for the inhomogeneous term of the right-hand side of the Poisson equation (see Eq. 30.1). This value is calculated based on the cells’ values and changes dynamically. The prefactor 0.1 is included to compensate for the units. The unit of Delta is mm s−1. If we insert the values for water, a stepwidth of 2.5 µ m, and a pressure drop of − 0.1 mbar mm−1 we obtain Delta = −0.0625 mm s−1.
Evaluated. In the last block, we evaluate the numerical output of the grid and calculate values we are interested in. The most obvious value is the total flow rate which is calculated according to Eq. 10.9. As we are working with discreet solutions, the integral can be written as
Q=h2ymax∑y=1zmax∑z=1F(y,z)
which yields the unit µl s−1. By multiplying this value with 60 we obtain the unit µl/min.
30.3.3 Implementing the Numerical Scheme
In the next step, we now implement Eq. 30.2. We use the working variable Delta to speed up the calculation slightly. If you look at Fig. 30.2a, you can see how the formula is supposed to look like for the first cell of the grid. If you type in the formula be sure to fix the reference to Delta. Microsoft Excel will most likely issue an error stating that you have just created a circular reference. We ignore this error for a moment. Now simply copy the value of the first cell across the whole grid. The result should look like Fig. 30.2b. Obviously, this does not look like a good solution.
Circular References. The problem is that Eq. 30.2 is an iterative approach. In terms of spreadsheet analysis tools (which Microsoft Excel obviously is), iterations can be translated to circular references. As an example, assume that the value in a spreadsheet in Microsoft Excel is dependent on a second value and v.v.. Let the first value be located in the cell A1, the second value in the cell A2. The formula for the value in A1 depends on A2 which in turn depends on A1. Usually if you try to implement such a calculation in Microsoft Excel, you would receive an error message indicating that you have just created a circular reference. Microsoft Excel would require to calculate A1 then update the value of A2 using the new value of A1, which in turn would require to reevaluate A1 and so forth. The number of calculations would be endless if a precise result is required. Due to this circular reference, Microsoft Excel will usually only calculate the value of cells containing circular references once.
Iterative Calculations. However, it is possible to activate iterative calculations and to tell Microsoft Excel to only perform a finite number of iterations. The circular reference would then still be there, but the program would know that it is only required to perform a finite number of iterations before the value can be accepted as being correct. This is done by selecting File–>Options–>Formulas and checking Enable iterative calculation. We can set the total number of required iterations to a reasonable value where 1000 is usually sufficient. We can also apply a second condition, which is the minimum change between two subsequent iterations. If the value of a cell changes for a value smaller than this threshold, Microsoft Excel will also abort the calculation and accept the result as being exact. This limit obviously depends on the required resolution, i.e., the number of decimal places required for the dependent variable. Usually a reasonable value is somewhere around 10−5.
If we now execute a spreadsheet which has circular references Microsoft Excel will perform the given number of maximum iterations for each cell.
Velocity Profiles. Once you activate iterative calculations, the result obtained should look like Fig. 30.3a. As you can see, the profile obtained resemble the profiles we have seen before for solutions to the Poisson equation. The result obtained with Microsoft Excel is reasonably exact. Fig. 30.3b shows the absolute error between the correct analytical solution given by Eq. 16.62 and the numerical result shown in Fig. 30.3a. As you can see, the error is largest in areas with steep gradients, i.e., near the edges and the corners of the profile. The result can be significantly improved if the grid is further refined, i.e., if smaller stepwidths are used.
Flow Rate. Besides the actual flow profile, the spreadsheet of Fig. 30.3a also outputs the flow rate which is given analytically by Eq. 16.65. The correct result for the chosen values would be Q = 2.11 µl/ min. Our numerical approximation is off slightly, yielding Q = 2.32 µl/ min which is an absolute error of around 10 %. The reason for this is that the flow rate, being the integral of the velocity profile, sums up the total errors. In order to achieve better precision, we would need to use smaller finite changes, i.e., a more refined grid.
30.3.4 Boundary Conditions
We have seen that Microsoft Excel calculates the flow profiles readily. It is straightforward to change the values of the calculation, i.e., use different pressure drops or fluid properties and see how the flow profiles change. The neat setup allows for a number of other interesting experiments.
Nonzero Dirichlet Boundary Conditions. As a first example, we will change the boundary conditions such that they are nonzero at one of the boundaries. If you recall our analytical solution, the sine or cosine series expansion resulted from an eigenfunction expansion in order to make the boundary conditions zero periodically (see Eq. 16.59). Accounting for a nonzero boundary condition at one of the two boundaries along one axis is a somewhat more challenging case analytically. However, for our spreadsheet, this poses no significant problem. As an example, set the top boundary condition to 10 mm s−1. This can be done by simply changing the value in cell AV25. After a few seconds, Microsoft Excel will output the profile shown in Fig. 30.4. As you can see, the velocity maximum moves out of the center and we obtain a somewhat stretched central flow.
Neumann Boundary Conditions. As a second example, we consider the case of Neumann boundary conditions where we do not indicate the value at a given cell of the grid, but the gradient. This is very easy to implement as we simply have to calculate the value of the respective cell as a function of the neighboring cell. Fig. 30.5 shows the example of Couette flow where we have used this boundary condition on the left- and right-hand side of the grid in order to indicate that we are not expecting any changes of the profile in this direction. As you can see, the grid is reused, simply removing the Dirichlet boundary conditions on the left right and right edges of the grid. However, we keep the top and the bottom boundary conditions, which are Dirichlet boundary conditions. We set the speed at which the top boundary is moved at 5 mm s−1. As you can see, we obtain the linear velocity profile expected from Couette flow.
Please note that in this case, we do not apply a pressure drop, as the flow in Couette flow is shear force driven (see section 15.2). The Navier-Stokes equation then becomes a Laplace equation because the right-hand side term is zero (see Eq. 30.1).
30.3.5 Different Profiles
Hagen-Poiseuille Flow Profile. Obviously Microsoft Excel gives us a large degree of freedom to modify the shape of the grid, i.e., to change the channel cross-section. As an example Fig. 30.6a shows the characteristic parabolic flow profile obtained for Hagen-Poiseuille flow in a circular tube (see section 16.3). This example uses a pressure drop of − 0.1 mbar mm−1, a stepwidth of 2.5 µ m, and the fluid properties of water. The diameter of the profile is around 47.5 µ m. As you can see, the profile seems “pixelated.” This is due to the fact that FDM is not very good at adapting to curved surfaces. However, you can see that the parabolic profile is obtained.
From Eq. 16.19 we know the analytical solution which states maximum center velocity of 5.64 mm s−1 according to Eq. 16.20.Microsoft Excel shows a maximum center velocity of around 5.84 mm s−1, which is only slightly off. Again, in order to obtain results with higher accuracy, we would need to increase the number of cells thus also decreasing the error introduced by the pixelation.
Channels With Complex Cross-Sections.Fig. 30.6 illustrates that the spreadsheet we just designed allows deriving the velocity profiles of arbitrarily shaped microfluidic channel cross-sections. Besides the classical Hagen-Poiseuille flow, we could also derive the profile in channels with semicircular-, cross- or triangular-shaped cross-sections. Such cross-sections are actually pretty common. Triangular-shaped channels are obtained, e.g., from laser micromachining [2]. Semicircular cross-sections are preferred for the manufacturing of microfluidic membrane valves and can be created, e.g., by tempering thermoplastics resists which causes reflow and curving of an originally rectangular channel profile [3].
For some of these profiles, analytical solutions exist while others are very challenging or even impossible to derive analytically. In comparison to this, deriving the flow profiles using the spreadsheet is a matter of a few minutes. The most challenging aspect is usually the “painting,”, i.e., the generation of the cross-section. Once this is in place, the values can be changed at will and the velocity profiles are calculated within a few seconds.
30.3.6 Alternative Use of the Spreadsheet
In this section we have used a Microsoft Excel spreadsheet for implementing a FDM in a convenient way. We did this using the Navier-Stokes equation for two-dimensional laminar flow given by Eq. 30.1. Using the example of the Couette flow we have also seen that this spreadsheet can (naturally) also be used if the Poisson equation becomes a Laplace equation because of a zero right-hand side. Interestingly, many phenomena in fluid mechanics, physics and electrical engineering can be described by Laplace and Poisson equations. As an example, take the following equation which is commonly referred to as Gauss’s law (of electrostatics) which relates the electrical field →E 
to the charge density ρel and is given by
Δ→E=ρel∈0
where ε0 is the vacuum permittivity. This equation is one of the famous Maxwell equation. The left-hand side of the equation describes the change in the electric field at all points in the domain. The right-hand side is a source/sink term and describes the electric charges inside of the domain. As you can see, this equation is a Poisson equation just like the simplified Navier-Stokes equation for the Poiseuille flow. In fact, we only need to modify our spreadsheet slightly in order to be usable to solve this equation. The result is shown in Fig. 30.7. As you can see, the output unit is now mV and the pressure term has been replaced by the charge density ρel in C mm−3. We also require 0. The unit of Delta is now likewise mV and the formula needs to be adapted slightly because (based on the units we use) we do not require a prefactor.
In Fig. 30.7 a capacitor with a planar top plate (held at a constant potential of 0 mV) and a non-planar bottom plate (held at a constant potential of 5 mV) is shown. As you can see, the example uses no charge density within the volume. The left- and right-hand side of the domain are assumed to be Neumann-type boundary condition. As you can see, the Microsoft Excel worksheet shows the resulting electrical.
30.4 Summary
In this section, we have introduced FDM as an extension of the classical mathematics of the Taylor series, which we have used extensively in section 25 when studying methods for solving differential equations numerically. As we have already used finite difference methods indirectly before while working on numerical methods, we focused on the implementation of an FDM using Microsoft Excel. For this we designed a spreadsheet which allowed us to solve the Poisson equation of the simplified Navier-Stokes equation for Poisson flow numerically with ease. We derived the velocity profiles of numerous cross-sections and compared the obtained results to the analytical solution wherever one was available. As we have seen, finite difference methods can be implemented conveniently using Microsoft Excel or comparable spreadsheet analysis tools.
We will conclude our discussion of FDM at this point, but return to these methods in section 33 where we will use them to solve time- and space-transient solutions to the Navier-Stokes equation numerically. We will also use FDM in section 34 where we will implement a complete three-dimensional solver based on FDM.