Combined Pressure-Driven and Electroosmotic Slip Flow through Elliptic Cylindrical Microchannels: The Effect of the Eccentricity of the Channel Cross-Section

Abstract

Electroosmotic force has been used extensively to manipulate fluid flow in a microfluidic system with various channel shapes, especially an elliptic cylinder. However, developing a computational domain and simulating fluid flow for a system involving an elliptic channel consumes a large amount of time. Moreover, the mathematical expression for the fluid velocity of electroosmotic flow in an elliptic channel may be given in the form of the Mathieu functions that have difficulty in achieving the numerical result. In addition, there is clear scientific evidence that confirms the slippage of fluid at the solid-fluid interface in a microscale system. In this study, we present the mathematical model of combined pressure-driven and electroosmotic flow through elliptic microchannels under the slip-fluid condition. From the practical point of view in fluidics, the effect of the eccentricity of the channel cross-section is investigated on the volumetric flow rate to overcome the difficulty. The results show that the substitution of the equivalent circular channel for an elliptic channel provides a valid flow rate under the situation that the areas of both channel cross-sections are equal and the eccentricity of the elliptic cross-section is less than 0.5. Additionally, the flow rate obtained from the substitution is more accurate when the slip length increases or the pressure-gradient-to-external-electric-field ratio decreases.

1. Introduction

Microfluidics is the study of fluid manipulation through a small-scale channel (microchannel) involving various fields of science such as engineering, biology, and chemistry. The research interest in microfluidics has been increasing over the past 10 years [1,2,3,4,5,6,7,8] in order to develop and create many microscale tools; for example, organs-on-a-chip [9,10], labs-on-a-chips [11], 3D printing [12,13], and RNA-delivery devices [14]. Some primary advantages of these technologies [15] are: (i) reducing the consumption of sample and reagent quantities, (ii) saving cost and time in operation and manufacturing products, and (iii) performing detection and separation of molecules at high resolution and sensitivity. An important role in driving fluid in a microchannel is played by electroosmosis.

Electroosmosis is the process of the transportation of fluid through a narrow channel under an electric field. This process occurs due to the effect of ionic charges at the fluid and solid interface of the wall channel [16]. A charged surface develops two parallel layers of ions known as the electrical double layers (EDL) as shown in Figure 1. The first layer contains the ions that are absorbed onto the inner wall surface by the chemical reactions. The second layer is the counterions, which are free ions, in the fluid that are attached to the channel wall due to the electrokinetic force from the ions in the first layer. As a result, the greater the distance of the wall surface, the less the counterion concentration. Therefore, the net charge density and the potential distribution decrease with a rise in the distance of the channel wall and reduce to zero at the bulk fluid. Moreover, the ion and potential distribution can be expressed by the Poisson–Boltzmann equation. An applied external electric field will lead to the movement of the ions in the second layer and result in the flow of the fluid. This type of fluid movement is also known as electroosmotic flow (EOF). In order to investigate the EOF’s behavior in microfluidic devices, many researchers have carried out numerous works [17,18,19,20,21,22,23] by both experimental and theoretical methods.

Various types of cross-sectional shapes were considered in the study of EOF, such as circles [24,25], ellipses [26,27], and rectangles [28,29]. One of the most important cross-sections of a microchannel existing in natural systems and microfluidic fabrication is in an elliptic shape. However, developing a computational domain for a system involving an elliptic channel consumes a large amount of time. Moreover, due to a lack of circular symmetry, a simulation in an elliptic channel entails a high computational cost compared to that in a circular channel. Additionally, in a theoretical study using a mathematical model, the velocity of EOF through an elliptic channel may be derived in the form of the special functions called the Mathieu and modified Mathieu functions. These special functions rely on the infinite series of the trigonometric or the hyperbolic function with coefficients. For example, the Mathieu function ${ce}_{2m}(z;q)$ is expressed by

$${ce}_{2m}(z;q)=\sum _{r=0}^{\infty }{A}_{2r}^{\left(2m\right)}cos2rz,$$

(1)

where $A_{2r}^{\left(2m\right)}$ denote the coefficients depending on the parameter $a_q$, which is a function of the parameter q. When the value of q is not quite small, a little error of $a_q$ yields a large error of $A_{2r}^{\left(2m\right)}$, especially when m and r are large. Therefore, the Mathieu and modified Mathieu functions may have difficulty in achieving the numerical result due to their stability issue [30]. To overcome this complication, the validation of using the circular channel instead has been proposed. In 2015, Liu et al. [31] investigated the influence of eccentricity on fully developed pure EOF through elliptic channels. The force by the pressure gradient was omitted. By fixing the perimeter of the channel cross-section, the fluid velocities in elliptic channels having different values of the eccentricity were evaluated. They concluded that considering the circular channel with the same perimeter instead of an elliptic channel is inappropriate. In 2021, Numpanviwat et al. [32] proposed the mathematical model of transient combined pressure-driven and electroosmotic flow in an elliptic microchannel under the no-slip condition. The effect of the eccentricity of the channel cross-section was studied in two situations by fixing either the area or the perimeter of the cross-section. The result is well consistent with Liu’s work when fixing the perimeter. Nonetheless, the result showed that the use of the cylindrical channel by fixing the area of the cross-section is validated compared with the result from the elliptic channel.

In a macroscale system, fluid flow problems are usually considered under the no-slip boundary condition on solid surfaces, which means that the velocity of the fluid at the surface will be zero. However, over the past few decades, there is much evidence indicating that the fluid slips at the wall surfaces [33,34,35,36]. Since the surface-area-to-volume ratio of the microscale system is very large, the slippage of fluid will significantly affect the interaction between the wall surface and the fluid during the fluid transport process. As a consequence, the slip-fluid condition will play a crucial role in fluid flow under a microscale or a smaller system.

Motivated by the aforementioned arguments, we present the mathematical model of combined pressure-driven and electroosmotic flow through elliptic cross-sectional microchannels under the slip fluid condition in order to examine the influence of the eccentricity of the elliptic-channel cross-section. The force by the EDL potential is derived from the Debye-Hückel approximation by employing the Mathieu and modified Mathieu functions. Then, the fluid velocity is obtained by solving the Navier–Stokes equations coupled with the continuity equation. Since a volumetric flow rate plays a more significant role in the microfluidic aspect [37], the effect of the eccentricity on the volumetric flow rate is investigated instead of the velocity. The investigations are performed in two situations: (i) the area of the elliptic cross-section of the channel is fixed and (ii) the perimeter of the elliptic cross-section of the channel is fixed. With regard to the slip flow phenomenon, the comparison of the effect of the eccentricity on various slip lengths is also carried out. To validate the consideration of the equivalent circular microchannel instead of the elliptic microchannel, the relative error of the flow rate with various eccentricities is demonstrated. Further, since the fluid velocity is expressed relating to the pressure gradient and the external electric field, it is important to consider the influence of the eccentricity relating to the ratio of the pressure gradient and the external electric field.

The rest of this paper is organized as follows. Section 2 provides background knowledge of the elliptic coordinate system and the Mathieu and modified Mathieu functions. The mathematical modeling for combined pressure-driven and electroosmotic flow is established in Section 3. Section 4 shows the boundary conditions used in this study. Section 5 presents the mathematical methodology for solving the boundary value problem. The numerical results are illustrated and discussed in Section 6. Finally, the conclusions are given in Section 7.

2. Preliminaries on the Elliptic Coordinate System

2.1. Elliptic Coordinates

According to an ellipse having two foci at $(c,0)$ and $(-c,0)$ in the 2-dimensional Cartesian coordinates, the elliptic coordinate system, as shown in Figure 2, is defined by

$$\begin{array}{c}\hfill x=ccosh\xi cos\eta ,\end{array}$$

(2)

$$\begin{array}{c}\hfill y=csinh\xi sin\eta ,\end{array}$$

(3)

where $\xi \in [0,\infty )$ denotes the confocal ellipses and $\eta \in [0,2\pi )$ denotes the asymptotic angle of the confocal hyperbolas.

The elongation of an ellipse along the major axis called the eccentricity e is measured by the ratio

$$e=\frac{c}{a},0<e<1,$$

(4)

where c denotes the focal length of an ellipse and a denotes the longest distance between the center and perimeter of an ellipse known as the semi-major axis.

The scale factors, also known as the basic vectors, for the elliptic coordinates $(\xi ,\eta )$ can be derived as follows:

$$h_{\xi }=h_{\eta }=c\sqrt{{cosh}^2\xi -{cos}^2\eta }.$$

(5)

As a consequence, the 2-dimensional Laplace operator in the elliptic coordinates $(\xi ,\eta )$ are scaled to

$${\nabla }^2=\frac{1}{h_{\xi }{h}_{\eta }}\left(\frac{{\partial }^2}{\partial {\xi }^2}+\frac{{\partial }^2}{\partial {\eta }^2}\right)=\frac{1}{c^2({cosh}^2\xi -{cos}^2\eta )}\left(\frac{{\partial }^2}{\partial {\xi }^2}+\frac{{\partial }^2}{\partial {\eta }^2}\right).$$

(6)

2.2. Mathieu and Modified Mathieu Functions

The solution of the 2-dimensional wave equation in the elliptic coordinates $(\xi ,\eta )$,

$$\frac{1}{c^2({cosh}^2\xi -{cos}^2\eta )}\left(\frac{{\partial }^{2}w}{\partial {\xi }^2}+\frac{{\partial }^{2}w}{\partial {\eta }^2}\right)+k^{2}w=0,$$

(7)

can be obtained by employing the separation of variables method $w(\xi ,\eta )=X\left(\xi \right)Y\left(\eta \right)$. The separation yields the two ordinary differential equations of second order, namely the Mathieu and modified Mathieu equations that are expressed as

$$\begin{array}{c}\hfill \frac{d^{2}Y}{d{\eta }^2}+(a_q-2qcos2\eta )Y=0,\end{array}$$

(8)

$$\begin{array}{c}\hfill \frac{d^{2}X}{d{\xi }^2}-(a_q-2qcosh2\xi )X=0,\end{array}$$

(9)

respectively. The parameter $a_q$ denotes the separation constant which varies on the parameter $q=k^2{c}^2/4$. With regard to the elliptic channel, the appropriate solutions of Equations (8) and (9) have to satisfy the symmetric properties of the major and minor axes of the ellipse. Hence, the solution of the 2-dimensional wave Equation (7) in the elliptic channel is

$$w(\xi ,\eta )=\sum _{m=0}^{\infty }{D}_{2m}{Ce}_{2m}(\xi ;q){ce}_{2m}(\eta ;q),$$

(10)

where ${ce}_{2m}(\eta ;q)$ denotes the periodic Mathieu functions of integral order with period $\pi$ and ${Ce}_{2m}(\xi ;q)$ denotes the periodic modified Mathieu functions corresponding to ${ce}_{2m}(\eta ;q)$.

3. Mathematical Modeling

In this section, the governing equations for combined pressure-driven electroosmotic slip flow through an elliptic microchannel are constructed. Deriving from the Cauchy momentum equation and the mass conservation law, the motion of an incompressible fluid can be described by the Navier–Stokes equations coupled with the continuity equation:

$$\begin{array}{cc}\hfill \rho \left(\frac{\partial \overrightarrow{v}}{\partial t}+(\overrightarrow{v}·\nabla )\overrightarrow{v}\right)& =-\nabla p+\mu {\nabla }^2\overrightarrow{v}+{\overrightarrow{f}}_{\mathrm{EOF}},\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \nabla ·\overrightarrow{v}& =0,\hfill \end{array}$$

(12)

where $\overrightarrow{v}$ denotes the 3-dimensional vector of the fluid velocity, $\rho$ denotes the fluid density, p denotes the pressure, $\mu$ denotes the fluid viscosity, and ${\overrightarrow{f}}_{\mathrm{EOF}}$ denotes the electroosmotic body force. In the electrostatic theory, the relationship between the ionic charge density of the fluid ${\rho }_e$, and the EDL potential inside the channel $\psi$ is defined in the form of the Poisson equation as ${\nabla }^2\psi =-{\rho }_e/\epsilon$. Thus, the electroosmotic body force can be described by

$${\overrightarrow{f}}_{\mathrm{EOF}}={\rho }_e\overrightarrow{E}=-\epsilon \left({\nabla }^2\psi \right)\overrightarrow{E},$$

(13)

where $\overrightarrow{E}$ denotes the vector of external electric field and $\epsilon$ denotes the fluid permittivity.

Considering that the elliptic channel is placed along the z-axis direction and the fluid flow through the channel is fully developed, we then have the fluid velocity, the external electric field, and the pressure gradient not depending on time and existing only in the z-direction; that is, $\overrightarrow{v}=(0,0,u)$, $\overrightarrow{E}=(0,0,E)$, and $\nabla p=(0,0,p_z)$. Therefore, combining Equations (11)–(13) yields the governing equation as follows:

$$\frac{1}{c^2({cosh}^2\xi -{cos}^2\eta )}\left(\frac{{\partial }^{2}u}{\partial {\xi }^2}+\frac{{\partial }^{2}u}{\partial {\eta }^2}\right)=\frac{p_z}{\mu }+\frac{\epsilon E}{\mu c^2({cosh}^2\xi -{cos}^2\eta )}\left(\frac{{\partial }^2\psi }{\partial {\xi }^2}+\frac{{\partial }^2\psi }{\partial {\eta }^2}\right),$$

(14)

where E and $p_z$ are constant and $(\xi ,\eta )$ denotes the elliptic coordinate system. Additionally, the continuity Equation (12) is now reduced to $\partial u/\partial z=0$, and hence the z-axis velocity u depends only on $\xi$ and $\eta$. Moreover, suppose that the gradient of the EDL potential exists only in the direction perpendicular to the channel surface; then $\partial \psi /\partial z=0$, which gives rise to $\psi$, is also the function only of $\xi$ and $\eta$.

4. Boundary Conditions

Since the governing equation for combined pressure-driven electroosmotic slip flow through an elliptic microchannel was derived, in this section, we state the boundary conditions used in the study.

In regard to the symmetry of an elliptic channel, the fluid flow and the EDL potential distribution are considered to be symmetric about the x- and y-axes, i.e.,

$$\begin{array}{cccc}\hfill \frac{\partial u}{\partial \xi }(0,\eta )& =\frac{\partial \psi }{\partial \xi }(0,\eta )=0,\hfill & \hfill & 0\le \eta \le 2\pi ,\hfill \end{array}$$

(15a)

$$\begin{array}{cccc}\hfill \frac{\partial u}{\partial \eta }(\xi ,0)& =\frac{\partial \psi }{\partial \eta }(\xi ,0)=0,\hfill & \hfill & 0\le \xi \le {\xi }_0,\hfill \end{array}$$

(15b)

$$\begin{array}{cccc}\hfill \frac{\partial u}{\partial \eta }(\xi ,\frac{\pi }{2})& =\frac{\partial \psi }{\partial \eta }(\xi ,\frac{\pi }{2})=0,\hfill & \hfill & 0\le \xi \le {\xi }_0.\hfill \end{array}$$

(15c)

According to the appearance of slip flow in the microscale system, the Navier slip condition is employed to the fluid on the channel boundary as

$$u({\xi }_0,\eta )+\frac{l}{c\sqrt{{cosh}^2{\xi }_0-{cos}^2\eta }}\frac{\partial u}{\partial \xi }({\xi }_0,\eta )=0,$$

(15d)

where l denotes the Navier slip length (see Figure 3), and ${\xi }_0=ln\left((1+\sqrt{1-e^2})e^{-1}\right)$ denotes the boundary interface of the channel. For the EDL potential on the surface of the channel, we suppose that the charge on the surface of the channel is uniform with a constant zeta potential $\zeta$. Hence,

$$\psi ({\xi }_0,\eta )=\zeta .$$

(15e)

All equations relating to the boundary conditions and the governing Equation (14) form the boundary value problem that will be used to find the solution for fluid velocity.

5. Solutions of the Boundary Value Problem

To find fluid velocity in an elliptic microchannel, we solve the boundary value problem consisting of the partial differential Equation (14) and the boundary conditions (15). By using the identities ${cosh}^2\xi =(1+cosh2\xi )/2$ and ${cos}^2\eta =(1+cos2\eta )/2$, the governing Equation (14) can be rewritten as

$$\frac{2}{c^2(cosh2\xi -cos2\eta )}\left(\frac{{\partial }^{2}u}{\partial {\xi }^2}+\frac{{\partial }^{2}u}{\partial {\eta }^2}\right)=\frac{p_z}{\mu }+\frac{2\epsilon E}{\mu c^2(cosh2\xi -cos2\eta )}\left(\frac{{\partial }^2\psi }{\partial {\xi }^2}+\frac{{\partial }^2\psi }{\partial {\eta }^2}\right).$$

(16)

The non-homogeneous Equation (16) can be transformed into a homogeneous equation by employing the change of variable technique. Let

$$u(\xi ,\eta )=w(\xi ,\eta )+\frac{\epsilon E}{\mu }\psi (\xi ,\eta )+\frac{c^2{p}_z}{8\mu }(cosh2\xi +cos2\eta ).$$

(17)

Then Equation (14) is homogenized to

$$\frac{2}{c^2(cosh2\xi -cos2\eta )}\left(\frac{{\partial }^{2}w}{\partial {\xi }^2}+\frac{{\partial }^{2}w}{\partial {\eta }^2}\right)=0.$$

(18)

The boundary conditions for the function $w(\xi ,\eta )$ can be derived from Equations (15a)–(15c) as

$$\frac{\partial w}{\partial \xi }(0,\eta )=0,\frac{\partial w}{\partial \eta }(\xi ,0)=0,\mathrm{and}\frac{\partial w}{\partial \eta }(\xi ,\frac{\pi }{2})=0.$$

(19)

Since Equation (18) is homogeneous, we can apply the method of separation of variables to find the solution by letting $w(\xi ,\eta )=X\left(\xi \right)Y\left(\eta \right)$, where $X\left(\xi \right)$ denotes a function depending solely on $\xi$ and $Y\left(\eta \right)$ denotes a function depending solely on $\eta$. As a result, the solution that satisfies Equation (18) and the boundary conditions (19) is

$$w(\xi ,\eta )=\sum _{m=0}^{\infty }{C}_{2m}cosh2m\xi cos2m\eta ,$$

(20)

where $C_{2m}$ are constant. Thus, we obtain the solution to the boundary value problem as follows:

$$u(\xi ,\eta )=\sum _{m=0}^{\infty }{C}_{2m}cosh2m\xi cos2m\eta +\frac{\epsilon E}{\mu }\psi (\xi ,\eta )+\frac{c^2{p}_z}{8\mu }(cosh2\xi +cos2\eta ).$$

(21)

The coefficients $C_{2m}$ can be calculated by substituting the velocity function u expressed in Equation (21) into the Navier slip condition (15d) and using the zeta potential condition (15e) to arrive at

$$\begin{array}{cc}\hfill 0=& \sum _{m=0}^{\infty }{C}_{2m}\left(cosh2m{\xi }_0+2mg(\eta ,l)sinh2m{\xi }_0\right)cos2m\eta \hfill \\ \hfill & +h(\eta ,l)+\frac{c^2{p}_z}{8\mu }(cosh2{\xi }_0+cos2\eta +2g(\eta ,l)sinh2{\xi }_0),\hfill \end{array}$$

(22)

where

$$\begin{array}{cc}\hfill h(\eta ,l)& =\frac{\epsilon E}{\mu }\left(\zeta +g(\eta ,l)\frac{\partial \psi }{\partial \xi }({\xi }_0,\eta )\right),\hfill \end{array}$$

(23)

and

$$\begin{array}{cc}\hfill g(\eta ,l)& =\frac{l}{c\sqrt{{cosh}^2{\xi }_0-{cos}^2\eta }}.\hfill \end{array}$$

(24)

For $l=0$, we have $g(\eta ,l)=0$, and the Navier slip condition is reduced to the no-slip condition, i.e., $u(0,\eta )=0$. Then the coefficients $C_{2m}$ can be derived by using the orthogonality of the function $cos2m\eta$ on an interval $[0,2\pi ]$ as

$$C_{2m}=-\frac{{\int }_0^{2\pi }\left(\epsilon E\zeta {\mu }^{-1}+c^2{p}_z{\left(8\mu \right)}^{-1}(cosh2{\xi }_0+cos2\eta )\right)cos2m\eta d\eta }{{\int }_0^{2\pi }cosh2m{\xi }_0{cos}^{2}2m\eta d\eta }.$$

(25)

For $l>0$, let k be a non-negative integer. Multiplying $cos2k\eta$ into Equation (22) and integrating from 0 to $2\pi$ with respect to $\eta$ yields

$$\sum _{m=0}^{\infty }{A}_{k,m}{C}_{2m}=B_k,k=0,1,2,\dots ,$$

(26)

where

$$\begin{array}{c}{A}_{k,m}={\int }_0^{2\pi }(cosh2m{\xi }_0+2mg(\eta ,l)sinh2m{\xi }_0)cos2m\eta cos2k\eta d\eta ,\hfill \end{array}$$

(27)

$$\begin{array}{c}{B}_k=-{\int }_0^{2\pi }\left(h(\eta ,l)+c^2{p}_z{\left(8\mu \right)}^{-1}(cosh2{\xi }_0+cos2\eta +2g(\eta ,l)sinh2{\xi }_0)\right)cos2k\eta d\eta .\hfill \end{array}$$

(28)

To find the coefficients $C_{2m}$, the fluid velocity u is approximated by

$$u(\xi ,\eta )=\sum _{m=0}^M{C}_{2m}cosh2m\xi cos2m\eta +\frac{\epsilon E}{\mu }\psi (\xi ,\eta )+\frac{c^2{p}_z}{8\mu }(cosh2\xi +cos2\eta ),$$

(29)

for a suitable fixed positive integer M. Hence, the infinite system (26) is reduced to the system with finite equations as

$$\sum _{m=0}^M{A}_{k,m}{C}_{2m}=B_k,k=0,1,2,\dots ,M,$$

(30)

and the coefficients $C_{2m}$ are given by solving the reduced finite system (30).

The solution $u(\xi ,\eta )$ in Equation (21) still contains the unknown EDL potential function $\psi$. To find $\psi$, we assume the fluid is the symmetric electrolyte. By omitting the convection from the fluid, the distribution of the ionic charge density ${\rho }_e$ inside the channel can be described by the Boltzmann equation:

$${\rho }_e=-2n_0{p}^{+}{z}_{\nu }sinh\frac{p^{+}{z}_{\nu }}{k_{B}T}\psi ,$$

(31)

where $n_0$ denotes the concentration of ions at the bulk, $p^{+}$ denotes the elementary proton charge, $z_{\nu }$ denotes the valence of the ions, $k_B$ denotes the Boltzmann constant, and T denotes the fluid absolute temperature. Substituting ${\rho }_e$ from the Boltzmann Equation (31) into the Poisson equation ${\nabla }^2\psi =-{\rho }_e/\epsilon$ yields the Poisson–Boltzmann equation:

$${\nabla }^2\psi (\xi ,\eta )=-\frac{{\rho }_e}{\epsilon }=\frac{2n_0{p}^{+}{z}_{\nu }}{\epsilon }sinh\frac{p^{+}{z}_{\nu }}{k_{B}T}\psi .$$

(32)

Suppose that zeta potential is low ($|p^{+}{z}_{\nu }\psi |<k_{B}T$) [38]; we can approximate

$$sinh\frac{p^{+}{z}_{\nu }}{k_{B}T}\psi \approx \frac{p^{+}{z}_{\nu }}{k_{B}T}\psi .$$

(33)

Therefore, the Poisson–Boltzmann Equation (32) is reduced to the Debye–Hückel approximation in elliptic coordinates $(\xi ,\eta )$:

$$\frac{1}{c^2({cosh}^2\xi -{cos}^2\eta )}\left(\frac{{\partial }^2\psi }{\partial {\xi }^2}+\frac{{\partial }^2\psi }{\partial {\eta }^2}\right)-{\kappa }^2\psi =0,$$

(34)

where $\kappa ={(2n_0{p}^{{+}^2}{z}_{\nu }^2/\epsilon k_{B}T)}^{1/2}$, which is known as the reciprocal of the EDL thickness.

Since Equation (34) is the 2-dimensional wave equation in the elliptic coordinate system, the solution that satisfies the boundary conditions (15a)–(15c) can be derived in the form of the Mathieu functions [30] as

$$\psi (\xi ,\eta )=\sum _{m=0}^{\infty }{D}_{2m}{Ce}_{2m}(\xi ;-q){ce}_{2m}(\eta ;-q),$$

(35)

where $D_{2m}$ is constant, ${Ce}_{2m}(\xi ;-q)$ denotes the periodic Mathieu functions with $q={\kappa }^2{c}^2/4$, and ${ce}_{2m}(\eta ;-q)$ denotes the modified Mathieu functions corresponding to ${Ce}_{2m}(\xi ;-q)$. Substituting the function $\psi$ into the boundary condition (15e), then we can use the orthogonality of the Mathieu functions ${ce}_{2m}(\eta ;-q)$ on the interval $[0,2\pi ]$ to determine the coefficients $D_{2m}$ as follows:

$$D_{2m}=\frac{{\int }_0^{2\pi }\zeta ·{ce}_{2m}(\eta ;-q)d\eta }{{\int }_0^{2\pi }{Ce}_{2m}({\xi }_0;-q){ce}_{2m}^2(\eta ;-q)d\eta }.$$

(36)

To non-dimensionalize Equations (29) and (35), we define the non-dimensional variables

$$u^{\ast }=\frac{u}{-\epsilon \zeta E{\mu }^{-1}},{\psi }^{\ast }=\frac{\psi }{\zeta },p_z^{\ast }=\frac{p_z}{8\epsilon {\zeta }^2{c}^{-3}},E^{\ast }=\frac{E}{\zeta c^{-1}}.$$

(37)

Thus, the dimensionless velocity profile $u^{\ast }$ is expressed in the form of

$$u^{\ast }(\xi ,\eta )=\sum _{m=0}^M{C}_{2m}^{\ast }cosh2m\xi cos2m\eta -{\psi }^{\ast }(\xi ,\eta )-\frac{p_z^{\ast }}{E^{\ast }}(cosh2\xi +cos2\eta ),$$

(38)

where $C_{2m}^{\ast }=C_{2m}/(-\epsilon \zeta E{\mu }^{-1})$ are given for $l=0$ by

$$C_{2m}^{\ast }=\frac{{\int }_0^{2\pi }\left(1+p_z^{\ast }{\left(E^{\ast }\right)}^{-1}(cosh2{\xi }_0+cos2\eta )\right)cos2m\eta d\eta }{{\int }_0^{2\pi }cosh2m{\xi }_0{cos}^{2}2m\eta d\eta },$$

(39)

and are obtained for $l>0$ by solving the system

$$\sum _{m=0}^M{A}_{k,m}{C}_{2m}^{\ast }=B_k^{\ast },k=0,1,2,\cdots ,M,$$

(40)

where

$$\begin{array}{c}{B}_k^{\ast }={\int }_0^{2\pi }\left(p_z^{\ast }{\left(E^{\ast }\right)}^{-1}(cosh2{\xi }_0+cos2\eta +2g(\eta ,l)sinh2{\eta }_0)+h^{\ast }(\eta ,l)\right)cos2k\eta d\eta ,\end{array}$$

(41)

$$\begin{array}{c}{h}^{\ast }(\eta ,l)=1+g(\eta ,l)\frac{\partial {\psi }^{\ast }}{\partial \xi }({\xi }_0,\eta ).\end{array}$$

(42)

The non-dimensional EDL potential ${\psi }^{\ast }$ can be derived as

$${\psi }^{\ast }(\xi ,\eta )=\sum _{m=0}^{\infty }{D}_{2m}^{\ast }{Ce}_{2m}(\xi ;-q){ce}_{2m}(\eta ;-q),$$

(43)

where

$$D_{2m}^{\ast }=\frac{{\int }_0^{2\pi }{ce}_{2m}(\eta ;-q)d\eta }{{\int }_0^{2\pi }{Ce}_{2m}({\xi }_0;-q){ce}_{2m}^2(\eta ;-q)d\eta }.$$

(44)

6. Effect of the Eccentricity on Slip Flow

From a fluidic aspect, a volumetric flow rate instead of velocity plays a more crucial role in flow manipulation [37]. As a consequence, to investigate the effect of the eccentricity, we consider the volumetric flow rate Q per unit area of the channel cross-section defined by

$$Q=\frac{1}{A}\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{{\xi }_0}{\int }}{c}^2({cosh}^2\xi -{cos}^2\eta )u(\xi ,\eta )d\xi d\eta ,$$

(45)

where A denotes the area of the cross-section and u is calculated via Equation (29) for the appropriate number $M=6$. Using the non-dimensional velocity $u^{\ast }$ defined in Equation (37), the dimensionless flow rate $Q^{\ast }$ can be calculated by the following equation:

$$Q^{\ast }=\frac{Q}{\epsilon \zeta E{\mu }^{-1}}=\frac{c^2}{A}\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{{\xi }_0}{\int }}({cosh}^2\xi -{cos}^2\eta )u^{\ast }(\xi ,\eta )d\xi d\eta .$$

(46)

The effect of the eccentricity is investigated in two aspects of the channel geometry: (i) the area of the elliptic cross-section of the channel is fixed and (ii) the perimeter of the elliptic cross-section of the channel is fixed. When the area of the cross-section is fixed, the focal length of the ellipse can be derived relating to the eccentricity e by

$$c=e\sqrt{\frac{A}{\pi \sqrt{1-e^2}}}.$$

(47)

On the other hand, the focal length can be calculated for each eccentricity by fixing the channel perimeter as follows:

$$c=\frac{eP}{4{\int }_0^{\pi /2}\sqrt{1-e^{2}sin\eta }d\eta },$$

(48)

where P denotes the perimeter of the elliptic cross-section. The properties of the fluid [38] and the other parameters used in the investigation are presented in

Table 1. Table of parameters used in the investigation.

Name	Symbol	Value	SI Unit
Fluid viscosity	$\mu$	$\phantom{-}9.00\times {10}^{-4}$	Pa s
Fluid permittivity	$\epsilon$	$\phantom{-}6.95\times {10}^{-10}$	F m−1
Pressure gradient in z-axis	$p_z$	$-2.00$	Pa m−1
Reciprocal of EDL thickness	$\kappa$	$\phantom{-}8.00\times {10}^4$	m−1
Zeta potential	$\zeta$	$-2.49\times {10}^{-4}$	$\mathrm{V}$
External electric field	E	$\phantom{-}5.00\times {10}^2$	V m−1

unless stated specifically otherwise.

Figure 4a presents the comparison of the dimensionless volumetric flow rates $Q^{\ast }$ through the elliptic channel with various eccentricities when the area of the cross-section is fixed to $4.5\pi \times {10}^3$ $\mathsf{\mu }{\mathrm{m}}^2$. The slip flow at the channel wall is considered with the slip length $l=0,10,100$, and 1000 $\mathsf{\mu }\mathrm{m}$. The result shows that when $l=0$, which is the no-slip wall, no significant difference in $Q^{\ast }$ exists in the channel with the eccentricity less than 0.5. However, when the eccentricity is greater than 0.5, we can see a considerable decrease in $Q^{\ast }$ as a rise in the eccentricity. When e increases, the shortest distance between the channel center and the perimeter (known as the semi-minor axis) of an ellipse decreases (see Figure 5a). The reduction of this distance results in the velocity profile of fluid. In particular, the less the distance, the lower the maximum velocity through the channel (see Figure 6). Hence, an increase in the eccentricity yields a lower volumetric flow rate. This result is well consistent with Numpanviwat’s work [32]. When the slip flow phenomenon is applied, i.e., $l\ne 0$, $Q^{\ast }$ increases as an increase in the slip length. In spite of that, for each slip length, there is no crucial difference in $Q^{\ast }$ among the channel with $e<0.5$, but $Q^{\ast }$ dramatically drops as an increase in the eccentricity when $e>0.5$. This result resembles that of the no-slip flow.

To investigate the effect of the eccentricity on slip flow when the perimeter is fixed, the variation of the non-dimensional flow rates $Q^{\ast }$ in the elliptic cross-sectional channel where the cross-section perimeter is fixed to $4.127\times {10}^2$ $\mathsf{\mu }\mathrm{m}$ is demonstrated in Figure 4b. The eccentricities and the slip lengths used in this investigation are the same as in the former comparison when the area of the cross-sections is fixed. The result shows that if the slip length increases, then so does $Q^{\ast }$. However, $Q^{\ast }$ fluctuates significantly when the eccentricity increases. Additionally, the flow rate increases gradually when $0<e<0.6$, and the flow rate decreases sharply when $e>0.6$. When the perimeter of an ellipse is fixed, the variation of the dimensionless flow rate $Q^{\ast }$ on the eccentricity e should be similar to when the area of the elliptic cross-section is fixed because the semi-minor axis also decreases as an increase in e (see Figure 5b). However, with a fixed perimeter, the area of an ellipse rises when $0<e<0.6$, as shown in Figure 7. Therefore, it yields a higher flow rate. On the other hand, when $e>0.6$, the area of the ellipse sharply drops as an increase in e. As a consequence, the flow rate decreases significantly.

In order to overcome the complexity of calculating the numerical result, one may treat an elliptic cross-sectional channel as a cylindrical channel. Consequently, it is worth examining here how the eccentricity affects the flow rate in the relative error aspect compared with the flow rate in the equivalent cylindrical channel. Figure 8 illustrates the comparison of the relative error of the non-dimensional flow rate in the elliptic channel with various eccentricities compared with the non-dimensional flow rate in the circular channel when either the area or the perimeter of the cross-sections is fixed. The relative error $\delta$ of the flow rate is defined by

$${\delta }_i=\left|\frac{Q_0^{\ast }-Q_i^{\ast }}{Q_i^{\ast }}\right|\times 100,$$

(49)

where $Q_0^{\ast }$ denotes the dimensionless flow rate in the circular channel and $Q_i^{\ast }$ denotes the dimensionless flow rate in the elliptic channel with eccentricity equal to i. The result in Figure 8a shows that, when the area is fixed, the relative errors of the no-slip flow in the channel with $e<0.5$ are less than $1\%$. When the slip length increases, the fluid velocity increases, but the velocity profile is still in the same pattern, as shown in Figure 6. Hence, the change in the volumetric flow rate on an increase in the eccentricity is practically identical among the channel with different slip lengths. This phenomenon can be observed in Figure 4a as the pattern of the reductions in $Q^{\ast }$ of the channel with four different slip lengths as e increases are strongly resembled. As a consequence, there is no crucial difference in the absolute error of the flow rate defined by $|Q_0^{\ast }-Q_i^{\ast }|$ between the channels having any two different slip lengths. Since the relative error is defined by

$${\delta }_i=\left|\frac{Q_0^{\ast }-Q_i^{\ast }}{Q_i^{\ast }}\right|\times 100=\frac{\mathrm{absolute}\mathrm{error}}{|Q_i^{\ast }|}\times 100,$$

(50)

a rise in slip length leads to a higher flow rate. However, there is no significant difference that occurs in the absolute error as the aforementioned argument. Thus, the relative error reduces when the slip length increases.

On the other hand, when the perimeter is fixed, Figure 8b demonstrates that the relative errors of the no-slip flow are notably high in the elliptic channel having $e>0.2$. The relative error is well consistent with the fluctuation of the volumetric flow rate due to the change in the cross-sectional area.

Therefore, when the slip flow is considered, the relative error decreases with an increase in the slip length. As a result, the simplification by using the equivalent circular channel for the elliptic channel can be done under the condition that the area of the cross-sections is fixed and the eccentricity is less than 0.5.

Furthermore, we investigate the relative error of the flow rate with various ratios of the non-dimensional pressure gradient $p_z^{\ast }$ and the non-dimensional external electric field $E^{\ast }$. According to the aforementioned result, we consider only the relative error when the area of the cross-sections is fixed. Figure 9 demonstrates the comparison of the relative errors of the flow rate with the ratio $p_z^{\ast }/E^{\ast }$ equal to 1/5, 1/2.5, 2.5, and 5 times that ratio used in the investigation in Figure 8. The results show that the relative errors enlarge when $p_z^{\ast }/E^{\ast }$ increases in all the channels with different values of the eccentricity. This phenomenon can be described as relating to the slip flow. When the ratio $p_z^{\ast }/E^{\ast }$ increases, the flow is more dominated by the pressure-gradient force compared to the EOF. It means, on the other hand, the EOF is less dominant in the flow. Since the EOF is acting on the EDL, which is the area adjacent to the wall, the weaker EOF will result similar to the smaller slip length. As a consequence, the relative error of the flow rate increases with an increase in the ratio $p_z^{\ast }/E^{\ast }$ in the same manner as a decrease in the slip length. Moreover, the effect of the ratio $p_z^{\ast }/E^{\ast }$ vanishes when the slip length is very large, i.e., $l>100$ $\mathsf{\mu }\mathrm{m}$.

7. Conclusions

We present the solution for the fluid velocity of slip flow driven by both pressure-gradient force and EOF through an elliptic cylindrical microchannel. The EOF is obtained by solving the Debye–Hückel approximation, and the fluid velocity is derived based on the Navier–Stokes and the continuity equations under the Navier slip condition. The volumetric flow rate is calculated numerically to investigate validation of the consideration of the elliptic cylindrical microchannel to be the equivalent circular cylindrical microchannel. The results show that the consideration of the equivalent circular channel provides a valid result of the flow rate under the condition that the area of the channel cross-section is fixed and the eccentricity is less than 0.5. Additionally, the consideration yields a more accurate flow rate in slip electroosmotic flow driven by an external electric field.