Next Article in Journal
Reduced Model for Properties of Multiscale Porous Media with Changing Geometry
Next Article in Special Issue
Origin of Irrational Numbers and Their Approximations
Previous Article in Journal
Modelling of Conditional Scalar Dissipation Rate in Turbulent Premixed Combustion
Previous Article in Special Issue
Improved Stability Criteria on Linear Systems with Distributed Interval Time-Varying Delays and Nonlinear Perturbations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Nattakarn Numpanviwat
1,† and
Pearanat Chuchard
2,*,†
1
Applied Mathematics Program, Faculty of Science and Technology, Valaya Alongkorn Rajabhat University under the Royal Patronage, Pathum Thani 13180, Thailand
2
Department of Mathematics, Faculty of Science, Ramkhamhaeng University, Bangkok 10240, Thailand
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Submission received: 31 January 2021 / Revised: 21 February 2021 / Accepted: 24 February 2021 / Published: 1 March 2021

Abstract

:
The semi-analytical solution for transient electroosmotic flow 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 field solution is obtained in the form of the Mathieu and modified Mathieu functions and it is capable of describing the flow behavior in the system when the boundary condition is either constant or varied. The fluid velocity is calculated numerically using the inverse Laplace transform in order to describe the transient behavior. Moreover, the flow rates and the relative errors on the flow 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 fixed, 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 difficulty in computing the Mathieu and modified Mathieu functions.

1. 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.

2. Preliminaries on the Elliptic Coordinate System

2.1. 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
x = c cosh ξ cos η ,
y = c sinh ξ sin η ,
where 0 ξ < and 0 η < 2 π . The identities of the elliptic coordinates related to the Cartesian coordinates are
x 2 c 2 cosh 2 ξ + y 2 c 2 sinh 2 ξ = 1 ,
x 2 c 2 cos 2 η y 2 c 2 sin 2 η = 1 ,
and the Laplacian in the elliptic coordinates is
2 = 1 c 2 ( cosh 2 ξ cos 2 η ) 2 ξ 2 + 2 η 2 .

2.2. Mathieu Functions

Consider the 2-dimensional wave equation in the elliptic coordinates
1 c 2 ( cosh 2 ξ cos 2 η ) 2 W ξ 2 + 2 W η 2 + k 2 W = 0 .
By substituting W ( ξ , η ) = F ( ξ ) G ( η ) and the identity
cosh 2 ξ cos 2 η = cosh 2 ξ cos 2 η 2 ,
we can split Equation (6) into two second order ordinary differential equations as follows:
d 2 G d η 2 + a 2 q cos 2 η G = 0 ,
d 2 F d ξ 2 a 2 q cosh 2 ξ F = 0 ,
where a is the separation constant and q = 0.25 k 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
W ( ξ , η ) = m = 0 C m 1 ( q ) Ce m ( ξ ; q ) ce m ( η ; q ) + C m 2 ( q ) Fe m ( ξ ; q ) ce m ( η ; q ) + m = 1 S m 1 ( q ) Se m ( ξ ; q ) se m ( η ; q ) + S m 2 ( q ) Ge m ( ξ ; q ) se m ( η ; q ) ,
where cem(η;q),sem(η;q) are the periodic Mathieu functions;
Cem(η;q),Sem(η;q) are the periodic modified Mathieu functions cem(η;q),sem(η;q), respectively;
Fem(η;q),Gem(η;q) are the non-periodic modified Mathieu functions cem(η;q),sem(η;q), respectively.
The coefficients C m 1 ( q ) , C m 2 ( q ) , S m 1 ( q ) and S m 2 ( 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].
According to the elliptic coordinate system ( ξ , η ) , where η represents the confocal hyperbolas defined on [ 0 , 2 π ) , the symmetric properties of the periodic Mathieu functions ce m ( η ; q ) and se m ( η ; q ) can be derived as follows:
  • ce 2 m ( η ; q ) is symmetrical about both x- and y-axes,
  • ce 2 m + 1 ( η ; q ) is symmetrical about x-axis but antisymmetrical about y-axis,
  • se 2 m + 1 ( η ; q ) is antisymmetrical about x-axis but symmetrical about y-axis,
  • se 2 m + 2 ( η ; q ) is antisymmetrical about both x- and y-axes,
for all integer m 0 . Moreover, the orthogonality of the periodic Mathieu functions is
0 2 π ce m ( η ; q ) ce n ( η ; q ) d η = 0 , for all integer m , n 0 and m n
0 2 π se m ( η ; q ) se n ( η ; q ) d η = 0 , for all integer m , n > 0 and m n
0 2 π ce m ( η ; q ) se n ( η ; q ) d η = 0 , for all integer m 0 and n > 0 .

3. Mathematical Modeling

3.1. Electroosmotic Force

The electroosmotic force is given by
F EOF = ρ e E ,
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:
2 ψ ( ξ , η ) = ρ e ε = 2 n 0 p + z ν ε sinh p + z ν ψ k B T ,
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
2 ψ = 1 c 2 ( cosh 2 ξ cos 2 η ) 2 ψ ξ 2 + 2 ψ η 2 = κ 2 ψ ,
where κ = ( 2 n 0 p + 2 z ν 2 ε 1 k B 1 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
F EOF = ε κ 2 E ψ .
Since Equation (16) is in the form of Equation (6), the general solution is obtained in the similar form as follows:
ψ ( ξ , η ) = m = 0 A m 1 ( q ) Ce m ( ξ ; q ) ce m ( η ; q ) + A m 2 ( q ) Fe m ( ξ ; q ) ce m ( η ; q ) + m = 1 A m 3 ( q ) Se m ( ξ ; q ) se m ( η ; q ) + A m 4 ( q ) Ge m ( ξ ; q ) se m ( η ; q ) ,
where A m 1 ( q ) , A m 2 ( q ) , A m 3 ( q ) and A m 4 ( q ) are constant with q = 0.25 κ 2 c 2 .
According to the symmetric property of an elliptic geometry, the potential ψ is supposed to be symmetrical about both x- and y-axes. Hence, we can express the boundary conditions of ψ in the elliptic coordinate system by the following equations:
ψ η ( ξ , 0 ) = 0 ,
ψ η ξ , π 2 = 0 ,
ψ ξ ( 0 , η ) = 0 .
The general solution shown in Equation (18) can be reduced, according to the symmetric property of the potential ψ , as follows:
ψ ( ξ , η ) = m = 0 A 2 m 1 ( q ) Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) + A 2 m 2 ( q ) Fe 2 m ( ξ ; q ) ce 2 m ( η ; q ) .
Since the partial derivative of the modified Mathieu function Fe 2 m ( ξ ; q ) with respect to ξ does not vanish when ξ = 0 [27], the constant A 2 m 2 ( q ) has to be zero in order to satisfy the boundary condition (21). Therefore, the general solution (22) is reduced to
ψ ( ξ , η ) = m = 0 A 2 m ( q ) Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) ,
where A 2 m ( q ) = A 2 m 1 ( 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
ψ ( ξ 0 , η ) = f ( η ) .
Now, we apply the boundary condition (24) to the general solution (23). Thus, the orthogonality of the Mathieu functions ce 2 m ( η ; q ) yields the constant A 2 m as
A 2 m = 0 2 π f ( η ) · ce 2 m ( η ; q ) d η 0 2 π Ce 2 m ( ξ 0 ; q ) ce 2 m 2 ( η ; q ) d η .

3.2. 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:
1 c 2 ( cosh 2 ξ cos 2 η ) 2 u ξ 2 + 2 u η 2 ρ μ u t = p z μ F EOF μ ,
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
1 c 2 ( cosh 2 ξ cos 2 η ) 2 u ξ 2 + 2 u η 2 ρ μ u t = p z μ + ε κ 2 E μ ψ .
To eliminate the time-derivative, we employ the Laplace transform L defined by
v ( ξ , η , s ) = L { u ( ξ , η , t ) } = 0 u ( ξ , η , t ) e s t d t ,
and the initial condition u ( ξ , η , 0 ) = 0 . Therefore, Equation (27) becomes
1 c 2 ( cosh 2 ξ cos 2 η ) 2 v ξ 2 + 2 v η 2 ρ s μ v = p z μ s + ε κ 2 E μ s ψ ,
which has the general solution written in the form of
v ( ξ , η , s ) = v c ( ξ , η , s ) + v p 1 ( ξ , η , s ) + v p 2 ( ξ , η , s ) ,
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
1 c 2 ( cosh 2 ξ cos 2 η ) 2 v c ξ 2 + 2 v c η 2 ρ s μ v c = 0 .
Since Equation (31) is in the form of Equation (6), the complementary solution v c is obtained in the similar form as follows:
v c = m = 0 B m 1 ( q s ) Ce m ( ξ ; q s ) ce m ( η ; q s ) + B m 2 ( q s ) Fe m ( ξ ; q s ) ce m ( η ; q s ) + m = 1 B m 3 ( q s ) Se m ( ξ ; q s ) se m ( η ; q s ) + B m 4 ( q s ) Ge m ( ξ ; q s ) se m ( η ; q s ) ,
where B m 1 ( q s ) , B m 2 ( q s ) , B m 3 ( q s ) , and B m 4 ( 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:
v p 1 = p z ρ s 2 .
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
v p 2 = m = 0 C 2 m Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) .
Since v p 2 is the particular solution of the corresponding non-homogeneous equation
1 c 2 ( cosh 2 ξ cos 2 η ) 2 v ξ 2 + 2 v η 2 ρ s μ v = ε κ 2 E μ s ψ ,
we have
2 m = 0 C 2 m Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) ρ s μ m = 0 C 2 m Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) = ε κ 2 E μ s m = 0 A 2 m Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) .
Moreover, substituting Equation (23) into Equation (16) and equating the term having the function A 2 m ( q ) with the same integer m, then we obtain
2 Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) = κ 2 Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) ,
for all integer m 0 . Observe that -4.6cm0cm
2 m = 0 C 2 m Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) = m = 0 C 2 m 2 Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) = m = 0 C 2 m κ 2 Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) ,
we can rewritten Equation (36) as
m = 0 C 2 m κ 2 Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) ρ s μ m = 0 C 2 m Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) = ε κ 2 E μ s m = 0 A 2 m ( q ) Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) .
By equating coefficients of the terms having the same order of the periodic Mathieu function ce 2 m ( η ; q ) , we can derive the constant C 2 m for the particular solution v p 2 from
κ 2 ρ s μ C 2 m = ε κ 2 E μ s A 2 m ( q ) or C 2 m = ε κ 2 E s ( μ κ 2 s ρ ) A 2 m ( q ) ,
for all integer m 0 .
Combining Equations (30), (32)–(34) yields the general solution of Equation (29) as follows:
v ( ξ , η , s ) = m = 0 B m 1 ( q s ) Ce m ( ξ ; q s ) ce m ( η ; q s ) + B m 2 ( q s ) Fe m ( ξ ; q s ) ce m ( η ; q s ) + m = 1 B m 3 ( q s ) Se m ( ξ ; q s ) se m ( η ; q s ) + B m 4 ( q s ) Ge m ( ξ ; q s ) se m ( η ; q s ) p z ρ s 2 + m = 0 C 2 m Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) .
Similarly to the potential ψ , we also suppose that v is symmetrical about both x- and y-axes. Therefore, the solution in Equation (41) is reduced to
v ( ξ , η , s ) = m = 1 B 2 m 1 ( q s ) Ce 2 m ( ξ ; q s ) ce 2 m ( η ; q s ) + m = 0 C 2 m Ce 2 m ( ξ ; q ) ce 2 m ( η ; q ) p z ρ s 2 .
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 2 m 1 can be determined using the orthogonality of the Mathieu function ce 2 m ( η ; q s ) by
B 2 m 1 = 0 2 π b ( ξ 0 , η , s ) ce 2 m ( η ; q s ) d η 0 2 π Ce 2 m ( ξ 0 ; q s ) ce 2 m 2 ( η ; q s ) d η ,
for all integer m 0 , where
b ( ξ , η , s ) = p z ρ s 2 n = 0 C 2 n Ce 2 n ( ξ ; q ) ce 2 n ( η ; q ) .
The fluid velocity u can be obtained by employing the inverse Laplace transform as follows:
u ( ξ , η , t ) = L 1 { v ( ξ , η , s ) } = 1 2 π i lim T γ i T γ + i T e s t v ( ξ , η , s ) d s ,
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 η ) 2 u ξ 2 + 2 u η 2 = p z μ + ε E μ 2 ψ ,
where u is now the function depending only on ξ and η . Let
u ( ξ , η ) = w ( ξ , η ) + ε E μ ψ ( ξ , η ) + c 2 p z ( cosh 2 ξ + cos 2 η ) 8 μ .
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 η ) 2 w ξ 2 + 2 w η 2 = 0 .
Using the method of separation of variables, we obtain
w ( ξ , η ) = D 0 + m = 1 D 2 m cosh ( 2 m ξ ) cos ( 2 m η ) ,
where D 2 m is constant. Substituting the general solution w ( ξ , η ) into Equation (48) yields
u ( ξ , η ) = D 0 + m = 1 D 2 m cosh ( 2 m ξ ) cos ( 2 m η ) + ε E μ ψ + c 2 p z 8 μ ( cosh 2 ξ + cos 2 η ) ,
where the constant D 2 m can be derived from the no-slip condition u ( ξ 0 , η , t ) = 0 with the orthogonality of the trigonometric functions as
D 2 m = 0 2 π [ ε E f ( η ) μ 1 + 0.125 c 2 p z μ 1 [ cosh ( 2 ξ 0 ) + cos ( 2 η ) ] ] cos ( 2 m η ) d η 0 2 π cosh ( 2 m ξ 0 ) cos 2 ( 2 m η ) d η .

4. 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 2 m ( z ; q ) , Ce 2 m ( z ; q ) , ce 2 m ( z ; q ) and Ce 2 m ( 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.
Figure 3 illustrates the velocity profiles of transient flow along both x- and y-axes. The results show that the flow develops only from the area adjacent to the wall at the beginning time t = 0.1 and t = 0.25 m s represented by the red and green lines, respectively. This is due to the effect of the moving ions in the fluid within the electrical double layers. However, at time t = 0.5 m s , represented by the cyan line, the flow at the center of the channel starts developing. The flow becomes fully developed at time t = 2.5 m s , which represented by the gray line.
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:
c = e A 1 e 2 ,
ξ 0 = ln 1 + 1 e 2 e .
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 x- and y-axes are shown in Figure 5a,b, respectively. The results show that there is no significant difference for the velocity profiles along both x- and 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 x- and 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 = e C 0 2 π 1 e 2 sin η d η .
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 x- and 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 x- and 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
Q ( t ) = 0 2 π 0 ξ 0 c 2 ( cosh 2 ξ cos 2 η ) u ( ξ , η , t ) d ξ d η .
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
δ i = Q i Q 0 Q 0 ,
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.

5. 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.

Author Contributions

Both authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

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 2 m ( z ; q ) , Ce 2 m ( z ; q ) , ce 2 m ( z ; q ) and Ce 2 m ( 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 2 m ( z ; q ) , Ce 2 m ( z ; q ) , ce 2 m ( z ; q ) and Ce 2 m ( z ; q ) can be expressed by
ce 2 m ( z ; q ) = r = 0 A 2 r ( 2 m ) cos 2 r z ,
Ce 2 m ( z ; q ) = r = 0 A 2 r ( 2 m ) cosh 2 r z ,
ce 2 m ( z ; q ) = ( 1 ) n r = 0 ( 1 ) r A 2 r ( 2 m ) cos 2 r z ,
Ce 2 m ( z ; q ) = ( 1 ) n r = 0 ( 1 ) r A 2 r ( 2 m ) cosh 2 r z .
To find the coefficients A 2 r ( 2 m ) , substituting Equation (A1) into Equation (8) and equating the coefficients of the term cos 2 r z for each r = 0 , 1 , 2 , yield the following recurrence relations:
a 2 m A 0 ( 2 m ) q A 2 ( 2 m ) = 0 ,
[ a 2 m 4 ] A 2 ( 2 m ) q A 4 ( 2 m ) + 2 A 0 ( 2 m ) = 0 ,
[ a 2 m 4 r 2 ] A 2 r ( 2 m ) q A 2 r + 2 ( 2 m ) + A 2 r 2 ( 2 m ) = 0 , r 2 ,
where a 2 m is the separation constant a corresponding to ce 2 m . Dividing Equations (A5) and (A6) by A 0 ( 2 m ) and Equation (A7) by A 2 r 2 ( 2 m ) , we obtain
a 2 m q v 0 ( 2 m ) = 0 ,
[ a 2 m 4 ] v 0 ( 2 m ) q v 0 ( 2 m ) v 2 ( 2 m ) + 2 = 0 ,
[ a 2 m 4 r 2 ] v 2 r 2 ( 2 m ) q v 2 r 2 ( 2 m ) v 2 r ( 2 m ) + 1 = 0 , r 2 ,
where
v 2 r ( 2 m ) = A 2 r + 2 ( 2 m ) A 2 r ( 2 m ) .
Rearranging Equations (A8)–(A10), we have
v 0 ( 2 m ) = a 2 m q ,
v 0 ( 2 m ) = 1 2 q 1 1 4 a 2 m q v 2 ( 2 m ) ,
v 2 r 2 ( 2 m ) = 1 4 r 2 q 1 1 4 r 2 a 2 m q v 2 r ( 2 m ) , r 2 .
Solving the recurrence relations (A13) and (A14) yields
v 0 = 1 2 q 1 1 4 a 2 m 1 64 q 2 1 1 16 a 2 m 1 576 q 2 1 1 36 a 2 m 1 2304 q 2 1 1 64 a 2 m .
Substituting v 0 from Equation (A12) into Equation (A15), we can calculate a 2 m by solving the following equation:
a 2 m = 1 2 q 2 1 1 4 a 2 m 1 64 q 2 1 1 16 a 2 m 1 576 q 2 1 1 36 a 2 m 1 2304 q 2 1 1 64 a 2 m .
Since v 2 r ( 2 m ) 0 as r 0 , we suppose v 2 r ( 2 m ) = 0 for all r R provided R is large enough. Then v 2 R 2 ( 2 m ) , v 2 R 4 ( 2 m ) , , v 0 ( 2 m ) can be obtained by solving the recurrence relations (A12)–(A14). The normalization of ce 2 m ( z ; q ) by
1 π 0 2 π ce 2 m 2 ( z ; q ) d z = 1 ,
yields
1 = 2 A 0 ( 2 m ) 2 + A 2 ( 2 m ) 2 + A 4 ( 2 m ) 2 + A 6 ( 2 m ) 2 + ,
or
1 A 0 ( 2 m ) = 2 A 0 ( 2 m ) A 0 ( 2 m ) 2 + A 2 ( 2 m ) A 0 ( 2 m ) 2 + A 4 ( 2 m ) A 0 ( 2 m ) 2 + A 6 ( 2 m ) A 0 ( 2 m ) 2 + = 2 + v 0 ( 2 m ) 2 + v 0 ( 2 m ) v 2 ( 2 m ) 2 + v 0 ( 2 m ) v 2 ( 2 m ) v 4 ( 2 m ) 2 + .
Therefore, A 0 ( 2 m ) , A 2 ( 2 m ) , , A 2 R ( 2 m ) can be calculated using Equations (A11) and (A19) and the values of v 0 ( 2 m ) , v 2 ( 2 m ) , , v 2 R 2 ( 2 m ) .
It should be noted that the values of A 2 r ( 2 m ) by the above procedure can be use for all of the functions ce 2 m ( z ; q ) , Ce 2 m ( z ; q ) , ce 2 m ( z ; q ) and Ce 2 m ( z ; q ) .

References

  1. Beskok, A.; Korniadakis, G.E. Report: A model for flows in channels, pipes, and ducts at micro and nano scales. Nanoscale Microscale Thermophys. Eng. 1999, 3, 43–77. [Google Scholar] [CrossRef]
  2. Araki, T.; Kim, M.S.; Suzuki, K. An experimental investigation of gaseous flow characteristics in microchannels. Nanoscale Microscale Thermophys. Eng. 2002, 6, 117–130. [Google Scholar] [CrossRef]
  3. Lee, H.B.; Yeo, I.W.; Lee, K.K. Water flow and slip on NAPL-wetted surfaces of a parallel-walled fracture. Geophys. Res. Lett. 2007, 34, L19401. [Google Scholar] [CrossRef]
  4. Duan, Z. Slip flow in doubly connected microchannels. Int. J. Therm. Sci. 2012, 58, 45–51. [Google Scholar] [CrossRef]
  5. Mao, Z.; Yoshida, K.; Kim, J.W. A droplet-generator-on-a-chip actuated by ECF (electro-conjugate fluid) micropumps. Microfluid. Nanofluid. 2019, 23, 130. [Google Scholar] [CrossRef]
  6. Xie, J.; Hu, G.H. Computational modelling of membrane gating in capsule translocation through microchannel with variable section. Microfluid. Nanofluid. 2021, 25, 17. [Google Scholar] [CrossRef]
  7. Hunter, R.J. Zeta Potential in Colloid Science: Principles and Applications, 1st ed.; Academic: London, UK, 1981; pp. 1–10. [Google Scholar] [CrossRef]
  8. Li, D. Electrokinetics in Microfluidics, 1st ed.; Elsevier: London, UK, 2004; pp. 8–28. [Google Scholar] [CrossRef]
  9. Kim, M.; Beskok, A.; Kihm, K. Electro-osmosis-driven micro-channel flows: A comparative study of microscopic particle image velocimetry measurements and numerical simulations. Exp. Fluids 2002, 33, 170–180. [Google Scholar] [CrossRef]
  10. Beddiar, K.; Fen-Chong, T.; Dupas, A.; Berthaud, Y.; Dangla, P. Role of pH in electro-osmosis: Experimental study on NaCl–water saturated kaolinite. Transp. Porous Med. 2005, 61, 93–107. [Google Scholar] [CrossRef]
  11. Bandopadhyay, A.; Chakraborty, S. Electrokinetically induced alterations in dynamic response of viscoelastic fluids in narrow confinements. Phys. Rev. E 2012, 85, 056302. [Google Scholar] [CrossRef]
  12. Chakraborty, J.; Ray, S. Chakraborty, S. Role of streaming potential on pulsating mass flow rate control in combined electroosmotic and pressure-driven microfluidic devices. Electrophoresis 2012, 33, 419–425. [Google Scholar] [CrossRef]
  13. Chinyoka, T.; Makinde, O.D. Analysis of non-Newtonian flow with reacting species in a channel filled with a saturated porous medium. J. Pet. Sci. Eng. 2014, 121, 1–8. [Google Scholar] [CrossRef]
  14. Shen, Y.; Shi, W.; Li, S.; Yang, L.; Feng, J.; Gao, M. Study on the electro-osmosis characteristics of soft clay from Taizhou with various saline solutions. Adv. Civ. Eng. 2020, 2020, 6752565. [Google Scholar] [CrossRef]
  15. Arulanandam, S.; Li, D. Liquid transport in rectangular microchannels by electroosmotic pumping. Colloids Surf. A 2000, 161, 89–102. [Google Scholar] [CrossRef]
  16. Reshadi, M.; Saidi, M.H. Firoozabadi, B.; Saidi, M.S. Electrokinetic and aspect ratio effects on secondary flow of viscoelastic fluids in rectangular microchannels. Microfluid. Nanofluid. 2016, 20, 117. [Google Scholar] [CrossRef]
  17. Siddiqui, A.A.; Lakhtakia, A. Debye-Hückel solution for steady electro-osmotic flow of micropolar fluid in cylindrical microcapillary. Appl. Math. Mech. Engl. Ed. 2013, 34, 1305–1326. [Google Scholar] [CrossRef] [Green Version]
  18. Tseng, S.; Tai, Y.H.; Hsu, J.P. Ionic current in a pH-regulated nanochannel filled with multiple ionic species. Microfluid. Nanofluid. 2014, 17, 933–941. [Google Scholar] [CrossRef]
  19. Tsao, H.K. Electroosmotic Flow through an Annulus. J. Colloid Interface Sci. 2000, 225, 247–250. [Google Scholar] [CrossRef] [PubMed]
  20. Na, R.; Jian, Y.; Chang, L.; Su, J.; Liu, Q. Transient electro-osmotic and pressure driven flows through a microannulus. Open J. Fluid Dyn. 2013, 3, 50–56. [Google Scholar] [CrossRef] [Green Version]
  21. Goswami, P.; Chakraborty, S. Semi-analytical solutions for electroosmotic flows with interfacial slip in microchannels of complex cross-sectional shapes. Microfluid. Nanofluid. 2011, 11, 255–267. [Google Scholar] [CrossRef]
  22. Chuchard, P.; Orankitjaroen, S.; Wiwatanapataphee, B. Study of pulsatile pressure-driven electroosmotic flows through an elliptic cylindrical microchannel with the Navier slip condition. Adv. Differ. Equ. 2017, 160. [Google Scholar] [CrossRef] [Green Version]
  23. Yun, J.H.; Chun, M.S.; Jung, H.W. The geometry effect on steady electrokinetic flows in curved rectangular microchannels. Phys. Fluids 2010, 22, 052004. [Google Scholar] [CrossRef] [Green Version]
  24. Vocale, P.; Geri, M.; Cattani, L.; Morini, G.L.; Spiga, M. Numerical analysis of electro-osmotic flows through elliptic microchannels. La Houille Blanche 2013, 3, 42–49. [Google Scholar] [CrossRef]
  25. Srinivas, B. Electroosmotic flow of a power law fluid in an elliptic microchannel. Colloids Surf. A Physicochem. Eng. Asp. 2016, 492, 144–151. [Google Scholar] [CrossRef]
  26. Parida, M.; Padhy, S. Electro-osmotic flow of a third-grade fluid past a channel having stretching walls. Nonlinear Eng. 2019, 8, 56–64. [Google Scholar] [CrossRef]
  27. McLachlan, N.W. Theory and Application of Mathieu Functions, 1st ed.; Clarendon: Oxford, UK, 1947; pp. 1–394. [Google Scholar]
  28. Liu, B.T.; Tseng, S.; Hsu, J.P. Effect of eccentricity on the electroosmotic flow in an elliptic channel. J. Colloid Interface Sci. 2015, 460, 81–86. [Google Scholar] [CrossRef]
  29. Cohen, Y.; Metzner, A.B. Apparent slip flow of polymer solutions. J. Rheol. 1985, 29, 67–102. [Google Scholar] [CrossRef]
  30. Tretheway, D.C.; Meinhart, C.D. Apparent fluid slip at hydrophobic microchannel walls. Phys. Fluids 2002, 14, 9–12. [Google Scholar] [CrossRef] [Green Version]
  31. Mathieu, E.L. Mémoire sur le mouvement vibratoire d’une membrane de forme elliptique. J. Math. Pures Appl. 1868, 13, 137–203. [Google Scholar]
Figure 1. The electrical double layer (red highlighted area) and electroosmotic flow: (a) Cross-section view; (b) Along the channel length view.
Figure 1. The electrical double layer (red highlighted area) and electroosmotic flow: (a) Cross-section view; (b) Along the channel length view.
Computation 09 00027 g001
Figure 2. Elliptic coordinate system ( ξ , η ) : the blue line represents the coordinate ξ and the red line represents the coordinate η .
Figure 2. Elliptic coordinate system ( ξ , η ) : the blue line represents the coordinate ξ and the red line represents the coordinate η .
Computation 09 00027 g002
Figure 3. Variation of the velocity profiles of transient electroosmotic flow on time: (a) Along x-axis; (b) Along y-axis.
Figure 3. Variation of the velocity profiles of transient electroosmotic flow on time: (a) Along x-axis; (b) Along y-axis.
Computation 09 00027 g003
Figure 4. Variation of the elliptic cross-sections of the channel on various eccentricities with: (a) Fixed area of the channel cross-sections; (b) Fixed circumference of the channel cross-sections.
Figure 4. Variation of the elliptic cross-sections of the channel on various eccentricities with: (a) Fixed area of the channel cross-sections; (b) Fixed circumference of the channel cross-sections.
Computation 09 00027 g004
Figure 5. Variation of the velocity profiles of steady electroosmotic flow, when the area of the channel cross-sections is fixed, with various eccentricities: (a) Along x-axis; (b) Along y-axis.
Figure 5. Variation of the velocity profiles of steady electroosmotic flow, when the area of the channel cross-sections is fixed, with various eccentricities: (a) Along x-axis; (b) Along y-axis.
Computation 09 00027 g005
Figure 6. Variation of the velocity profiles of steady electroosmotic flow, when the circumference of the channel cross-sections is fixed, with various eccentricities: (a) Along x-axis; (b) Along y-axis.
Figure 6. Variation of the velocity profiles of steady electroosmotic flow, when the circumference of the channel cross-sections is fixed, with various eccentricities: (a) Along x-axis; (b) Along y-axis.
Computation 09 00027 g006
Figure 7. Variation of the flow rates of steady electroosmotic flow on various eccentricities with: (a) Fixed area of the channel cross-sections; (b) Fixed circumference of the channel cross-sections.
Figure 7. Variation of the flow rates of steady electroosmotic flow on various eccentricities with: (a) Fixed area of the channel cross-sections; (b) Fixed circumference of the channel cross-sections.
Computation 09 00027 g007
Figure 8. Variation of the relative errors with three different values of external electric field on various eccentricities with: (a) Fixed area of the channel cross-sections; (b) Fixed circumference of the channel cross-sections.
Figure 8. Variation of the relative errors with three different values of external electric field on various eccentricities with: (a) Fixed area of the channel cross-sections; (b) Fixed circumference of the channel cross-sections.
Computation 09 00027 g008
Table 1. Table of parameters used in the transient flow investigation.
Table 1. Table of parameters used in the transient flow investigation.
NameSymbol ValueSI Unit
Focal lengthc 4.50 × 10 5 m
Eccentricitye 0.60 -
Fluid density ρ 1.00 × 10 3 k g m 3
Fluid viscosity μ 9.00 × 10 4 Pa s
Fluid permittivity ε 6.95 × 10 10 F m 1
Pressure gradient in z-axis p z 2.00 Pa m 1
Reciprocal of EDL thickness κ 8.00 × 10 4 m 1
Surface potential f ( η ) 2.49 × 10 3 V
External electric fieldE 5.00 × 10 2 V m 1
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Numpanviwat, N.; Chuchard, P. Transient Pressure-Driven Electroosmotic Flow through Elliptic Cross-Sectional Microchannels with Various Eccentricities. Computation 2021, 9, 27. https://0-doi-org.brum.beds.ac.uk/10.3390/computation9030027

AMA Style

Numpanviwat N, Chuchard P. Transient Pressure-Driven Electroosmotic Flow through Elliptic Cross-Sectional Microchannels with Various Eccentricities. Computation. 2021; 9(3):27. https://0-doi-org.brum.beds.ac.uk/10.3390/computation9030027

Chicago/Turabian Style

Numpanviwat, Nattakarn, and Pearanat Chuchard. 2021. "Transient Pressure-Driven Electroosmotic Flow through Elliptic Cross-Sectional Microchannels with Various Eccentricities" Computation 9, no. 3: 27. https://0-doi-org.brum.beds.ac.uk/10.3390/computation9030027

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop