16.1 Elliptical and Circular Profiles

16.1.1 Analytical Solution: An Approach

Circular profiles are a special case of elliptical profiles, which is why we will treat them together. Obviously, these profiles are very relevant to microfluidics as they represent the case of a circular tube or a squeezed capillary in the case of the elliptical profile.

16.1.1.1 Approach for v_x

We do not have an intuitive approach for v_x, which is why we first default to the ellipses equation that is given by

${\left(\frac{x}{a}\right)}^2+{\left(\frac{y}{b}\right)}^2=1$

This function simply describes the contour of our cross-section, where a and b are the ellipse’s half-axis along the y-axis and the half-axis along the z-axis, respectively (see Fig. 16.1a). Obviously, for a = b = r we obtain the equation for a circle.

As stated, we assume no-slip condition for Poiseuille flows, which is why a good first guess would be to choose a function that would ensure that the velocity is in fact zero on the contour ∂C. In this case our approach would read

$v_x\left(y,z\right)=v_0\left(1-{\left(\frac{y}{a}\right)}^2+{\left(\frac{z}{b}\right)}^2\right)$

  (Eq. 16.1)

This function would at least ensure that for all values (x, y) on the contour, we would obtain a velocity of 0. However, we introduced an unknown variable here, namely v₀. If our approach of v_x is correct we should be able to determine this unknown variable after solving the PDE.

16.1.2 Velocity Profile for Elliptical Cross-Sections

We are now able to describe the velocity profile using our approach and Eq. 16.2. We find

$v_x\left(y,z\right)=-\frac{1}{2\eta }\frac{dp}{dx}\frac{a^2{b}^2}{a^2+b^2}\left(1-{\left(\frac{y}{a}\right)}^2+{\left(\frac{z}{b}\right)}^2\right)$

Eq. 16.3 is the equation for the velocity profile of the ellipse. Please note that this equation suggests that the velocity profile is actually negative. This is not the case: The pressure drop $\frac{dp}{dx}$ is also negative as the flow is driven from areas of higher pressure to areas of lower pressure. In many textbooks you may see the pressure gradient being replaced by a pressure drop over a given length of the channel $-\frac{\Delta p}{l}$, in which case Eq. 16.3 can be rewritten to

$v_x\left(y,z\right)=\frac{\Delta p}{2\eta l}\frac{a^2{b}^2}{a^2+b^2}\left(1-{\left(\frac{y}{a}\right)}^2+{\left(\frac{z}{b}\right)}^2\right)$

16.1.3 Flow Rate for Elliptical Cross-Sections

After having derived the velocity profiles we will now derive the flow rate for the Poiseuille flow in cylindrical coordinates according to Eq. 10.9 as

$\frac{dV}{dt}={\int }_{A}vdA$

As the cross-section is elliptical we will carry out this integral in polar coordinates (see Fig. 16.1b). For this we need to transform Eq. 16.4 into polar coordinates. We will express y and z as function of the cylindrical radius 0 ≤ ρ ≤ 1 and the angle 0 ≤ φ ≤ 2π, resulting in

$\begin{array}{lll}y\left(\rho ,\phi \right)\hfill & =\hfill & \rho a\cos \phi \hfill \\ z\left(\rho ,\phi \right)\hfill & =\hfill & \rho b\sin \phi \hfill \\ {\left(\frac{y}{a}\right)}^2+{\left(\frac{z}{b}\right)}^2\hfill & =\hfill & {\rho }^2{\cos }^2\phi +{\rho }^2{\sin }^2\phi ={\rho }^2\hfill \\ v_x\left(\rho ,\phi \right)\hfill & =\hfill & -\frac{1}{2\eta }\frac{dp}{dx}\frac{a^2{b}^2}{a^2+b^2}\left(1-{\rho }^2\right)\hfill \end{array}$

In order to transform the integral from Cartesian to polar coordinates, we need to find the Jacobian determinant (see section 7.5) given by

$\begin{array}{lll}\det J\hfill & =\hfill & \left|\begin{array}{cc}\frac{\partial y}{\partial \rho }& \frac{\partial y}{\partial \phi }\\ \frac{\partial z}{\partial \rho }& \frac{\partial z}{\partial \rho }\end{array}\right|=\left|\begin{array}{cc}\frac{\partial }{\partial \rho }\left(\rho a\cos \phi \right)& \frac{\partial }{\partial \rho }\left(\rho a\cos \phi \right)\\ \frac{\partial }{\partial \rho }\left(\rho b\sin \phi \right)& \frac{\partial }{\partial \rho }\left(\rho b\sin \phi \right)\end{array}\right|\hfill \\ \hfill & =\hfill & \left|\begin{array}{cc}a\cos \phi & -\rho a\sin \phi \\ b\sin \phi & \rho b\cos \phi \end{array}\right|\hfill \\ \hfill & =\hfill & a\cos \phi \rho b\cos \phi +\rho a\sin \phi b\sin \phi \hfill \\ \hfill & =\hfill & \rho ab\left({\sin }^2\phi +{\cos }^2\phi \right)\hfill \\ \hfill & =\hfill & \rho ab\hfill \end{array}$

Having derived the determinant we can now carry out the integration:

$\begin{array}{lll}\frac{dV}{dt}\hfill & =\hfill & {\displaystyle {\int }_0^1{\displaystyle {\int }_0^{2\pi }v\left(\rho ,\phi \right)\rho abd\phi d\rho }}\hfill \\ \hfill & =\hfill & -\frac{1}{2\eta }\frac{dp}{dx}\frac{a^3{b}^3}{a^2+b^2}{\displaystyle {\int }_0^1\left(1-{\rho }^2\right)\rho d\rho {\displaystyle {\int }_0^{2\pi }d\phi }}\hfill \\ \hfill & =\hfill & -\frac{1}{2\eta }\frac{dp}{dx}\frac{a^3{b}^3}{a^2+b^2}{\left[\frac{{\rho }^2}{2}-\frac{{\rho }^4}{4}\right]}_0^{1}2\pi \hfill \\ Q\hfill & =\hfill & -\frac{\pi }{4\eta }\frac{dp}{dx}\frac{a^3{b}^3}{a^2+b^2}\hfill \end{array}$

For details on how the Jacobian determinant is used for integral transformations refer to section 7.3.

16.1.4 Velocity Profile for Circular Cross-Sections

As stated, by simply replacing a = b = r we obtain the velocity profile for the circle from Eq. 16.3:

$\begin{array}{lll}{v}_x\left(y,z\right)\hfill & =\hfill & -\frac{1}{2\eta }\frac{dp}{dx}\frac{r^2}{2}\left(1-{\left(\frac{y}{r}\right)}^2+{\left(\frac{z}{r}\right)}^2\right)\hfill \\ \hfill & =\hfill & \frac{1}{4\eta }\frac{dp}{dx}\left(r^2-\left(y^2+z^2\right)\right)\hfill \end{array}$

or

$v_x\left(y,z\right)=-\frac{\Delta }{4\eta l}\frac{dp}{dx}\left(r^2-\left(y^2+z^2\right)\right)$

if you prefer the notation with the pressure drop over the length of the channel.

16.1.5 Flow Rate for Circular Cross-Sections

Applying a = b = r to Eq. 16.5 yields the flow rate for cylindrical cross-sections given by

$\begin{array}{lll}\frac{dV}{dt}\hfill & =\hfill & -\frac{\pi }{4\eta }\frac{dp}{dx}\frac{r^2}{2r^2}\hfill \\ Q\hfill & =\hfill & -\frac{\pi r^4}{8\eta }\frac{dp}{dx}\hfill \end{array}$

16.1.6 Visualization

In the following we will visualize the equations we found using the Maple listing shown in listing 16.1.

1   restart:read "Core.txt":
2
3   #profile for the ellipsis
4   velx_ellipse:=(y,z,dpdx,a,b,eta)-> -1/(2*eta)*dpdx*(a^2*b^2)/(a^2+b^2)*(1 - ((y/a)^2 + (z/b)^2));
5
6   #we define a helper function that takes care of the units
7   velx_ellipse2:=(y,z,dpdx,a,b)-> simplify(velx_ellipse(y*Unit(um),z*Unit(um),dpdx*Unit(mbar/mm),a*Unit(um),b*Unit(um),1*Unit(mPa*s))/Unit(mm/s));
8
9   #now plot the profile for an elliptical channel
10     quickplot3d([velx_ellipse2(y,z,-0.001,100,50),velx_ellipse2(y,z,-0.01,100,50),velx_ellipse2(y,z,-0.02,100,50)],y=-100..100,z=0..50,0..2,"y [um]","z [um]",typeset(v[x]^empty," [mm/s]"));
11
12  #now plot the profile for a round channel
13     quickplot3d([velx_ellipse2(y,z,-0.001,100,100),velx_ellipse2(y,z,-0.01,100,100),velx_ellipse2(y,z,-0.02,100,100)],y=-100..100,z=0..100,0..2,"y [um]","z [um]",typeset(v[x]^empty," [mm/s]"));

In line 1 we restart the Maple server and include Core.txt. In line 4 we define the flow profile. This is Eq. 16.3, which is called velx_ellipse in the listing. In order to simplify working with units, we define a second function velx_ellipse2, which passes on the independent variables with suitable units. We choose µ m for the dimensions and mbar mm⁻¹ for the pressure drop. Also, we pass along the viscosity value of water. If you want to describe the flow for different fluids, change the variable eta.

Finally in line 10 we calculate the velocity profiles for an elliptical channel with a = 100 µ m and b = 50 µ m. We use the function quickplot3d for creating the plot. Refer to section 2.3.3 for details about this function. We supply different values for the pressure drop, i.e., − 0.001 mbar mm⁻¹, − 0.01 mbar mm⁻¹, and − 0.02 mbar mm⁻¹. The output of the function is shown in Fig. 16.2a. As you can see, for the last value, the flow velocity in the center of the channel approaches 2 mm s⁻¹.

The second case we study is the flow in a circular tube. Here, we use the same values for the different pressure drops and assume a channel of a = b = r = 50 µ m. For the highest pressure drop, the center flow velocity is around 1.25 mm s⁻¹. The calculated flow velocity profiles are shown in Fig. 16.2b.

16.2.1 Derivation

We again start with Eq. 15.14, reducing the equation to the one-dimensional case for which we find

$\begin{array}{lll}\frac{dp}{dx}\hfill & =\hfill & \eta \frac{{\partial }^2{v}_x}{\partial z^2}\hfill \\ \frac{dp}{dx}\hfill & =\hfill & \eta \frac{d^2{v}_x}{d{z}^2}\hfill \end{array}$

Here we have also neglected the volume forces setting k_x = 0. Eq. 16.9 is a second-order ODE for which we need two boundary conditions, which are given by the no-slip condition at the two channel walls: v_x (h/2) = 0, v_x (−h/2) = 0. We integrate Eq. 16.9 twice finding

$\begin{array}{lll}\frac{d{v}_x}{dz}\hfill & =\hfill & \frac{1}{\eta }\frac{dp}{dx}z+c_1\hfill \\ v_x\left(z\right)\hfill & =\hfill & \frac{1}{2\eta }\frac{dp}{dx}{z}^2+c_{1}z+c_2\hfill \end{array}$

Inserting the boundary conditions results in

$\begin{array}{lll}0\hfill & =\hfill & 0\frac{1}{2\eta }\frac{dp}{dx}\frac{h^2}{4}+c_1\frac{h}{2}+c_2\hfill \\ 0\hfill & =\hfill & 0\frac{1}{2\eta }\frac{dp}{dx}\frac{h^2}{4}-c_1\frac{h}{2}+c_2\hfill \end{array}$

from which can find

$c_2=-\frac{1}{2\eta }\frac{dp}{dx}\frac{h^2}{4}$

$c_1=0$

Inserting Eq. 16.11 and Eq. 16.12 into Eq. 16.10 yields

$\begin{array}{lll}{v}_x\left(z\right)\hfill & =\hfill & \frac{1}{2\eta }\frac{dp}{dx}\left(z^2-\frac{h^2}{4}\right)\hfill \\ v_x\left(z\right)\hfill & =\hfill & \frac{1}{2\eta }\frac{dp}{dx}\left(z^2-{\left(\frac{h}{4}\right)}^2\right)\hfill \\ v_x\left(z\right)\hfill & =\hfill & -\frac{h^2}{8\eta }\frac{dp}{dx}\left(1-{\left(\frac{z}{\frac{h}{2}}\right)}^2\right)\hfill \end{array}$

16.2.2 Velocity Profile

Eq. 16.13 described the velocity profile for the planar infinitesimally extended channel. It is a parabola with the maximum velocity

$v_{x,\max }=-\frac{h^2}{8\eta }\frac{dp}{dx}$

Again, the pressure difference is negative since we expect a pressure drop along the channel: $\frac{dp}{dx}<0$. In this case, the overall velocity becomes positive.

16.2.3 Flow Rate

Having derived the velocity profile we can calculate the flow rate using Eq. 10.9 as

$\begin{array}{lll}Q\hfill & =\hfill & \underset{A}{\int }vdA\hfill \\ \hfill & =\hfill & -\frac{h^2}{8\eta }\frac{dp}{dx}{\displaystyle {\int }_0^h\left(1-{\left(\frac{z}{\frac{h}{2}}\right)}^2\right)\omega dz}\hfill \\ \hfill & =\hfill & -\frac{\omega h^2}{8\eta }\frac{dp}{dx}{\left[\left(z-\frac{4}{3H^2}{z}^3\right)\right]}_{{-}^{\frac{h}{2}}}^{{}^{\frac{h}{2}}}\hfill \\ \hfill & =\hfill & -\frac{\omega h^2}{8\eta }\frac{dp}{dx}\left(\frac{h}{2}-\frac{h}{6}-\left(-\frac{h}{2}-\frac{h}{6}\right)\right)\hfill \\ \hfill & =\hfill & -\frac{\omega h^3}{12\eta }\frac{dp}{dx}\hfill \end{array}$

16.2.4 Shear Stress Profile

The shear stress profile can be obtained from Eq. 9.5 and Eq. 16.13. It is given by

$\begin{array}{lll}\tau \left(z\right)\hfill & =\hfill & \frac{h^2}{8}\frac{dp}{dx}\frac{2z}{\frac{h}{2}}\hfill \\ \hfill & =\hfill & \frac{dp}{dx}z\hfill \end{array}$

$=\frac{dp}{dx}\left|z\right|$

The profile needs to be symmetric although Eq. 16.15 does not suggest it: For negative values of z you will obtain negative shear forces, which is physically not meaningful. This can be accounted for, e.g., by only considering the absolute value (see Eq. 16.16).

16.2.5 Visualization

Fig. 16.4 shows the velocity and the shear force profile in an infinitesimally extended channel with 10 µ m height. As you can see, the velocity profile is a parabola, whereas the shear force profile is linear. The maximum velocity in the channel center is about 0.125 mm s⁻¹. Fig. 16.5 shows the calculated profiles for a channel with 50 µ m height. Here, the maximum flow velocity can reach up to above 1.5 mm s⁻¹.

16.3 Flows in Circular Cross-Sections: Hagen-Poiseuille Flow

16.3.1 Introduction

Obviously, circular cross-sections are a very commonly occurring in technical flow systems. Macroscopic examples include pipes and the tap; microfluidic examples include capillaries and circular tubes. The fluid mechanics of circular tubes or channels is commonly referred to as Hagen-Poiseuille¹flow. As we will see, the flow profile is parabolic in shape and the shear stress is linear (see Fig. 16.6). Calculating the flow in circular cross-sections is not difficult as it exploits the fact that the Navier-Stokes equation in cylindrical coordinates is one-dimensional. There are two approaches deriving the same equation. We will first look at the direct solution of the Navier-Stokes equation in cylindrical coordinates which is the Euler equation.

16.3.2 Derivation Using the Euler Equation

Derivation. The first method for approaching the fluid mechanics of circular tube can be directly derived from the Euler equation (see Eq. 14.8), which we derived for the circular flow tube. As such, a circular tube or pipe is nothing else than a circular flow tube, which makes this approach very intuitive. Again we only consider the fully-developed and stationary case and neglect gravity. In this case the Euler equation simplifies to

$\frac{\tau }{r}+\frac{\partial \tau }{\partial r}=\frac{dp}{dz}$

where we assume a pressure drop along the x-axis, which is independent of r; this allows us to treat $\frac{dp}{dx}$ as a constant. If we reduce the fluids to incompressible Newtonian fluids we can use Eq. 9.5 to obtain a dependency between the shear forces τ and the velocity v. Please note that, as we will be working with cylindrical coordinates, we use z as the coordinate axis along the channel, in which case Eq. 16.17 can be written as

$\begin{array}{ll}\frac{\eta }{r}\frac{dv}{dr}+\eta \frac{d^{2}v}{d{r}^2}\hfill & =\frac{dp}{dz}\hfill \\ \frac{1}{r}\frac{dv}{dr}+\frac{d^{2}v}{d{r}^2}\hfill & =\frac{1}{\eta }\frac{dp}{dz}\hfill \end{array}$

Analytical Solution. Eq. 16.18 is a second-order inhomogeneous ODE for which we need to find a suitable solution. Listing 16.2 shows the Maple code we can use to obtain the required solution.

1    restart:read "Core.txt":
2    #we define the differential equation and the boundary values
3    ode1:=1/r*diff(v(r),r)+diff(diff(v(r),r),r) = 1/eta*dpdz;
4    bc1:=v(R)=0;bc2:=D(v)(0)=0;
5
6    #and solve the ODE
7    v:=unapply(rhs(dsolve({ode1,bc1,bc2},v(r))),r,R,eta,dpdz);
8
9    #we define a helper function that will take care of the units
10   v2:=(r,R,dpdz)->simplify(v(r*Unit(um),R*Unit(um),1*Unit(mPa*s),dpdz*Unit(mbar/mm))/Unit(mm/s)):
11
12   #we plot for a 10 um diameter tube and varying pressure drops
13   quickplot([r->v2(r,5,-0.01),r->v2(r,5,-0.05),r->v2(r,5,-0.1)],-5..5,"radius [um]","velocity [mm/s]");
14   #we plot for a 50 um diameter tube and varying pressure drops
15   quickplot([r->v2(r,25,-0.01),r->v2(r,25,-0.05),r->v2(r,25,-0.1)],-25..25,"radius [um]","velocity [mm/s]");
16
17   #calculate the shear forces
18   tau:=unapply(eta*diff(v(r,R,eta,dpdz),r),r,dpdx);
19
20   #again we define a little helper function
21   tau2:=(r,dpdz)->simplify(abs(tau(r*Unit(um),dpdz*Unit(mbar/mm))/Unit(N/m^2))):
22
23   #we plot for a 10 um diameter tube and varying pressure drops
24   quickplot([r->tau2(r,-0.01),r->tau2(r,-0.05),r->tau2(r,-0.1)],-5..5,"radius [um]",typeset("shear stress ",N/m^2));
25
26   #we plot for a 50 um diameter tube and varying pressure drops
27   quickplot([r->tau2(r,-0.01),r->tau2(r,-0.05),r->tau2(r,-0.1)],-25..25,"radius [um]",typeset("shear stress ",N/m^2));

The code is not very complicated. In line 1 we restart the Maple server and include Core.txt. In line 3 we define the ODE and the two boundary conditions. The first boundary condition is given by the fact that the velocity at the radius r = R must be zero (no-slip boundary condition). Furthermore, we expect a rotationally symmetric flow profile in which case the first derivative of the velocity profile at the center (r = 0) must also be zero. Please note that we used the variable dpdz to describe the pressure drop, which will make it easier to work with the resulting function that takes dpdz as an independent variable.

Finally, in line 7 we instruct Maple to solve the ODE. The function unapply allows us to directly assign the right-hand side of our resulting equation to the function v, which uses the independent variables r, the tube radius R, the fluid’s viscosity η, and the pressure drop $\frac{dp}{dx}$. Maple returns the following function:

$\begin{array}{ll}v\left(r\right)\hfill & =\frac{1}{4\eta }\frac{dp}{dz}\left(r^2-R^2\right)\hfill \\ \hfill & =\frac{R^2}{4\eta }\frac{dp}{dz}\left({\left(\frac{r}{R}\right)}^2-1\right)\hfill \end{array}$

$=v_{\max }\left({\left(\frac{r}{R}\right)}^2-1\right)$

Velocity Profile. Eq. 16.20 is the Hagen-Poiseuille flow profile for the rotationally symmetric tube or pipe. It is a parabolic profile with the maximum velocity at the center of the tube. In order to visualize this equation, we define a helper function v2 in line 10. This new function will take the r and R in mm and $\frac{dp}{dz}$ in mbar/mm⁻¹. The viscosity value is set for the value of water. Finally, v2 is divided by mm s⁻¹ and a simplification is applied to it. This results in the function outputting numerical values, which is a prerequisite to plotting functions in Maple.

Starting from line 13 the flow profile of Eq. 16.20 is calculated for two tube diameters using different values for $\frac{dp}{dz}$: − 0.01 mbar/mm⁻¹, − 0.05 mbar/mm⁻¹, and − 0.1 mbar/mm⁻¹. Obviously all of these values must be negative since we assume the pressure to decrease as we move along z due to the losses caused by the fluid’s internal friction. Fig. 16.7 shows the two sets of profiles calculated for circular tubes with diameters 10 µ m and 50 µ m, respectively. As you can see, the flow velocities can become significant for a higher tube diameter R and increasing values of $\frac{dp}{dz}$.

Flow Rate. Having derived the velocity profile we can now derive the flow rate according to

$Q=\underset{A}{\int }vdA$

$\begin{array}{l}=\frac{R^2}{4\eta }\frac{dp}{dz}{\int }_0^R\left({\left(\frac{r}{R}\right)}^2-1\right)2\pi rdr\hfill \\ =\frac{\pi R^2}{2\eta }\frac{dp}{dz}\pi {\int }_0^R\left(\frac{r^2}{R^3}-r\right)dr\hfill \\ =\frac{\pi R^2}{2\eta }\frac{dp}{dz}\pi {\left[\frac{r^4}{4R^2}-\frac{r^2}{2}\right]}_0^R\hfill \\ =-\frac{\pi R^4}{8\eta }\frac{dp}{dz}\hfill \end{array}$

where we have used the infinitesimal area of the circle dA = 2 π r dr (see Eq. 3.104). Eq. 16.22 is sometimes referred to as the Hagen-Poiseuille law. It is a very important function for calculating the flow in technical system, such as pipe and hydraulic systems. From all the equations commonly encountered in microfluidics, Eq. 16.22 may well be the most widely known. Please note that we have (indirectly) already derived this equation when discussing the flow rate through circular cross-section. Thus Eq. 16.8 and Eq. 16.22 are identical.

Shear Stress Profile. Next we will calculate the shear force profile. The shear stress is related to the velocity profile by Eq. 9.5 and can therefore be calculated from Eq. 16.19 as

$\begin{array}{lll}\tau \hfill & =\hfill & \eta \frac{R^2}{4\eta }\frac{dp}{dz}\frac{d}{dr}\left({\left(\frac{r}{R}\right)}^2-1\right)\hfill \\ \hfill & =\hfill & \frac{1}{2}\frac{dp}{dz}r\hfill \end{array}$

Compared to the velocity profile, this function is not symmetrical although it must be in order to yield physically reasonable results. The calculation is correct as the radius cannot have negative values in reality. However, when plotting the profiles, we may want to plot the parabola shape symmetrically in which case we insert negative values for the radius. In listing 16.2, line 18 we let Maple find this derivative as tau and define again a helper function tau2 in line 21. This function ensures that the shear force is positive using the function abs, which returns the absolute value of a function. Furthermore tau2 also ensures that the output of the function is in the correct units, i.e., in N m⁻². Starting from line 24 we plot this function using the two tube diameters used earlier. These plots are shown in Fig. 16.8.

16.3.3 Derivation Using a Force Balance Approach

The second potential approach to finding the Hagen-Poiseuille flow profile is by applying a force balance to a volume segment of the fluid. As we will see, this method is appealing because we do not have to solve a second-order but only a first-order ODE, which is very easy to solve.

Derivation. The derivation follows the schematic shown in Fig. 16.9. The forces acting on a cylindrical infinitesimal volume of length dx and radius dr are the pressure forces acting on the circular surfaces (the difference in which originates from the pressure drop) and the shear stress forces that represent the fluid friction acting on the surface of the volume. The force balance is given by

$\left(p_1-p_2\right)A_{\text{front}}-\tau d{A}_{\text{cylinder}}\underset{\_}{\underset{\_}{!}}0$

where

$p_2=p_1+\frac{\partial p}{\partial z}dz$

$A_{\text{front}}=\pi r^2$

$d{A}_{\text{cylinder}}=2\pi rdz$

$\tau =\eta \frac{\partial v}{\partial r}$

where we have used Eq. 9.5 to convert the shear stress into a function of the velocity and thus restrict the approach to Newtonian fluids. Inserting Eq. 16.25, Eq. 16.26, Eq. 16.27, and Eq. 16.28 into Eq. 16.24 yields

$\begin{array}{lll}-\frac{\partial p}{\partial z}\pi r^{2}dz-\eta \frac{\partial v}{\partial r}2\pi rdz\hfill & \underset{\_}{\underset{\_}{!}}\hfill & 0\hfill \\ \frac{\partial p}{\partial z}r\hfill & =\hfill & -2\eta \frac{\partial v}{\partial r}\hfill \end{array}$

Analytical Solution. Eq. 16.29 is in fact not a PDE but a simple first-order ODE that can be solved by separation of variables to result in

$\begin{array}{lll}-2\eta dv\hfill & =\hfill & \frac{dp}{dz}rdr\hfill \\ -2\eta v\hfill & =\hfill & \frac{dp}{dz}\frac{1}{2}{r}^2+c\hfill \end{array}$

where we can derive the integration constant from the no-slip condition at the wall, i.e., v (R) = 0 from which we find

$\begin{array}{lll}v\left(r\right)\hfill & =\hfill & -\frac{1}{4\eta }\frac{dp}{dz}\left(r^2-R^2\right)\hfill \\ \hfill & =\hfill & -\frac{R^2}{4\eta }\frac{dp}{dz}\left({\left(\frac{r}{R}\right)}^2-1\right)\hfill \end{array}$

As we can see Eq. 16.31 and Eq. 16.20 are identical. We have thereby succeeded in deriving the Hagen-Poiseuille flow profile using this alternative approach.

16.3.4 Average Flow Velocity

As we have seen, the Hagen-Poiseuille flow profile takes the form of a parabola (see Eq. 16.20) with the maximum flow rate in the center of the tube. In many applications it may be advantageous not to work with this velocity profile but with an average velocity v_av that could, e.g., be simply multiplied with the cross-section area A = πR² to obtain the flowrate Q. This average velocity is encountered fairly often because it makes the calculation significantly simpler. We will therefore shortly derive this velocity. The simplest derivation stems from the fact that the flowrate through the cross-section must be identical:

$Q_{\text{av}}=\pi R^2{v}_{\text{av}}$

$Q_{\text{effect}.}={\int }_0^{R}v\left(r\right)2\pi rdr$

where we have used the infinitesimal area of the circle (see Eq. 3.104) for calculating the effective flowrate using the velocity profile. Please note that this is the approach we also took when deriving the flowrate in section 16.3.2. Obviously Eq. 16.32 and Eq. 16.33 must be equal in which case we find

$\begin{array}{lll}{v}_{av}\hfill & =\hfill & \frac{2}{R^2}{\int }_0^{R}v\left(r\right)rdr\hfill \\ \hfill & =\hfill & \frac{1}{2\eta }\frac{dp}{dx}{\displaystyle {\int }_0^R\frac{r^3}{R^2}-rdr}\hfill \\ \hfill & =\hfill & \frac{1}{2\eta }\frac{dp}{dx}{\left[\frac{r^4}{4R^2}-\frac{r^2}{2}\right]}_0^R\hfill \\ \hfill & =\hfill & \frac{1}{2\eta }\frac{dp}{dx}\left(\frac{R^2}{4}-\frac{R^2}{2}\right)\hfill \\ \hfill & =\hfill & -\frac{R^2}{8\eta }\frac{dp}{dx}\hfill \end{array}$

where we note that

$v_{\text{av}}=\frac{v_{\max }}{2}$

which allows us to write the velocity profile also as

$\begin{array}{lll}v\left(r\right)\hfill & =\hfill & -2v_{\text{av}}\left({\left(\frac{r}{R}\right)}^2-1\right)\hfill \\ \hfill & =\hfill & 2v_{\text{av}}\left(1-{\left(\frac{r}{R}\right)}^2\right)\hfill \end{array}$

We will use this notation on several occasions as it is slightly more compact.

16.3.5 Darcy Friction Factor

Darcy Friction Factor for Laminar Hagen-Poiseuille Flow. You may occasionally come across a notation of Eq. 16.34 that can be derived by rearrangement and is given as

$\begin{array}{lll}\frac{dp}{dz}\hfill & =\hfill & -v_{av}\frac{8\eta }{R^2}\frac{\mathrm{Re}}{\mathrm{Re}}\hfill \\ \hfill & =\hfill & -v_{av}\frac{32\eta }{d^2}\frac{2\mathrm{Re}}{2\mathrm{Re}}\hfill \\ \hfill & =\hfill & -v_{av}\frac{64\eta }{d^2}\frac{\mathrm{Re}}{2\mathrm{Re}}\hfill \\ \hfill & =\hfill & -v_{av}f\frac{\eta }{2R^2}\mathrm{Re}\hfill \end{array}$

where we have introduced a factor that is referred to as the Darcy friction factor f for Hagen-Poiseuille flow. It is given by

$f_{\text{Hagen-Poiseuille}}=\frac{64}{\mathrm{Re}}$

The Darcy friction factors are experimentally determined coefficients that allow a correlation between average flow velocity (and therefore volume flowrate) and pressure drop expressed solely by constants and the Reynolds number. This you can see from Eq. 16.37. Eq. 16.38 is the friction factor for laminar Hagen-Poiseuille flow. It is only valid for Reynolds numbers in laminar flow. If the Reynolds number increases to the transition regime or even to turbulent flow Eq. 16.37 can still be applied, but a more complicated equation for the Darcy friction factor would need to be found. As stated, these are usually derived empirically from experimental data. By using only complicated equations for this one factor the overall equations remain very simple. Eq. 16.37 is generally referred to as the Darcy-Weisbach^1,2equation.

Higher Reynolds Numbers: Colebrook-White Equation. If Reynolds numbers are above 4000 different approximations for the Darcy friction factor have to be found. Most common is the following implicit formula that was first given by Colebrook¹[5] based on earlier work by him and White [6], which is why it is commonly referred to as the Colebrook-White equation:

$\frac{1}{\sqrt{f}}=-2{\log }_{10}\left(\frac{\epsilon }{3.7d_H}+\frac{2.51}{\mathrm{Re}\sqrt{f}}\right)$

which uses the roughness height ε and the hydraulic diameter (see section 16.3.6). Eq. 16.39 is an implicit formula, meaning that it needs to be solved numerically as there is no analytical solution for f. However, this equation can be approximated conveniently using algebra tools or even numerical calculators. However, several attempts to finding explicit analytical approximations to this equation for certain fluid mechanical cases have been made, which resulted in a wide range of approximations for f from Eq. 16.39. This work is still ongoing with recent contributions [7].

16.3.6 Hydraulic Diameter

A commonly used method for more complicated geometries is the usage of the hydraulic diameter. This virtual diameter is defined as

$d_H=\frac{4A}{l_c}=\frac{4\cdot \text{cross-section area}}{\text{circumference}}$

It essentially “converts” an arbitrary cross-section to a diameter of an (fluid mechanically more-or-less) equivalent circular cross-section. Usually, the equations derived for the Hagen-Poiseuille flow are then simply applied. However, for most profiles, the errors one makes when using the hydraulic diameter instead of the correct equations for the respective cross-sections are quite severe.

16.4 Flows in Rectangular Cross-Sections: Solution to Poisson and Laplace Equations

16.4.1 Introduction

Until now, we have considered mainly cross-sections that can be reduced to one-dimensional problems with the exception of the elliptical cross-section, where we have used a trial-and-error approach that proved to be successful. Now we will turn to the case of a rectangular channel cross-section. As we will see, this cross-section is significantly more complicated. Again, we begin with Eq. 15.14 from which we find

$\frac{dp}{dx}=\eta \left(\frac{{\partial }^2{v}_x}{\partial z^2}+\frac{{\partial }^2{v}_x}{\partial y^2}\right)+k_x$

In the first step, we will neglect volume forces which reduces the equation to

$\frac{{\partial }^2{v}_x}{\partial z^2}+\frac{{\partial }^2{v}_x}{\partial y^2}=\frac{1}{\eta }\frac{dp}{dx}$

The right-hand side of Eq. 16.40 is a constant. Before attempting to solve this equation, we will spend a moment to consider what we actually learn from it. It tells us that for any point in the cross-section the change of tangents to the velocity profile in the direction of the y-axis plus the change of the tangents to the velocity profile in the direction of the z-axis will always sum up to a constant that is dependent on the viscosity and the pressure drop driving the flow. Let us refer to these “change of the tangents” as the “change in steepness of the velocity profile”. The lower the viscosity and/or the higher the pressure drop, the bigger this sum of changes will become. In other words, with decreasing viscosity and increasing pressure drop, the flow profile will change steepness faster, i.e., the relative changes in the steepness of the velocity profiles will become bigger.

Obviously these observations do not help us in finding a solution directly. However, understanding what the equation tells us makes it somewhat more intuitive to understand the problem. We need to find an equation for which, at any given point (y, z) in the cross-section the sum of the changes along the y-axis and z-axis are the given constant. Therefore, if we expect a large change along the y-axis, the change along the z-axis must be small and v.v. We have already calculated the flow profile in an infinitesimally elongated channel to be of parabolic shape (see section 16.2). In this case the change along the y-axis was neglected.

16.4.2 Naive Approach

We can approach this problem naively and assume that Maple may be able to solve this problem for us. The following line will ask Maple to solve Eq. 16.40 for us using the command pdsolve, which is used to solve PDEs:

equ1:=diff(diff(u[x](y,z),y),y)+diff(diff(u[x](y,z),z),z) = -1/eta*dpdx;
pdsolve(equ1,u[x](y,z));

The output we will obtain will most likely look somewhat like this:

u[x](y, z) = _F1(z - I y) + _F2(z + I y) - (dxdy*y^2)/(eta)

What Maple is telling is that it expects the solution to consist of two functions that take complex argument expressions involving y and z. It also expects this function to contain a term only containing y, which mainly takes care of the constant. Obviously, this solution is not very helpful. Besides being rather vague, this solution also does not account for our boundary conditions. All together, this was not a very successful approach.

16.4.3 General Approach: Substitution and Separation of Variables

We have derived a number of solutions to commonly encountered (partial) differential equations in section 8.2 and section 8.3. We will now apply these methods for finding the solution to Eq. 16.40. In the more general terms, we will solve the Poisson equation (see section 8.1.5), which is (in general) given by Eq. 8.17 as

$\Delta f\left(y,z\right)=g\left(y,z\right)$

The Poisson equation is an inhomogeneous partial differential equation. We therefore first default to solving the homogeneous case, i.e., the Laplace equation (see section 8.1.6) given by Eq. 8.18 as

$\Delta f\left(y,z\right)=0$

We will solve this function in a two-dimensional case on a rectangular boundary 0 ≤ y ≤ w and 0 ≤ z ≤ h with the boundary conditions f (0, z) = f (w, z) = f (y, 0) = f (y, h) = 0. In fluid mechanical terms, this solution will give us the Poiseuille flow profile in channels with rectangular cross-sections of height h and width w. We assume no-slip condition which gives rise to our boundary conditions.

We proceed using a separation of variables approach using f (y, z) = Y (y) Z (z) as detailed in section 8.3.3.5

$\begin{array}{lll}f\left(y,z\right)\hfill & =\hfill & Y\left(y\right)Z\left(z\right)\hfill \\ \Delta f\left(y,z\right)\hfill & =\hfill & \frac{{\partial }^{2}Y}{\partial y^2}Y+\frac{{\partial }^{2}Z}{\partial z^2}Z=0\hfill \end{array}$

There are two main ways of approaching Eq. 16.43, both of which you are likely to encounter in fluid mechanics textbooks. It is a question of personal favor which one to choose. We will briefly outline both.

16.4.4 Single λ Approach

The first method of going about it could be called a “single λ approach”. In this case we rewrite Eq. 16.43 to

$\begin{array}{lll}\frac{{\partial }^{2}Y}{\partial y^2}Z\hfill & =\hfill & -\frac{{\partial }^{2}Z}{\partial z^2}Y\hfill \\ -\frac{1}{Y}\frac{{\partial }^{2}Y}{\partial y^2}\hfill & =\hfill & \frac{1}{Z}\frac{{\partial }^{2}Z}{\partial z^2}=-{\text{λ}}^2\hfill \end{array}$

where we again state that in order for this equation to be fulfilled, both derivatives must equal a constant, which we term −λ². We have taken this approach in section 8.3.3.5 while solving the one-dimensional heat equation. We use λ as a squared term here in order to be exact about the signs.

Coordinate System and Boundary Conditions. Depending on the textbook you may find different versions of the solution to Eq. 16.42. Fig. 16.10 lists the most commonly encountered coordinate systems used.

Solution on the Bounds$-\frac{w}{2}\le y\le +\frac{w}{2}$ and $0\le \text{z }\le \text{h}$. This case is depicted by the coordinate system in Fig. 16.10a. We use the following Dirichlet boundary conditions:

$\begin{array}{lll}f\left(-\frac{w}{2},z\right)\hfill & =\hfill & 0\to Y\left(-\frac{w}{2}\right)=0\hfill \\ f\left(+\frac{w}{2},z\right)\hfill & =\hfill & 0\to Y\left(\frac{w}{2}\right)=0\hfill \\ f\left(y,0\right)\hfill & =\hfill & 0\to Z\left(0\right)=0\hfill \\ f\left(y,h\right)\hfill & =\hfill & 0\to Z\left(h\right)=0\hfill \end{array}$

We have already done most of the work in finding the solution to this case. Referring to Tab. 8.2 we find the following eigenvalue functions as a solution for Z (Z) with λ < 0:

$Z(z)={\displaystyle \sum _{n=1}^{\infty }{c}_n}\sin \left(n\pi \frac{z}{n}\right)$

which also gives us

$\lambda =\frac{n\pi }{h}$

The case for Y (y) is less straightforward. Referring to Tab. 8.2 we find for λ > 0 the solution

Y (y) = c₀ sinh (λ y) + c₁ cosh (λ y)

However, after applying the boundary condition, we are only left with the trivial solution. This case happens surprisingly often when solving differential equations and it is due to one simple fact: We know too much of the solution. Remember, the boundary conditions are used to determine λ. We already applied two boundary conditions and solved for λ. If we apply two boundary conditions again, we will end up with a situation similar to having more equations than unknowns to solve for. In applying more boundary conditions than we actually need, we are left with only the trivial solution as a result.

However, we can also apply a different type of boundary condition. We know that the system we are investigating is symmetric along the z-axis. Therefore, it is reasonable to assume that the solution function will also be symmetric along the z-axis and thus the resulting function will have its maximum at y = 0. In this case, we can assume the first derivative of Y (y) to be zero at the origin: ${\frac{dY}{dy}|}_{y=\frac{w}{0}}=0$. This is a Neumann boundary condition. Again referring to Tab. 8.2 we find the corresponding solution to be

$\begin{array}{lll}Y\left(y\right)\hfill & =\hfill & c_1\cosh \left(\lambda y\right)\hfill \\ \hfill & =\hfill & c_1\cosh \left(n\pi \frac{y}{h}\right)\hfill \end{array}$

where we have already inserted the solution for λ from Eq. 16.46. Now, Y (y) does not yet fulfill our boundary conditions $Y\left(-\frac{w}{2}\right)=Y\left(+\frac{w}{2}\right)=0$. Unfortunately plugging in these boundary conditions will leave us only with the trivial solution for c₁ = 0, and we need to do a little more thinking. If we recall the shape of the hyperbolic cosine function (see Fig. 3.4b) we know that this function is never 0. In order to satisfy our boundary conditions, we need to use the negative hyperbolic sine function and shift it such that it will be zero at $Y\left(\frac{w}{2}\right)\text{and}Y\left(+\frac{w}{2}\right)$. As the function is symmetric, we only need to account for one of the points, which makes life a little easier. We assume our function to look like the following:

$\begin{array}{lll}Y(y)\hfill & =\hfill & c_2-c_1\cosh \left(n\pi \frac{y}{h}\right)\hfill \\ Y\left(\frac{w}{2}\right)\hfill & =\hfill & c_2-c_1\cosh \left(n\pi \frac{w}{2h}\right)\underset{\_}{\underset{\_}{!}}0\to c_2=c_1\cosh \left(n\pi \frac{w}{2h}\right)\hfill \\ Y\left(y\right)\hfill & =\hfill & c_1\left(\cosh \left(n\pi \frac{w}{2h}\right)-\cosh \left(n\pi \frac{y}{h}\right)\right)\hfill \\ \hfill & =\hfill & c_1\cosh \left(n\pi \frac{w}{2h}\right)\left(1-\frac{\cosh \left(n\pi \frac{y}{h}\right)}{\cosh \left(n\pi \frac{w}{2h}\right)}\right)\hfill \\ \hfill & =\hfill & c_3\left(1-\frac{\cosh \left(n\pi \frac{y}{h}\right)}{\cosh \left(n\pi \frac{w}{2h}\right)}\right)\hfill \end{array}$

where the function now satisfies the boundary condition. Having derived Y (y) and Z (z) we can now assemble the homogeneous solution to Eq. 16.42:

$\begin{array}{lll}{f}_h\left(y,z\right)\hfill & =\hfill & Y(y)Z(z)\hfill \\ \hfill & =\hfill & c_1{c}_3\left(1-\frac{\cosh \left(n\pi \frac{y}{h}\right)}{\cosh \left(n\pi \frac{w}{2h}\right)}\right){\displaystyle \sum _{n=1}^{\infty }{c}_n\sin \left(n\pi \frac{z}{h}\right)}\hfill \\ \hfill & =\hfill & {\displaystyle \sum _{n=1}^{\infty }{c}_n\left(1-\frac{\cosh \left(n\pi \frac{y}{h}\right)}{\cosh \left(n\pi \frac{w}{2h}\right)}\right)\sin \left(n\pi \frac{z}{h}\right)}\hfill \end{array}$

This is the solution to the homogeneous Poisson (i.e., Laplace equation) given by Eq. 16.42. Obviously, the constants c_n still need to be determined. We will do this by taking into account the term right-hand side of Eq. 16.41, i.e., g (y, z). This will provide us with the inhomogeneous solution we are seeking.

Inhomogeneous Solution to the Navier-Stokes Equation. We first start by forming the second-order partial differentials of Eq. 16.49 that we need in order to take into account the right-hand side equation g (y, z). In our case we find

$\begin{array}{lll}\frac{{\partial }^2{f}_h}{\partial y^2}\hfill & =\hfill & -{\displaystyle \sum _{n=1}^{\infty }{c}_n{\left(\frac{n\pi }{h}\right)}^2\frac{\cosh \left(n\pi \frac{y}{h}\right)}{\cosh \left(n\pi \frac{w}{2h}\right)}\sin \left(n\pi \frac{z}{h}\right)}\hfill \\ \frac{{\partial }^2{f}_h}{\partial z^2}\hfill & =\hfill & -{\displaystyle \sum _{n=1}^{\infty }{c}_n}\left(1-\frac{\cosh \left(n\pi \frac{y}{h}\right)}{\cosh \left(n\pi \frac{w}{2h}\right)}\right)+{\left(\frac{n\pi }{h}\right)}^2\sin \left(n\pi \frac{z}{h}\right)\hfill \\ \frac{{\partial }^2{f}_h}{\partial y^2}+\frac{{\partial }^2{f}_h}{\partial z^2}\hfill & =\hfill & -{\displaystyle \sum _{n=1}^{\infty }{c}_n{\left(\frac{n\pi }{h}\right)}^2\sin \left(n\pi \frac{z}{h}\right)\left(\frac{\cosh \left(n\pi \frac{y}{h}\right)}{\cosh \left(n\pi \frac{w}{2h}\right)}+1-\frac{\cosh \left(n\pi \frac{y}{h}\right)}{\cosh \left(n\pi \frac{w}{2h}\right)}\right)}\hfill \\ \hfill & =\hfill & {\displaystyle \sum _{n=1}^{\infty }{c}_n{\left(\frac{n\pi }{h}\right)}^2\sin \left(n\pi \frac{z}{h}\right)}\hfill \end{array}$

From Eq. 16.41 we know that this must equal the right-hand side g (y, z). If you look carefully at Eq. 16.50 you will see that the sine series that we already had in our initial result (Eq. 16.49) is still there. It therefore makes perfect sense to expand the right-hand side of Eq. 16.41 to a sine series as described in section 4.3.3.8. We can then determine c_n simply by comparing the coefficients. Referring to Eq. 16.40 you see that the right-hand side is not a function of y and z:

$g\left(y,z\right)=\frac{1}{\eta }\frac{dp}{dx}=\text{const}.$

We therefore merely need to Fourier expand this constant to a sine function as detailed in Eq. 4.34. The expansion should be along the z-axis since this is the axis along which our sine function was expanded (see Eq. 16.49). In this case we find

$\frac{1}{\eta }\frac{dp}{dx}={\displaystyle \sum _{n=0}^{\infty }\frac{4}{\eta \pi }\frac{dp}{dx}\frac{1}{\left(2n+1\right)}\sin \left(\left(2n+1\right)\pi \frac{z}{h}\right)}$

Now we need to match the coefficients of Eq. 16.50 and Eq. 16.51. First we find that we only need to take into account all terms for odd values of n since the even terms vanish when expanding a constant to a Fourier series (see section 4.3.3.3). We therefore find

${\displaystyle \sum _{n=0}^{\infty }\frac{4}{\eta \pi }\frac{dp}{dx}\frac{1}{\left(2n+1\right)}\sin \left(\left(2n+1\right)\pi \frac{z}{h}\right)\underset{\_}{\underset{\_}{!}}-{\displaystyle \sum _{n=0}^{\infty }{c}_n}{\left(\frac{\left(2n+1\right)\pi }{h}\right)}^2\sin \left(\left(2n+1\right)\pi \frac{z}{h}\right)}$

from which we find

$\begin{array}{lll}-\frac{4}{\eta \pi }\frac{dp}{dx}\frac{1}{\left(2n+1\right)}\hfill & \underset{\_}{\underset{\_}{!}}\hfill & c_n{\left(\frac{\left(2n+1\right)\pi }{h}\right)}^2\hfill \\ c_n\hfill & =\hfill & -\frac{4h^2}{\eta {\pi }^3}\frac{dp}{dx}\frac{1}{{\left(2n+1\right)}^3}\hfill \end{array}$

which gives us the missing constants c_n. We are now able to formulate the overall solution to the Poisson equation in rectangular channel cross-sections given by Eq. 16.49 and Eq. 16.52

$\begin{array}{lll}f\left(y,z\right)\hfill & =\hfill & -{\displaystyle \sum _{n=0}^{\infty }\frac{4h^2}{\eta {\pi }^3}\frac{dp}{dx}\frac{1}{{\left(2n+1\right)}^3}\left(1-\frac{\cosh \left(\left(2n+1\right)\pi \frac{y}{h}\right)}{\cosh \left(\left(2n+1\right)\pi \frac{w}{2h}\right)}\right)\sin \left(\left(2n+1\right)\pi \frac{z}{h}\right)}\hfill \\ \hfill & =\hfill & \frac{4h^2}{\eta {\pi }^3}\frac{dp}{dx}{\displaystyle \sum _{n=0}^{\infty }\frac{1}{{\left(2n+1\right)}^3}\left(1-\frac{\cosh \left(\left(2n+1\right)\pi \frac{y}{h}\right)}{\cosh \left(\left(2n+1\right)\pi \frac{w}{2h}\right)}\right)\sin \left(\left(2n+1\right)\pi \frac{z}{h}\right)}\hfill \end{array}$

Solutions on the Bounds$-\frac{w}{2}\le y\le +\frac{w}{2}and-\frac{h}{2}\le z\le +\frac{h}{2}$. You may come across different solutions for the Poiseuille flow profile in rectangular cross-sections if a different coordinate system is used. The case we will now discuss is depicted in Fig. 16.10b. We use the following Dirichlet boundary conditions:

$\begin{array}{l}f\left(-\frac{w}{2},z\right)=0\to Y\left(-\frac{w}{2}\right)=0\hfill \\ f\left(+\frac{w}{2},z\right)=0\to Y\left(\frac{w}{2}\right)=0\hfill \\ f\left(y,-\frac{h}{2}\right)=0\to Z\left(-\frac{h}{2}\right)=0\hfill \\ f\left(y,+\frac{h}{2}\right)=0\to Z\left(+\frac{h}{2}\right)=0\hfill \end{array}$

We will start by finding the solution to Z (z), which we can simply look up in Tab. 8.2 after applying the boundary conditions:

$Z(z)={\displaystyle \sum _{n=0}^{\infty }{c}_n\cos \left(\left(2n+1\right)\frac{\pi }{h}z\right)\to \lambda =\left(2n+1\right)\frac{\pi }{h}}$

Now we can use the same derivation for Y (y) we took for the coordinate system of Fig. 16.10a. Again we find

$\begin{array}{lll}Y\left(y\right)\hfill & =\hfill & c_2-c_1\cosh \left(\left(2n+1\right)\pi \frac{y}{h}\right)\hfill \\ Y\left(\frac{w}{2}\right)\hfill & =\hfill & c_2-c_1\cosh \left(\left(2n+1\right)\pi \frac{w}{2h}\right)\underset{\_}{\underset{\_}{!}}0\to c_2=c_1\cosh \left(\left(2n+1\right)\pi \frac{w}{2h}\right)\hfill \\ Y\left(y\right)\hfill & =\hfill & c_1\left(\cosh \left(2n+1\right)\pi \frac{w}{2h}\right)-\cosh \left(\left(2n+1\right)\pi \frac{y}{h}\right)\hfill \\ \hfill & =\hfill & c_1\cosh \left(\left(2n+1\right)\pi \frac{w}{2h}\right)\left(1-\frac{\cosh \left(\left(2n+1\right)\pi \frac{y}{h}\right)}{\cosh \left(\left(2n+1\right)\pi \frac{w}{2h}\right)}\right)\hfill \\ \hfill & =\hfill & c_3\left(1-\frac{\cosh \left(\left(2n+1\right)\pi \frac{y}{h}\right)}{\cosh \left(\left(2n+1\right)\pi \frac{w}{2h}\right)}\right)\hfill \end{array}$

which is almost the same as Eq. 16.48 but uses a different λ. Therefore we find as the homogeneous solution:

$\begin{array}{lll}{f}_h\left(y,z\right)\hfill & =\hfill & Y\left(y\right)Z\left(z\right)\hfill \\ \hfill & =\hfill & c_1{c}_3\left(1-\frac{\cosh \left(\left(2n+1\right)\pi \frac{y}{h}\right)}{\cosh \left(\left(2n+1\right)\pi \frac{w}{2h}\right)}\right){\displaystyle \sum _{n=0}^{\infty }{c}_n\cos \left(\left(2n+1\right)\pi \frac{z}{h}\right)}\hfill \\ \hfill & =\hfill & {\displaystyle \sum _{n=0}^{\infty }{c}_n\left(1-\frac{\cosh \left(\left(2n+1\right)\pi \frac{y}{h}\right)}{\cosh \left(\left(2n+1\right)\pi \frac{w}{2h}\right)}\right)\cos \left(\left(2n+1\right)\pi \frac{z}{h}\right)}\hfill \end{array}$

Obviously, we now have a cosine series, which is why we expand the right-hand side of Eq. 16.40 to a cosine Fourier series according to Eq. 4.42, in which case we find

$\frac{1}{\eta }\frac{dp}{dx}={\displaystyle \sum _{n=0}^{\infty }{\left(-1\right)}^n\frac{4}{\left(2n+1\right)\pi }\frac{4}{\eta \pi }\frac{dp}{dx}\cos \left(\frac{\left(2n+1\right)\pi z}{h}\right)}$

By comparison with Eq. 16.55 we can now determine the missing coefficient c_n. However, we must first calculate the left-hand side of Eq. 16.55 for which we need the partial derivatives of Eq. 16.55:

$\begin{array}{lll}\frac{{\partial }^2{f}_h}{\partial y^2}\hfill & =\hfill & -{\displaystyle \sum _{n=0}^{\infty }{c}_n}{\left(\frac{\left(2n+1\right)\pi }{h}\right)}^2\frac{\cosh \left(\left(2n+1\right)\pi \frac{y}{h}\right)}{\cosh \left(\left(2n+1\right)\pi \frac{w}{2h}\right)}\cos \left(\left(2+1\right)\pi \frac{z}{h}\right)\hfill \\ \frac{{\partial }^2{f}_h}{\partial z^2}\hfill & =\hfill & -{\displaystyle \sum _{n=0}^{\infty }{c}_n\left(1-\frac{\cosh \left(\left(2n+1\right)\pi \frac{y}{h}\right)}{\cosh \left(\left(2n+1\right)\pi \frac{w}{2h}\right)}\right)+{\left(\frac{\left(2n+1\right)\pi }{h}\right)}^2\cos \left(\left(2n+1\right)\pi \frac{z}{h}\right)}\hfill \\ \frac{{\partial }^2{f}_h}{\partial y^2}+\frac{{\partial }^2{f}_h}{\partial z^2}\hfill & =\hfill & -{\displaystyle \sum _{n=0}^{\infty }{c}_n{\left(\frac{\left(2n+1\right)\pi }{h}\right)}^2\cos \left(\left(2n+1\right)\pi \frac{z}{h}\right)\left(\frac{\cosh \left(\left(2n+1\right)\pi \frac{y}{h}\right)}{\cosh \left(\left(2n+1\right)\pi \frac{w}{2h}\right)}+1-\frac{\cosh \left(\left(2n+1\right)\pi \frac{y}{h}\right)}{\cosh \left(\left(2n+1\right)\pi \frac{w}{2h}\right)}\right)}\hfill \\ \hfill & =\hfill & -{\displaystyle \sum _{n=0}^{\infty }{c}_n{\left(\frac{\left(2n+1\right)\pi }{h}\right)}^2\cos \left(\left(2n+1\right)\pi \frac{z}{h}\right)}\hfill \\ \hfill & \hfill & \hfill \end{array}$