Transient Pressure-Driven Electroosmotic Flow through Elliptic Cross-Sectional Microchannels with Various Eccentricities

: The semi-analytical solution for transient electroosmotic ﬂow through elliptic cylindrical microchannels is derived from the Navier-Stokes equations using the Laplace transform. The electroosmotic force expressed by the linearized Poisson-Boltzmann equation is considered the external force in the Navier-Stokes equations. The velocity ﬁeld solution is obtained in the form of the Mathieu and modiﬁed Mathieu functions and it is capable of describing the ﬂow behavior in the system when the boundary condition is either constant or varied. The ﬂuid velocity is calculated numerically using the inverse Laplace transform in order to describe the transient behavior. Moreover, the ﬂow rates and the relative errors on the ﬂow rates are presented to investigate the effect of eccentricity of the elliptic cross-section. The investigation shows that, when the area of the channel cross-sections is ﬁxed, the relative errors are less than 1% if the eccentricity is not greater than 0.5. As a result, an elliptic channel with the eccentricity not greater than 0.5 can be assumed to be circular when the solution is written in the form of trigonometric functions in order to avoid the difﬁculty in computing the Mathieu and modiﬁed Mathieu functions. parameters from the ﬁrst and the second parts of the investigation on the ﬂuid velocity, respectively. The results show that, when the area is ﬁxed, there is no signiﬁcant difference for the ﬂow rate in the channels with the eccentricity between 0 and 0.5. When the eccentricity is greater than 0.5; however, the ﬂow rate decreases rapidly as the eccentricity increases. On the other hand, when the circumference is ﬁxed, the ﬂow rate increases signiﬁcantly when the eccentricity is between 0.2 and 0.6 but the ﬂow rate decreases dramatically when the eccentricity is greater than 0.6.


Introduction
A manipulation on fluid flow under a microscale system has been studied intensively worldwide over the past decades [1][2][3][4][5][6] due to the utility of its applications in various fields; for example, lab-on-the-chips, portable medical devices, drug delivery devices and computer chips. One key to success in controlling fluid behavior through a very small scale channel is to employ an application of the well-known electrokinetic phenomenon, namely electroosmosis. Electroosmosis is a phenomenon of the fluid movement through a channel by the electrokinetic force, which relies on two layers of ions [7]. The first layer is the ions on the inner wall surface of the channel due to the chemical reaction between the channel surface and the fluid. The second layer is the counterions in the fluid adjacent to the channel wall as a result of the electrokinetic force from the first layer. The formation of these two electrical layers is called the electrical double layer (EDL) shown by red highlighted area in Figure 1. An application of an external electric field to the channel leads to a movement of the counterions in the fluid by electrokinetic force called as electroosmotic force. Therefore, the fluid flow in consequence of the drag force from the moving counterions induced by electroosmotic force is named as electroosmotic flow (EOF).
Since EDL locates on the interface of the channel and the fluid, the cross-sectional shape of cylindrical channel plays an important role on the fluid behavior [8]. Recently, understanding of fluid behavior through microchannels has been achieved using both experimental and theoretical studies [9][10][11][12][13][14]. In literature, the channels are considered as cylindrical microtubes with various cross-sections including rectangles [15,16], circles [17,18], annuli [19,20], ellipses [21,22] and other complex shapes [21,23], according to the occurrences of microfluidic applications.  With regard to the shape of many microfluidic devices and natural systems such as the vascular system, elliptic cylindrical geometry is considered to be the channel in many researches [21,22,[24][25][26] in order to study the behavior of electroosmotic flow. However, the solution for fluid velocity through an elliptic cylindrical channel using the mathematical methodology is obtained in the form of Mathieu functions [22] which have difficulty in numerically computing its value due to the convergence of the function [27]. In 2015, Liu et al. [28] studied the effect of eccentricity on pure electroosmotic flow through an elliptic channel. Pressure-driven force was neglected and the flow was considered to be fully developed which means the transient behavior was omitted. Furthermore, the effect of eccentricity was investigated only when the circumference of the channel cross-section was fixed.
Therefore, in this study, we derive the solution for transient combined electroosmotic and pressure-driven flow through an elliptic cylindrical microchannel. The governing equation for fluid velocity is written in the homogeneous linear partial differential equation depending on both ξ and η in the (ξ, η) elliptic coordinates. By employing the Laplace transform to eliminate the time-derivative term in the governing equation, the solution is obtained as the function of ξ and η. As a consequence, the solution is capable of describing the electroosmotic flow in the system when the boundary condition is either constant or varied on η, which normally occurs in an elliptic cross-sectional microchannel [4,21,22,29,30]. However, the solution of the Laplace transform which is in the form of the special functions namely the Mathieu and modified Mathieu functions has complexity to find the inverse Laplace transform analytically. Thus, the numerical inverse Laplace transform is employed to yield the semi-analytical solution. This semi-analytical solution for transient electroosmotic flow allows us to investigate behavior of electroosmotic flow directly via the expression of the solution and takes an advantage of making the problem tractable for variations of domain and boundary condition which provide more efficiency (speed and accuracy) in computing the numerical results compared to the traditional numerical method. Moreover, the velocity profiles of the fluid flow are calculated numerically from the semi-analytical solution in order to demonstrate the transient behavior of electroosmotic flow. The effect of eccentricity of the elliptic channel on flow characteristics including the velocity profile, the volumetric flow rate and the relative error on the volumetric flow rate is investigated with two aspects: the circumference of the channel cross-section is fixed and the area of the channel cross-section is fixed.

Elliptic Coordinate System
The elliptic coordinate system (ξ, η) for an elliptic geometry with two foci at (c, 0) and (−c, 0) of the Cartesian coordinate system (x, y) is an orthogonal coordinate system in which ξ represents the confocal ellipses and η represents the confocal hyperbolas. The elliptic coordinates can be illustrated as shown in Figure 2. The elliptic coordinates (ξ, η) are defined by where 0 ≤ ξ < ∞ and 0 ≤ η < 2π. The identities of the elliptic coordinates related to the Cartesian coordinates are x 2 c 2 cos 2 η − y 2 c 2 sin 2 η = 1, (4) and the Laplacian in the elliptic coordinates is

Mathieu Functions
Consider the 2-dimensional wave equation in the elliptic coordinates 1 c 2 (cosh 2 ξ − cos 2 η) By substituting W(ξ, η) = F(ξ)G(η) and the identity we can split Equation (6) into two second order ordinary differential equations as follows: where a is the separation constant and q = 0.25k 2 c 2 . The solutions of the ordinary differential Equations (8) and (9) are written in the form of two special functions introduced by French mathematician, Émile Léonard Mathieu in 1868 [31]. The special functions are called as the Mathieu functions and modified Mathieu functions. Combining the solutions of Equations (8) and (9), we obtain the solution of Equation (6) as where ce m (η; q), se m (η; q) are the periodic Mathieu functions; Ce m (η; q), Se m (η; q) are the periodic modified Mathieu functions corresponding to ce m (η; q), se m (η; q), respectively; Fe m (η; q), Ge m (η; q) are the non-periodic modified Mathieu functions corresponding to ce m (η; q), se m (η; q), respectively.
The coefficients C 1 m (q), C 2 m (q), S 1 m (q) and S 2 m (q) are functions depending solely on q, and the integer m denotes the order of the Mathieu and modified Mathieu functions. The formula of the above special functions can be found in [27].

Electroosmotic Force
The electroosmotic force is given by where E denotes the applied external electric field and ρ e denotes the ionic charge density which can be expressed by the Poisson-Boltzmann equation as follows: where ψ denotes the potential inside the channel, n 0 denotes the ionic concentration at the bulk, p + denotes the charge of a proton, z v denotes the valence of the ions in the fluid, ε denotes the permittivity of the fluid, k B is the Boltzmann constant, and T denotes the absolute temperature.
If the potential ψ has low value, we can simplify Equation (15) by using sinh x ≈ x to be the well-known Debye-Hückel approximation in the elliptic coordinates as where κ = (2n 0 p + 2 z 2 ν ε −1 k −1 B T −1 ) 1/2 represents the reciprocal of the EDL thickness. Combining Equations (14)-(16), we can write the electroosmotic force related to the potential ψ as Since Equation (16) is in the form of Equation (6), the general solution is obtained in the similar form as follows: According to the symmetric property of an elliptic geometry, the potential ψ is supposed to be symmetrical about both xand y-axes. Hence, we can express the boundary conditions of ψ in the elliptic coordinate system by the following equations: The general solution shown in Equation (18) can be reduced, according to the symmetric property of the potential ψ, as follows: Since the partial derivative of the modified Mathieu function Fe 2m (ξ; −q) with respect to ξ does not vanish when ξ = 0 [27], the constant A 2 2m (q) has to be zero in order to satisfy the boundary condition (21). Therefore, the general solution (22) is reduced to where A 2m (q) = A 1 2m (q). Moreover, if the potential at the wall surface is determined by the function f (η), we can write the boundary of the potential ψ at the boundary of the elliptic channel ξ = ξ 0 as Now, we apply the boundary condition (24) to the general solution (23). Thus, the orthogonality of the Mathieu functions ce 2m (η; −q) yields the constant A 2m as

Fluid Velocity
Suppose that the fluid flow does not swirl, we can describe the fluid behavior of electroosmotic flow through an elliptic cylindrical microchannel along z-axis by the Navier-Stokes equations with the electroosmotic force in the elliptic cylindrical coordinate system (ξ, η, z) as shown: where u(ξ, η, t) denotes the fluid velocity in z-direction, ρ denotes the fluid density, µ denotes the fluid viscosity, and p z denotes the constant pressure gradient in z-direction. Substituting the term of electroosmotic force in Equation (26) with the relation in Equation (17) yields To eliminate the time-derivative, we employ the Laplace transform L defined by and the initial condition u(ξ, η, 0) = 0. Therefore, Equation (27) becomes which has the general solution written in the form of where v c is the complementary solution (homogeneous solution), v p 1 is the particular solution corresponding to non-homogeneous term p z µ −1 s −1 , and v p 2 is the particular solution corresponding to non-homogeneous term εκ 2 Eψµ −1 s −1 .
For the complementary solution v c , we consider the homogeneous equation Since Equation (31) is in the form of Equation (6), the complementary solution v c is obtained in the similar form as follows: where B 1 m (q s ), B 2 m (q s ), B 3 m (q s ), and B 4 m (q s ) are constant with q s = 0.25ρc 2 sµ −1 . For the particular solution v p 1 , since the non-homogeneous term p z µ −1 s −1 does not depend on both ξ and η, we can assume that v p 1 is constant. By the method of undetermined coefficients, we obtain v p 1 as follows: For the particular solution v p 2 , since the non-homogeneous term εκ 2 Eψµ −1 s −1 is in the form of the Mathieu and modified Mathieu functions shown in Equation (23), we assume that the particular solution v p 2 is also in the same formula as Since v p 2 is the particular solution of the corresponding non-homogeneous equation 1 c 2 (cosh 2 ξ − cos 2 η) we have Moreover, substituting Equation (23) into Equation (16) and equating the term having the function A 2m (q) with the same integer m, then we obtain for all integer m ≥ 0. Observe that we can rewritten Equation (36) as By equating coefficients of the terms having the same order of the periodic Mathieu function ce 2m (η; −q), we can derive the constant C 2m for the particular solution v p 2 from for all integer m ≥ 0. Combining Equations (30), (32)-(34) yields the general solution of Equation (29) as follows: Similarly to the potential ψ, we also suppose that v is symmetrical about both xand y-axes. Therefore, the solution in Equation (41) In this study, we consider the no-slip condition of the fluid, that is, u(ξ 0 , η, t) = 0 which implies v(ξ 0 , η, s) = 0.
As a result, the constant B 1 2m can be determined using the orthogonality of the Mathieu function ce 2m (η; −q s ) by The fluid velocity u can be obtained by employing the inverse Laplace transform as follows: where γ is a number that Re(γ) is greater than the real part of all singularities of v(ξ, η, s) in the complex s-plane.
For steady flow, the time-derivative of the fluid velocity is vanished. Therefore, Equation (26) with the relation in Equations (14) and (15) can be reduced to 1 c 2 (cosh 2 ξ − cos 2 η) where u is now the function depending only on ξ and η. Let We can rewrite Equation (47) by using the identity in Equation (7) to be the Laplace equation of the function w(ξ, η) as 1 c 2 (cosh 2 ξ − cos 2 η) Using the method of separation of variables, we obtain where D 2m is constant. Substituting the general solution w(ξ, η) into Equation (48) yields where the constant D 2m can be derived from the no-slip condition u(ξ 0 , η, t) = 0 with the orthogonality of the trigonometric functions as

Numerical Results
In this section, we employ the numerical inverse Laplace transform on the general solution of fluid velocity as shown in Equation (42) in order to investigate the behavior of transient electroosmotic flow through an elliptic cylindrical microchannel. The calculation of ce 2m (z; q), Ce 2m (z; q), ce 2m (z; −q) and Ce 2m (z; −q) is presented in Appendix A. The dimensions of the elliptic cross-section, the fluid properties [8] and the other parameters used in the investigation are presented in Table 1.    After fully developing, the flow becomes steady. As a result, we use the velocity of steady flow shown in Equation (51) to investigate the effect of eccentricity on fluid flow through the elliptic channel. The eccentricity of ellipses varies from 0 (circle) to 0.9 by increment of 0.1. The variation of the elliptic cross-sections of the channel on various eccentricities is presented in Figure 4.  The investigation on the effect of eccentricity consists of two parts. For the first part, we fix the area A of all ellipses to be 4500π µm 2 (see Figure 4a). This value is obtained by calculating the area of the ellipse having the parameters as shown in Table 1. Then, for each eccentricity, we can calculate the focal length c and the boundary ξ 0 of the ellipse as follows: Figure 5 presents the comparison of the velocity profiles of electroosmotic flow through the elliptic channel with various eccentricities when the area is fixed. The velocity profiles along xand y-axes are shown in Figure 5a,b, respectively. The results show that there is no significant difference for the velocity profiles along both xand y-axes in the channels with the eccentricity between 0 and 0.6. However, when the eccentricity is greater than 0.6, the velocity drops along the entire xand y-axes as an increase of the eccentricity. For the second part, we fix the circumference C of all ellipses equal to about 412.7 µm (see Figure 4b). This value is also obtained by calculating the circumference of the ellipse having the parameters as shown Table 1. Then, for each eccentricity, we can calculate the boundary ξ 0 of the ellipse by Equation (54) and the focal length c as c = eC 2π 0 1 − e 2 sin η dη .
(55) Figure 6 presents the comparison of the velocity profiles of electroosmotic flow through the elliptic channel with various eccentricities when the circumference is fixed. The velocity profiles along xand y-axes are shown in Figure 6a,b, respectively. The results show that the pattern of the velocity profiles is quite similar compared to the one when the area is fixed. Indeed, the velocity profiles are not significantly different when the eccentricity is between 0 and 0.6. However, the velocity decreases along the entire xand y-axes as an increase of the eccentricity when the eccentricity is greater than 0.6. Since the efficiency of flow control in microfluidics mainly relies on the accuracy of flow rate, we also investigate the effect of eccentricity on the volumetric flow rate Q of the fluid through the elliptic cylindrical microchannel. The flow rates can be calculated by integrating the fluid velocity of steady flow, shown in Equation (51), over the cross-sectional area. Therefore, the formula of the flow rate is expressed as Figure 7 demonstrates the comparison of the flow rates calculated when the area (see Figure 7a) and the circumference (see Figure 7b) of ellipse are fixed by using the parameters from the first and the second parts of the investigation on the fluid velocity, respectively. The results show that, when the area is fixed, there is no significant difference for the flow rate in the channels with the eccentricity between 0 and 0.5. When the eccentricity is greater than 0.5; however, the flow rate decreases rapidly as the eccentricity increases. On the other hand, when the circumference is fixed, the flow rate increases significantly when the eccentricity is between 0.2 and 0.6 but the flow rate decreases dramatically when the eccentricity is greater than 0.6.  Furthermore, we define the relative error δ of the flow rate by where Q i is the flow rate with the eccentricity equal i. The comparison of the relative errors of the flow rates is presented in Figure 8 when the area and the circumference of the ellipse are fixed (Figure 8a,b, respectively). In addition, to investigate the effect of external electric field on the relative error, we consider the variation of the external electric field on the electroosmotic flow to be 0.1 and 10 times of the value shown in Table 1. Indeed, we investigate the relative error using the three different values of the external electric field E = 50, 500 and 5000 V m −1 , respectively.  The results show that, when the area is fixed, the relative errors of all three cases are less than 1% in the channel with the eccentricity not greater than 0.5. However, when the eccentricity is greater than 0.5, the relative errors increase sharply as an increase of the eccentricity. On the other hand, when the circumference is fixed, the relative errors are less than 1% if the eccentricity is not greater than 0.2. However, there is a fluctuation in the relative errors with the value greater than 1% when the eccentricity is greater than 0.2.
Moreover, the results show that, in all three cases, the relative errors increase as a decrease of the external electric field.
Since the relative error defined by Equation (57) can be used to describe the error of the flow rate in the elliptic channel when we consider it as the circular channel. Therefore, the results in Figure 8 also show that the consideration of an elliptic channel as the circular channel with the same circumference of the cross-section may cause high error in the flow rate which is inappropriate to use for any investigation. This influence of eccentricity agrees with that of Liu [28].
However, this study reveals a significant result that, when the area of cross-section is fixed, we can assume that the elliptic channel with the eccentricity not greater than 0.5 to be the circular channel. This assumption can be used to avoid the difficulty in computing the velocity solution for electroosmotic flow in the elliptic channel, which is in the form of the Mathieu and modified Mathieu functions. Moreover, it provides a simpler way to further the numerical investigation of electroosmotic flow in some elliptic systems.

Conclusions
In this study, we present the semi-analytical solution for transient pressure-driven electroosmotic flow through elliptic cylindrical microchannels based on the Mathieu and modified Mathieu functions. The velocity profiles of transient pressure-driven electroosmotic flow are calculated numerically in order to investigate the transient behavior. Moreover, the variations of the fluid velocities, the flow rates and the relative errors of the flow rates with various eccentricities of the channel cross-section are presented to investigate the effect of eccentricity of the elliptic cross-sectional microchannel. The results show that the elliptic channel with the eccentricity not greater than 0.5 can be considered as the circular channel with the same area of the cross-section to avoid the difficulty in numerically computing the solution for fluid velocity in the elliptic channel, which is based on the Mathieu and modified Mathieu functions. In addition, this assumption provides a simpler way of furthering numerical investigation of electroosmotic flow in some elliptic systems. Data Availability Statement: Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The procedure to calculate the Mathieu and modified Mathieu functions including ce 2m (z; q), Ce 2m (z; q), ce 2m (z; −q) and Ce 2m (z; −q), which are in the formula of the solution, is presented as follows [27]: For q > 0 and integer m ≥ 0, the functions ce 2m (z; q), Ce 2m (z; q), ce 2m (z; −q) and Ce 2m (z; −q) can be expressed by ce 2m (z; q) = To find the coefficients A (2m) 2r , substituting Equation (A1) into Equation (8) and equating the coefficients of the term cos 2rz for each r = 0, 1, 2, . . . yield the following recurrence relations: where a 2m is the separation constant a corresponding to ce 2m . Dividing Equations (A5) and (A6) by A