Analytical solutions for electroosmotic core–annular flow of power-law fluids

1 Introduction

Micro- and nanofluidic devices have been increasingly used for various biomedical, chemical, and environmental applications, such as disease diagnostics, drug delivery, and particle sorting, where an accurate control of fluid flow is often necessary for precise sample manipulation (Song et al., 2023; Shi et al., 2024; Ohshima, 2024; Khanjani et al., 2025; David and Cai, 2026). Electroosmosis is an effective method for fluid transport in micro- and nanochannels (Li, 2004; Chang and Yeo, 2009; Sasmal et al., 2026; Yu et al., 2026). It is the fluid motion induced by the movement of excess counter-ions in the electric double layer (EDL) when an external electric field is applied across the channel length (Hunter, 1981; Masliyah and Bhattacharjee, 2006). Compared with traditional pressure-driven flow, electroosmotic flow possesses a nearly plug-like velocity profile in microchannels, which can substantially reduce the hydrodynamic resistance and sample dispersion (Whitesides and Stroock, 2001; Ghosal, 2004; Ghosal, 2006). To date, the work on electroosmosis, and more generally, electrokinetic phenomena, has been mainly focused on Newtonian fluids (Zhao and Yang, 2012; Xuan, 2019; Alizadeh et al., 2021; Sun et al., 2023), although micro- and nanofluidic devices applications frequently involve complex fluids (e.g., blood, colloidal suspension, and polymer solution) that exhibit viscoelastic behavior, shear-dependent behavior (shear thinning or shear thickening), or both (Lu et al., 2017; Browne et al., 2020; Kang et al., 2025). These non-Newtonian rheological properties have been reported to strongly affect electrokinetic applications, such as fluid pumping (Berli, 2010; Roy and Dhar, 2026), sample mixing (Kumar et al., 2025; Islam et al., 2026; Sasmal et al., 2026), energy conversion (Xie and Jian, 2022; Paul et al., 2020; Roy and Poddar, 2026; Fan and Liu, 2026), and particle separation (Khair et al., 2012; Bentor and Xuan, 2024; Bentor et al., 2026). A fundamental understanding of electroosmotic flow in non-Newtonian fluids is critical for the design and optimization of these micro- and nanochannels devices.

As reviewed by Zhao and Yang (2013a), a considerable number of theoretical studies have been conducted to investigate the effects of fluid rheology on the velocity of electroosmotic flow, either alone or combined with pressure-driven flow, in non-Newtonian fluids. The fluid shear-thinning or shear-thickening effect is characterized by the power-law or Carreau model (Das and Chakraborty, 2006; Chakraborty, 2007; Zhao et al., 2008; Tang et al., 2009; Zhao and Yang, 2011a; Yang, 2011b; Babaie et al., 2011; Zhao and Yang, 2013b), while the fluid elasticity effect is characterized by the upper convected Maxwell or Oldroyd-B model (Zhao and Yang, 2009; Zhao et al., 2013; Mahapatra and Bandopadhyay, 2021). The combination of both effects is described using the Phan–Thien–Tanner (PTT), Giesekus, or FENE-P model (Afonso et al., 2009; Choi et al., 2012; Ferrás et al., 2016; Natu and Ghosh, 2021). It is noteworthy that a generalized exponential-structure rheological constitutive equation was recently proposed (Herrera-Valencia et al., 2024) to predict complex non-Newtonian behaviors, such as shear thinning, shear thickening, thixotropy, yield stress, shear banding, elasticity, and normal stress differences (Ramírez-Torres et al., 2024). This model has been utilized to study the electroosmotic flow of viscoelastic fluids under both a pulsating electric field (Herrera-Valencia et al., 2024) and a combined pulsating electric field and pressure gradient (Herrera-Valencia et al., 2025). All the above studies considered the electroosmotic flow of single non-Newtonian fluids. However, recent studies have indicated the inadequacy of these single-fluid models in predicting the experimentally measured electroosmotic velocity of viscoelastic or shear-thinning fluids (Chang and Tsao, 2007; Berli and Olivares, 2008; Leung et al., 2025). Several experiments have suggested that interactions between polymer molecules and channel walls lead to the formation of a polymer depletion or adsorption layer, altering the local fluid rheology (de Gennes, 1981; Aubert and Tirrell, 1982; Ma and Graham, 2005; Park et al., 2019). This near-wall rheological heterogeneity is important in electroosmotic flow because the electrical body force is concentrated primarily within the EDL adjacent to the channel wall (Hunter, 1981; Li, 2004). Therefore, even when the wall-modified layer is thin relative to the channel dimension, it may still substantially overlap with the electroosmotic forcing region and thereby exert a strong influence on the velocity of the bulk fluid.

The concept of a polymer depletion layer (PDL) has been proposed by multiple groups to qualitatively explain the discrepancy between experimentally and theoretically obtained electroosmotic velocities (Chang and Tsao, 2007; Berli and Olivares, 2008). In particular, Mukherjee et al. (2017) developed a PTT model-based theory that incorporates a Newtonian PDL, yielding a closed-form solution for the electroosmotic velocity of viscoelastic fluids in a parallel-plate microchannel. The predictions of their two-layer model were demonstrated to agree better with the experimental velocity of polyethylene oxide solutions reported by Huang et al. (2016) than those of the single-fluid model. Moschopoulos et al. (2020) extended this two-layer model by considering the PDL and ionic steric effects on the electroosmotic flow of viscoelastic fluids. Their revised numerical model provided better agreement with the experimental results of Huang et al. (2016) than the analytical model of Mukherjee et al. (2017). Recently, Leung et al. (2025) developed a similar two-layer numerical model of power-law fluids to understand the experimentally measured electroosmotic velocity of shear-thinning xanthan gum solutions reported by Bentor et al. (2024). Their model could predict with good accuracy the nonlinear electroosmotic velocity in xanthan gum solutions with varying polymer and buffer concentrations across the tested electric fields. More recently, Leung et al. (2026) presented a two-layer theoretical model and obtained an approximate analytical solution for the electroosmotic velocity of power-law fluids in a parallel-plate microchannel. The authors also derived an analytical formula that could accurately predict the nonlinear electric field index of the experimental electroosmotic velocity (Bentor et al., 2024). In addition, several other theoretical papers have used the two-layer model to investigate the electroosmotic pumping and heat transfer of non-Newtonian fluids in slit microchannels (Sousa et al., 2011; Afonso et al., 2013; Huang et al., 2014; Mahapatra and Bandopadhyay, 2022; Roychowdhury et al., 2022; Deng and Xiao, 2022).

In contrast, there are only a limited number of studies on the use of a two-layer model for the electroosmotic flow of non-Newtonian fluids in cylindrical microchannels. Olivares et al. (2009) derived an analytical formula for the electroosmotic velocity of polymer solutions in a round capillary, in which the polymer depletion or adsorption layer was described using a power-law model with parameters different from those of the bulk. The near-wall fluid layer was assumed to be much thicker than the EDL in their model, whose predictions could explain the experimental data in viscoelastic carboxymethyl cellulose solutions with fitted rheological properties. Deng et al. (2021) revised this two-layer model without imposing the thin EDL requirement and obtained semi-analytical solutions for the electroosmotic flow (Deng, 2022) and the combined electroosmotic and pressure-driven flow (Deng et al., 2021) of power-law fluids.

This work builds upon our earlier two-layer theoretical model for a slit microchannel (Leung et al., 2026) to obtain analytical solutions for the generalized electroosmotic core–annular flow of power-law fluids in a cylindrical microchannel. The choice of the power-law constitutive equation over the more accurate Carreau (Bird et al., 1987; Raihan et al., 2021) and Carreau–Yasuda (Yasuda et al., 1981; Jagdale et al., 2020) models is due to its mathematical simplicity, which has been found to be sufficient in our recent paper (Leung et al., 2025) to predict the fluid shear-thinning effect on the experimentally measured electroosmotic velocity of xanthan gum solutions (Bentor et al., 2024). Similar to the models of Deng et al. (2021) and Deng (2022), our model solves the EDL potential field from the Poisson–Boltzmann equation under the Debye–Hückel approximation (Hunter, 1981; Li, 2004) and the electroosmotic velocity field from the Cauchy momentum equation (Olivares et al., 2009; Zhao and Yang, 2013b). It neglects the potential difference at the interface of the two layers and assumes a uniform permittivity and ionic concentration across the core and annular fluids. These simplifications enable us to derive approximate and simplified analytical solutions for the electroosmotic velocity and flow rate under appropriate assumptions, avoiding the need for numerical calculation of the integral equations used in previous studies (Deng et al., 2021; Deng, 2022; Deng and Xiao, 2022). Moreover, either layer of the fluid in our model can have arbitrary power-law model parameters, ranging from shear-thinning to Newtonian and shear-thickening behavior. Therefore, the presented analytical solutions for the electroosmotic velocity and flow rate will be useful for studying the wall-induced polymer depletion or adsorption effect, along with electroosmotic pumping and electrokinetic energy conversion in generalized core–annular flows of power-law fluids.

2 Theory

2.1 EDL potential and electroosmotic velocity

We consider the electroosmotic core–annular flow of power-law fluids in a cylindrical microchannel of radius R subjected to an axial electric field, E, as illustrated in Figure 1. The core fluid in the channel center has a flow consistency index, m I, and a flow behavior index (or power-law index for convenience), n I, in the power-law model with the local velocity, u I. The annular fluid of radial thickness δ along the channel wall has m II, n II, and u II, accordingly. Note that n = 1 corresponds to a Newtonian fluid, while n < 1 and n > 1 indicate shear-thinning and shear-thickening fluids, respectively (Bird et al., 1987). We assume that the core and annular fluids are incompressible and immiscible. The channel wall has a uniform zeta potential, ζ, meeting the small surface potential requirement, ζ < 25mV, for the Debye–Hückel approximation (Hunter, 1981; Li, 2004), which enables the analytical solution for the EDL potential field and, in turn, the analytical solution for the electroosmotic velocity field. The flow is treated as steady and fully developed in an infinitely long cylindrical microchannel with only the axial (i.e., z-direction) velocity being non-zero and dependent solely on the radial coordinate r.

The electroosmotic velocities in the core (reflected by the subscript I) and annular (reflected by the subscript II) fluids are each obtained by solving the Cauchy momentum equation (Equations 1a, 1b) (Olivares et al., 2009; Zhao and Yang, 2013b): where ρ e is the local free-charge density and associated with the EDL potential. We assume that the core and annular fluids are prepared in the same buffer solution, wherein the ions follow the Boltzmann distribution regardless of fluid rheology in straight, uniform micro- and nanochannels (Zhao and Yang, 2013a). Thus, the EDL potential ψ can be determined from the linearized Poisson–Boltzmann equation for small zeta potentials (Equations 2a, 2b) (Hunter, 1981; Masliyah and Bhattacharjee, 2006): where ε is the fluid permittivity, which is assumed to be uniform throughout the channel and independent of the fluid rheology for dilute polymer solutions, and κ is the inverse of the so-called Debye length (Equation 3) (Li, 2004; Chang and Yeo, 2009):

1 r d d r r m I − d u I d r n I − 1 d u I d r + ρ e E = 0 , (1a)

1 r d d r r m II − d u II d r n II − 1 d u II d r + ρ e E = 0 , (1b)

ε r d d r r d ψ d r = − ρ e , (2a)

ρ e = − ε κ 2 ψ , (2b)

κ = ∑ c i z i 2 e 2 ε k B T . (3)

In this definition, c i and z i are the concentration and valence of ionic species i, respectively, e is the elementary charge, T is the fluid temperature, and k B is the Boltzmann constant. With the boundary conditions (Equations 4a, 4b), the EDL potential can be obtained by solving Equations 2a,b, , as given in Equation 5:where I 0 x is the zeroth-order modified Bessel function of the first kind.

ψ r = R = ζ , (4a)

d ψ d r r = 0 = 0 , (4b)

ψ r = ζ I 0 κ r I 0 κ R , (5)

The Cauchy momentum Equations 1a,b, are subject to the no-slip boundary condition at the channel wall and the symmetry condition in the center, as shown in Equations 6a, 6b:

d u I d r r = 0 = 0 , (6a)

u II r = R = 0 . (6b)

Additionally, the velocity and shear stress are required to be continuous at the interface of immiscible fluids, as shown in Equations 6c, 6d:

u I r = R − δ = u II r = R − δ , (6c)

m I − d u I d r n I − 1 d u I d r r = R − δ = m II − d u II d r n II − 1 d u II d r r = R − δ . (6d)

The final expressions for the electroosmotic velocity profiles in the core and annular fluids are presented in Equations 7a, 7b, respectively:

u I r = u II R − δ + κ 1 − n I n I − ε ζ E m I 1 n I ∫ κ r κ R − δ I 1 κ r I 0 κ R 1 n I d κ r for 0 ≤ r ≤ R − δ , (7a)

u II r = κ 1 − n II n II − ε ζ E m II 1 n II ∫ κ r κ R I 1 κ r I 0 κ R 1 n II d κ r for R − δ ≤ r ≤ R , (7b)

where I 1 κ r is the first-order modified Bessel function of the first kind. If the two fluids are represented by the same power-law model, with n I = n II = n and m I = m II = m being their rheological parameters, we can obtain the following electroosmotic velocity of power-law fluids (Equation 8) in cylindrical microchannels by combining Equations 7a,b, :

u PL r = κ 1 − n n − ε ζ E m 1 n ∫ κ r κ R I 1 κ r I 0 κ R 1 n d κ r . (8)

This expression is identical to that derived by Zhao and Yang (2013b) and has a closed-form expression only when n is equal to 1 or 1 / 2. Specifically, for n = 1 and m = η, Equation 8 further reduces to Equation 9:which corresponds to the well-known electroosmotic velocity of Newtonian fluids in cylindrical microchannels under the Debye–Hückel approximation (Rice and Whitehead, 1965).

u N r = − ε ζ E η 1 − I 0 κ r I 0 κ R , (9)

2.2 Approximate analytical solutions

To obtain analytical solutions for the electroosmotic core–annular flow of power-law fluids with arbitrary rheological parameters, we employ the following approximation to the Bessel function (Abramowitz and Stegun, 1965; Zhao and Yang, 2013b), I 1 κ r, in Equations 7a,b, :

I 1 κ r = κ r 2 , 0 ≤ κ r < 1 e κ r 2 π κ r 1 − 3 8 κ r , κ r ≥ 1 . (10)

Substituting Equation 10 into Equations 7a,b, yields the analytical solutions for the electroosmotic velocity in the core and annular fluids for the following three different cases. It is important to note that the piecewise approximation of the Bessel function in Equation 10 yields an increasing error as κ r approaches unity (Abramowitz and Stegun, 1965; Zhao and Yang, 2013b).

2.2.1 Cases 1: κ R − δ ≥ 1 and κ R > 1

This case corresponds to microchannels, wherein the electroosmotic velocity of the annular fluid for κ R − δ ≤ κ r ≤ κ R can be obtained from Equation 7b, as shown in Equation 11a:

u II r = U PL , II I 0 1 n II κ R { − 1 2 π n II 1 2 n II [ Γ 1 − 1 2 n II , − κ R n II + 3 8 n II 2 Γ − 1 2 n II , − κ R n II − Γ 1 − 1 2 n II , − κ r n II − 3 8 n II 2 Γ − 1 2 n II , − κ r n II ] } . (11a)

The electroosmotic velocity of the core fluid for 0 ≤ κ r < 1 can be obtained from Equation 7a, as shown in Equation 11b:while that for 1 ≤ κ r < κ R − δ is given as follows (Equation 11c):

u I r = u II R − δ + U PL , I I 0 1 n I κ R { 1 − κ r 1 + n I n I n I + 1 2 1 n I + − 1 2 π n I 1 2 n I [ Γ 1 − 1 2 n I , − κ R − δ n I + 3 8 n I 2 Γ − 1 2 n I , − κ R − δ n I − Γ 1 − 1 2 n I , − 1 n I − 3 8 n I 2 Γ − 1 2 n I , − 1 n I ] } , (11b)

u I r = u II R − δ + U PL , I I 0 1 n I κ R { − 1 2 π n I 1 2 n I [ Γ 1 − 1 2 n I , − κ R − δ n I + 3 8 n I 2 Γ − 1 2 n I , − κ R − δ n I − Γ 1 − 1 2 n I , − κ r n I − 3 8 n I 2 Γ − 1 2 n I , − κ r n I ] } . (11c)

In the above, Γ a , b denotes the upper incomplete gamma function (Abramowitz and Stegun, 1965) and U PL is the generalized Smoluchowski velocity for power-law fluids (Zhao et al., 2008) given as follows:

U PL = n κ 1 − n n − ε ζ E m 1 n . (12)

Note that for n = 1 and m = η, Equation 12 further reduces to the classical Helmholtz–Smoluchowski equation for Newtonian fluids (Hunter, 1981; Masliyah and Bhattacharjee, 2006):

U HS = − ε ζ E η . (13)

2.2.2 Cases 2: κ R − δ < 1 and κ R ≤ 1

This case corresponds to nanochannels, wherein the electroosmotic velocity of the annular fluid for κ R − δ ≤ κ r < κ R can be obtained from Equation 7b, as shown in Equation 14a:

u II r = U PL , II I 0 1 n II κ R κ R 1 + n II n II − κ r 1 + n II n II n II + 1 2 1 n II . (14a)

The electroosmotic velocity of the core fluid for 0 ≤ κ r < κ R − δ can be obtained from Equation 7a, as shown in Equation 14b:

u I r = u II R − δ + U PL , I I 0 1 n I κ R κ R − δ 1 + n I n I − κ r 1 + n I n I n I + 1 2 1 n I . (14b)

2.2.3 Cases 3: κ R − δ < 1 and κ R > 1

This case corresponds to mesoscopic or submicron channels, wherein the electroosmotic velocity of the annular fluid for κ R − δ ≤ κ r < 1 can be obtained from Equation 7b, as shown in Equation 15a:while that for 1 ≤ κ r < κ R is given as Equation 15b:

u II r = U PL , II I 0 1 n II κ R { 1 − κ r 1 + n II n II n II + 1 2 1 n II + − 1 2 π n II 1 2 n II [ Γ 1 − 1 2 n II , − κ R n II + 3 8 n II 2 Γ − 1 2 n II , − κ R n II − Γ 1 − 1 2 n II , − 1 n II − 3 8 n II 2 Γ − 1 2 n II , − 1 n II ] } , (15a)

u II r = U PL , II I 0 1 n II κ R { − 1 2 π n II 1 2 n II [ Γ 1 − 1 2 n II , − κ R n II + 3 8 n II 2 Γ − 1 2 n II , − κ R n II − Γ 1 − 1 2 n II , − κ r n II − 3 8 n II 2 Γ − 1 2 n II , − κ r n II ] } . (15b)

The velocity profile for 0 ≤ κ r < κ R − δ is given as Equation 15c:

u I r = u II R − δ + U PL , I I 0 1 n I κ R κ R − δ 1 + n I n I − κ r 1 + n I n I n I + 1 2 1 n I . (15c)

2.3 Simplified analytical solutions

Under the assumption of a thin EDL, i.e., κ R ≫ 1, which is satisfied in typical microfluidic applications (Li, 2004; Chang and Yeo, 2009; Alizadeh et al., 2021), we can view the EDL as locally flat with no overlap in cylindrical microchannels. Thus, the EDL potential in Equation 5 can be simplified as Equation 16 (Ghosal, 2004; Olivares et al., 2009):

ψ r ≈ ζ e − κ R − r . (16)

Consequently, the integral terms in Equations 7a,b, can both be simplified, leading to Equations 17a, 17b:

u I r = u II R − δ − U PL , I e − κ R − r n I − e − κ δ n I for 0 ≤ r ≤ R − δ , (17a)

u II r = − U PL , II e − κ R − r n II − 1 for R − δ ≤ r ≤ R . (17b)

The centerline electroosmotic velocity of the core–annular flow of power-law fluids, U EO = u I 0, is presented in Equation 18:where we have assumed e − κ R / n I = e − κ R / n II = 0 for κ R ≫ 1 and U PL , I and U PL , II are the generalized Smoluchowski velocities for the power-law core and annular fluids, respectively, as defined in Equation 12. We note that Equation 18 reduces to the same formula as in our prior study for U EO of power-law fluids in a parallel-plate microchannel (Leung et al., 2026), indicating the independence of the centerline electroosmotic velocity (or electroosmotic slip velocity) on the channel geometry under the thin EDL assumption. It, therefore, can also be used to predict the experimentally measured U EO of shear-thinning xanthan gum solutions in a rectangular microchannel (Bentor et al., 2024), as demonstrated in our previous study (Leung et al., 2026), for which the process is not provided here. The volumetric flow rates of the core and annular fluids can each be obtained as follows (Equations 19a, 19b):

U EO = 1 − e − κ δ n II U PL , II + e − κ δ n I U PL , I , (18)

Q I = 2 π ∫ 0 R − δ u I r r d r = π R − δ 2 1 − e − κ δ n II U PL , II + π R − δ − n I κ 2 + n I 2 κ 2 e − κ δ n I U PL , I , (19a)

Q II = 2 π ∫ R − δ R u II r r d r = 2 π U PL , II δ 2 R − δ 2 − n II κ R − n II κ − R − δ − n II κ e − κ δ n II . (19b)

One can observe that both U EO and Q I of the core fluid have contributions from U PL (see Equation 12) in both the core (denoted by the subscript I) and annular (denoted by the subscript II) power-law fluids, the extent of which depends on the relative thickness of the annular fluid to the EDL, κ δ (see Figure 1).

3 Results and discussion

We first assess the accuracy of the above-derived approximate and simplified analytical solutions for electroosmotic core–annular flow of power-law fluids by comparing the predictions with those from the exact solution. The validated analytical solutions are then used to examine the effects of the core and annular fluid parameters on the electroosmotic velocity and flow rate. The following material properties and constants were employed in the analysis: Debye length, κ − 1 = 100 nm; fluid permittivity regardless of the rheology, ε = 6.91 × 10 − 10 F / m; Newtonian fluid viscosity, η = 1 mPa · s; flow consistency index of power-law fluids, m = 50 mPa · s 1 / 2; electric field, E = 100 V / cm; and wall zeta potential, ζ = − 25 mV.

3.1 Validation of the approximate and simplified analytical solutions

To validate the approximate and simplified analytical solutions, we obtain the exact solution from Equations 7a,b, for the special case of n I = 1 / 2 for the core fluid (i.e., shear thinning) and n II = 1 for the annular fluid (i.e., Newtonian), where the electroosmotic velocity profile in the core fluid (0 ≤ r ≤ R − δ) is given in Equation 20a:while that in the annular fluid (R − δ ≤ r ≤ R) is given by Equation 9, i.e., Equation 20b:

u I r = − ε ζ E η 1 − I 0 κ R − δ I 0 κ R + κ − ε ζ E m I 2 12 I 0 2 κ R κ R − δ 3 F 3 2 3 2 , 3 2 ; 2 , 5 2 , 3 ; κ R − δ 2 − κ r 3 F 3 2 3 2 , 3 2 ; 2 , 5 2 , 3 ; κ r 2 , (20a)

u II r = u N r = − ε ζ E η 1 − I 0 κ r I 0 κ R . (20b)

Here, F 3 2 a 1 , a 2 ; b 1 , b 2 , b 3 ; z is a generalized hypergeometric function (Abramowitz and Stegun, 1965).

Figure 2A shows the comparison of the predicted electroosmotic velocity profiles from the exact solution in Equations 20a,b, with those from our approximate analytical solution in Equations 11a–c for the case of κ R − δ ≥ 1 and κ R > 1, which are valid in typical electrokinetic microchannel flows (Li, 2004; Chang and Yeo, 2009). For convenience, the predicted velocity data were all normalized by the centerline electroosmotic velocity U EO at κ R = 50. The thickness of the Newtonian annular fluid was fixed at κ δ = 0.5. We observe that the predictions of the approximate analytical solution show good agreement with those of the exact solution when κ R ≥ 5. A visible deviation appears at κ R = 2, where the approximate solution predicts a 7.6% higher U EO than the exact solution. This is mainly attributed to the piecewise approximation method used to derive Equations 11a–c, which, as noted above, has been reported to introduce the largest error when κ R approaches unity (Zhao and Yang, 2013b). As the electrostatic driving force is present only in the EDL with a typical thickness of three times the Debye length, i.e., 3 κ − 1 (Hunter, 1981; Masliyah and Bhattacharjee, 2006), the electroosmotic flow exhibits a nearly plug-like velocity profile in the bulk of channels with κ R ≥ 5. In smaller channels, the EDL begins to overlap, and the electroosmotic velocity profile transitions to a parabolic shape at κ R = 2, with a significantly lowered magnitude. Figure 2B compares the predicted electroosmotic velocity profiles from the exact solution in Equations 20a,b, with those from our approximate analytical solution in Equations 14a,b, for the case of κ R − δ < 1 and κ R < 1, which is a condition typically encountered in nanochannels (Chang and Yeo, 2009; Alizadeh et al., 2021). Here, κ δ = 0.1 was chosen to prevent δ from exceeding R. The predictions from the approximate and exact solutions are indistinguishable at both κ R = 0.3 and 0.5, while the former predicts a 6.8% smaller U EO than the exact solution at κ R = 0.8 because of a similar reason to that for the κ R = 2 case in Figure 2A.

Figure 3 shows a comparison of the electroosmotic velocity profiles (normalized by the centerline electroosmotic velocity U EO at κ R = 50) obtained from the simplified analytical solutions in Equations 17a,b, with those from the exact solutions in Equations 20a,b . The same shear-thinning core and Newtonian annular flow as shown in Figure 2 is considered here, and the thickness of the Newtonian annular fluid is again fixed at κ δ = 0.5. One can observe that at κ R = 10, the simplified analytical solution noticeably overestimates the exact electroosmotic velocity in the bulk of the core fluid by 3.4%. This deviation arises from the thin-EDL approximation in Equation 16, which has been reported not fully capture the effects of channel curvature on the EDL potential (Ghosal, 2004; Petsev and Lopez, 2006). For the case of κ R = 50, the simplified analytical solution provides a much more accurate prediction of the electroosmotic velocity profile, with a relative error from the exact solution of no more than 0.6%. Since this κ R value is commonly encountered in electrokinetic microchannel flows (Ghosal, 2004; Li, 2004; Chang and Yeo, 2009; Alizadeh et al., 2021), our simplified analytical solutions, including Equations 17a,b, for the electroosmotic velocity profile, Equation 18 for the centerline electroosmotic velocity U EO, and Equations 19a,b, for the electroosmotic flow rate under the thin EDL assumption, should be sufficient for the prediction of electroosmotic core–annular flow of power-law fluids.

3.2 Use of the simplified analytical solutions for flow characterization

3.2.1 Effect of the annular fluid thickness

Figure 4A shows the electroosmotic velocity profiles obtained from the simplified analytical solutions in Equations 17a,b, for a shear-thinning core and Newtonian annular flow when the annular fluid thickness δ varies from 0 for an entirely shear-thinning fluid to R for an entirely Newtonian fluid in a cylindrical microchannel with κ R = 50. The flow behavior index of the core fluid was fixed at n I = 0.75. The predicted velocities are normalized by the centerline electroosmotic velocity in the purely shear-thinning fluid, U EO , δ = 0. We observe a rapid increase in core fluid velocity with increasing δ because of the much smaller viscosity of the Newtonian annular fluid than that of the shear-thinning core fluid. In other words, the less-viscous annular fluid acts as a lubrication layer that reduces the effective viscosity within the EDL, thereby enhancing the core fluid flow. In contrast, the velocity profile of the annular fluid follows that for δ = R for a purely Newtonian fluid flow because u II r in Equation 17b is independent of both δ and the rheology of the core fluid. As a result, with increasing δ, the overall electroosmotic velocity profile of the core–annular flow approaches that of a purely Newtonian fluid flow and becomes nearly identical to the latter at κ δ = 5 (data not shown in Figure 4A). This phenomenon takes place because electroosmotic flow is induced by the electrostatic force confined within the EDL (Hunter, 1981) and hence becomes independent of the core fluid when κ δ ≫ 1.

Figure 4B shows the normalized electroosmotic velocity profiles of for a shear-thinning core and Newtonian annular flow with n I = 1.25 at different annular fluid thicknesses, κ δ, for κ R = 50. The flow consistency index of the core fluid m I was set to be identical to that of the shear-thinning fluid in Figure 4A. We observe from Figure 4B that increasing κ δ also rapidly increases the electroosmotic velocity of the shear-thickening core fluid because of its much higher viscosity than the Newtonian annular fluid. Moreover, the rate of increment is much higher than that for the case of shear-thinning core fluid in Figure 4A, Specifically, the normalized centerline electroosmotic velocity U EO / U EO , δ = 0 increases from 1 at δ = 0 for a purely shear-thickening fluid flow to 81.2 at δ = R for a purely Newtonian fluid flow, compared to the corresponding change from 1 to 20.5 in the for shear-thinning core and Newtonian annular flows. This is because the electroosmotic flow increases the viscosity of the shear-thickening fluid while decreasing that of the shear-thinning fluid, such that the difference in viscosity between the core and annular fluids becomes even larger in the former, leading to a stronger influence of κ δ. Meanwhile, the absolute U EO in the shear-thickening core fluid is only approximately one-quarter of that in the shear-thinning core fluid because of the discrepancy in their viscosity values.

It is important to note that the annular fluid thickness δ may be experimentally controlled differently depending on the physical realization of the core–annular configuration. When the annular fluid represents a wall-induced polymer-depletion or adsorption layer, δ can be influenced by the polymer formulation, polymer–wall interactions, wall surface properties, buffer conditions, and imposed flow conditions (de Gennes, 1981; Ma and Graham, 2005; Park et al., 2019). Its relative effects may also be regulated through the channel radius and EDL thickness, which determine δ / R and κ δ, respectively. Alternatively, a prescribed core–annular configuration may be produced by coaxial injection of separate core and sheath fluids, for which the inlet flow-rate ratio controls the core radius and hence δ (Rodriguez-Trujillo et al., 2020). Our analysis here treats δ as an independently specified or experimentally measured parameter, rather than predicting it from the underlying molecular interactions.

3.2.2 Effect of the core fluid rheology

Figure 5A shows the predicted centerline electroosmotic velocity U EO (normalized by the corresponding U EO , δ = 0) from Equation 18 for the power-law core and Newtonian annular flows at different κ δ. The flow behavior index n I of the core fluid was varied from 0.6 to 1.4, i.e., transitioning from shear-thinning to Newtonian and then shear-thickening behavior. Representative shear-thinning fluids include aqueous polymer solutions, such as xanthan gum (Leung et al., 2025) and carboxymethyl cellulose (Benchabane and Bekkour, 2008). Representative shear-thickening fluids include concentrated suspensions, such as cornstarch in water and glass spheres in oil (Brown and Jaeger, 2009). A similar range of flow behavior index values has frequently been used in previous studies for illustration purposes (Olivares et al., 2009; Huang et al., 2014; Deng et al., 2021; Deng, 2022; Deng and Xiao, 2022). In contrast, the flow consistency index of the core fluid was fixed at m I = 50 mPa · s 1 / 2 to isolate the influence of n I on U EO in theory. It is acknowledged that the flow behavior and consistency indices of non-Newtonian fluids often vary simultaneously in practice. Nevertheless, the influence of m I is retained explicitly in the analytical solutions through the generalized Smoluchowski velocity U PL in Equation 12. Specifically, at a fixed n I, U PL , I ∝ m I − 1 / n I in Equation 18. Hence, a larger m I suppresses the electroosmotic contribution from the core fluid, thereby enhancing the lubrication effect of the Newtonian annular fluid.

Consistent with the observation in Figure 4, U EO increases with increasing κ δ for each n I because the power-law core fluid has a greater viscosity than the Newtonian annular fluid or, more accurately, U PL , I < U PL , II in Equation 18. The rate of this enhancement, as reflected by the slope of U EO with respect to κ δ, apparently decreases at larger values of κ δ. It becomes essentially negligible for κ δ > 3, wherein the annular fluid is thicker than the EDL, and hence, the electric driving force acting on the core fluid approaches zero. We also observe from Figure 5A that the normalized U EO increases with increasing n I at the same κ δ because a more viscous core fluid leads to an increased viscosity difference (or more accurately, the difference between U PL , I and U PL , II in Equation 18) between the core and annular fluids. In other words, the lubrication effect of the less viscous and, in turn, more electroosmotic Newtonian annular fluid is enhanced when the power-law core fluid becomes more viscous and hence less electroosmotic with increasing n I. For example, at κ δ = 5 , U EO / U EO , δ = 0 = 97.6 for the shear-thickening core fluid with n I = 1.4, which is more than one order of magnitude larger than U EO / U EO , δ = 0 = 7.81 for the shear-thinning core fluid with n I = 0.6. It is, however, important to note that the absolute U EO becomes smaller at larger n I due to the increased viscosity of the core fluid as it becomes more shear-thickening.

Figure 5B shows the predicted volumetric flow rate of the core fluid Q I (normalized by the corresponding Q I , δ = 0) from the simplified analytical solution in Equation 19a for the power-law core and Newtonian annular flows considered in Figure 5A. We observe that increasing κ δ causes a non-monotonic change in Q I for each n I of the core fluid because it enhances U EO, as demonstrated in Figure 5A, while reducing the cross-sectional area of the core fluid passage. Moreover, since U EO barely increases for κ δ > 3, one can expect that the maximum Q I should appear at approximately κ δ = 3, as demonstrated in Figure 5B. Specifically, the optimum κ δ is 3.2 for all tested n I values because the core-fluid-associated term can be safely neglected in Equation 19a when κ δ > 3. At this κ δ, we can obtain Q I / Q I , δ = 0 = 87.4 for the shear-thickening core fluid with n I = 1.4, compared to 6.77 for the shear-thinning core fluid with n I = 0.6. Both of these normalized Q I values are approximately 10% smaller than the corresponding normalized U EO in Figure 5A because the former must account for the reduced area of the core-fluid flow, as noted above. However, the ratio of the maximum Q I between the cases with n I = 1.4 and n I = 0.6, i.e., 87.4 / 6.77 = 12.9, is approximately identical to the corresponding ratio of U EO, which is 97.6 / 7.81 = 12.5. This similar trend between Q I and U EO further justifies the enhanced lubrication effect of the more electroosmotic Newtonian annular fluid when the power-law core fluid becomes more viscous and, in turn, less electroosmotic with increasing n I. However, similar to U EO, the absolute Q I also decreases when the core fluid becomes less shear-thinning or more shear-thickening with increasing n I.

3.2.3 Effect of annular fluid rheology

As an extension of the power-law core and Newtonian annular flows shown in Figure 5A, we present in Figure 6A the predicted centerline electroosmotic velocity U EO (normalized by the corresponding U EO , δ = 0) for the more general shear-thickening core and power-law annular flows at different κ δ. The flow behavior index of the core fluid was fixed at n I = 1.2, while that of the annular fluid was varied from n II = 0.6 to 1.6. The core and annual fluids were assumed to have an equal flow consistency index, m I = m II = 50 mPa · s n. We note that n II > n I corresponds to the situation where polymers adsorb to the channel wall such that the annular fluid becomes more viscous and shear-thickening than the core fluid (de Gennes, 1981; Olivares et al., 2009). Therefore, the annular fluid acts to impede the electroosmotic flow of the less viscous and hence more electroosmotic core fluid, in contrast to the lubrication effect in the power-law core and Newtonian annular flows presented in Figure 5. Increasing κ δ, as shown in Figure 6A, lowers the value of U EO, with the rate of decrease also becoming smaller and essentially negligible for κ δ > 3 because of the same reason as noted above for Figure 5A. The case n II = n I is equivalent to a purely shear-thickening electroosmotic flow, and hence, U EO remains unchanged regardless of the value of κ δ. In contrast, 1 ≤ n II < n I can be understood as the situation of polymer depletion adjacent to the channel wall, wherein the annular fluid becomes less viscous and shear-thickening (de Gennes, 1981; Olivares et al., 2009). Consequently, increasing κ δ enhances the electroosmotic core flow, leading to U EO / U EO , δ = 0 > 1, whose rate, however, becomes smaller at larger κ δ, similar to the trend observed in Figure 5A for n II < n I = 1. Further reducing n II such that n II < 1 < n I leads to a stronger increase in the normalized U EO because of the increasing viscosity difference (or more accurately, the difference between U PL , I and U PL , II in Equation 18 as noted above) between the shear-thickening core and shear-thinning annular fluids.

Figure 6B presents the variations in the normalized flow rate of the shear-thickening core fluid, Q I / Q I , δ = 0, with κ δ for varying n II of the power-law annular fluid. One can observe that for cases with n II < n I, Q I first increases and then decreases with increasing κ δ because of the same reason as noted above for the electroosmotic power-law core/Newtonian annular flow in Figure 5B. Interestingly, the optimum κ δ to achieve the maximum Q I also exhibits a weak non-monotonic trend with respect to n II. As highlighted in Figure 6B, the optimum κ δ increases slightly from 2.01 at n II = 0.6 to 2.17 at n II = 0.8, followed by a slight decrease to 1.84 at n II = 1.0. For all the other cases with n II ≥ n I, Q I continuously decreases with increasing κ δ because U EO becomes smaller (for n II > n I, see Figure 6A or the core-fluid passage shrinks. Overall, the presence of the annular fluid enhances the electroosmotic flow of the power-law core fluid in terms of U EO and Q I if the former is more electroosmotic, with a greater generalized Smoluchowski velocity, U PL, in Equation 12. A larger difference between the two flow behavior index values or generalized Smoluchowski velocities leads to a stronger increase in the electroosmotic core–annular flow of power-law fluids.

4 Conclusion

We have developed a theoretical model for the generalized electroosmotic core–annular flow of power-law fluids in cylindrical microchannels. Approximate analytical solutions are obtained for the electroosmotic velocity with arbitrary power-law model parameters by employing an approximation to the modified Bessel function. They are verified to provide good accuracy for both κ R > 1 and κ R < 1, which are typically encountered in electrokinetic microfluidic and nanofluidic devices, respectively. We also derive simplified analytical solutions under the thin EDL assumption and validate them by benchmarking against the exact solution. Moreover, these simplified analytical solutions are used to examine the effects of the annular fluid thickness κ δ and the flow behavior indices of the core and annular fluids, n I and n II, on the electroosmotic velocity and flow rate. We demonstrate that the presence of a Newtonian annular fluid acts to increase U EO of the power-law core fluid with increasing κ δ or n I because of the lower viscosity of the former. It also enhances Q I of the core fluid with increasing n I at the same κ δ, while the optimum κ δ to maximize Q I remains unchanged with n I. For the more general electroosmotic core–annular flow, the power-law annular fluid increases U EO of the shear-thickening core fluid with increasing κ δ as long as n II < n I or, more generally, U PL , II > U PL , I. Under this same condition, however, the annular fluid does not always enhance Q I of the core fluid, and the optimum κ δ to maximize Q I shows a weak dependence on n II.

In summary, our derived analytical solutions clearly reveal the respective influences of the core and annular fluid rheology on the electroosmotic flow in terms of the generalized Smoluchowski velocity of power-law fluids, U PL in Equation 12. The annular fluid acts as either a lubricant or a resistant layer to the electroosmotic flow of the core fluid, depending on whether the former has a larger or smaller U PL than the latter. Our analysis here can help understand the effect of the wall-induced polymer depletion or adsorption on the electroosmotic velocity and flow rate of shear-thinning or shear-thickening fluids. It can also assist in the design and optimization of electrokinetic microfluidic devices containing two or more layers of power-law fluids. Moreover, the presented analytical solutions complement numerical simulations by providing explicit parameter dependencies, rapid evaluation of limiting and optimal conditions, and useful benchmarks for numerical models. However, we acknowledge that this analytical tractability is achieved at the cost of idealized assumptions, whereas numerical methods offer greater flexibility for complex geometries, nonlinear electrostatics, transient flows, and other conditions beyond the present model. In addition, we acknowledge the inability of the power-law model to reproduce Newtonian viscosity plateaus at low and high shear rates, elastic effects, stress relaxation, memory effects, and other important rheological phenomena commonly observed in complex fluids (Bird et al., 1987). In future work, we will extend the present analysis to more accurate constitutive equations for the prediction of electroosmotic flow in viscoelastic fluids.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Soumya K Srivastava, West Virginia University, United States

Reviewed by: Edtson Emilio Herrera Valencia, National Autonomous University of Mexico, Mexico

Sumit Kumar Mehta, SRM University AP, India

Jaime Ernesto Munoz Sanchez, National Polytechnic Institute (IPN), Mexico

Footnotes

c, ionic concentration; δ, thickness of the annular fluid; e, elementary charge; E, electric field; ε, fluid permittivity; κ, inverse of Debye length; k B, Boltzmann constant; m, flow consistency index; n, flow behavior index; ψ, EDL potential; Q, volumetric flow rate; R, microchannel radius; ρ e, free charge density; T, fluid temperature; u, electroosmotic velocity; U EO, centerline electroosmotic velocity; U HS, Helmholtz–Smoluchowski of Newtonian fluids; U PL, generalized Smoluchowski velocity of power-law fluids; z, valence of ions; ζ, wall zeta potential.