Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Interaction of Negative Bias Instability and Self-Heating Effect on Threshold Voltage and SRAM (Static Random-Access Memory) Stability of Nanosheet Field-Effect Transistors
Previous Article in Journal
IC Packaging Material Identification via a Hybrid Deep Learning Framework with CNN–Transformer Bidirectional Interaction
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Simulation of the Influence of Non-Uniform ζ Potential on Interfacial Flow

State Key Laboratory of Photon-Technology in Western China Energy, International Collaborative Center on Photoelectric Technology and Nano Functional Materials, Laboratory of Optoelectronic Technology of Shaanxi Province, Institute of Photonics & Photon Technology, Northwest University, Xi’an 710127, China
*
Author to whom correspondence should be addressed.
Micromachines 2024, 15(3), 419; https://doi.org/10.3390/mi15030419
Submission received: 9 February 2024 / Revised: 18 March 2024 / Accepted: 19 March 2024 / Published: 21 March 2024
(This article belongs to the Section A:Physics)

Abstract

:
Zeta potential (ζ potential) is a significant parameter to characterize the electric property of the electric double layer (EDL), which is important at the solid–liquid interface. Non-uniform ζ potential could be developed on a chemically uniform solid–liquid interface due to external flow. However, its influence on the flow has never been concerned. In this investigation, we numerically studied the influence of non-uniform 2D ζ potential on the flow at the solid–liquid interface. It is found, that even without any external electric field and only considering the influence of 2D ζ potential distribution, swirling flow can be generated near EDL, according to the rotational electric volume force. The streamwise vortices, which are important in the turbulent boundary layer, are theoretically predicted in this laminar flow model when considering the 2D distribution of ζ potential, implying the necessity of considering the origin of streamwise vortices of the turbulent boundary layer from the perspective of electrokinetic flow. In addition, the ζ potential distribution can promote the wall shear stress. Therefore, more attention must be paid to shear-sensitivity circumstances, like biomedical, medical devices, and in vivo. We hope that the current investigation can help us to better understand the effect of charge distribution on interfacial flow and provide theoretical guidance for the development of related applications in the future.

1. Introduction

Electric double layer (EDL) is a kind of physical structure that exists generally in the two-phase interface or even three-phase interface (including solid, liquid, and gas). Taking the solid–liquid interface as an example, when the electrolyte solution flows through the channel, charges will be redistributed on the channel wall. Under the influence of the surface charge, ions in the solution with the opposite charge of the surface charge, i.e., counter-ions, will gather on the surface of the solid, causing the corresponding positive or negative charges to be reordered. Two thin layers with different electric properties can be formed on the interface, including a Stern layer and diffuse layer, constituting the important and common interface electric structure—EDL.
Due to the widespread existence of interfacial phenomena [1,2,3,4,5], EDL and its dynamics play important roles in many fields, such as microfluidic technology [6,7,8], electrochemical analysis [9], energy systems [10,11], and biomedical applications (e.g., pharmaceutical analysis [12] and cell detection [13]). Relying on the electric properties of EDL, the electrokinetic (EK) flow (like electroosmotic flow) under the action of electric field have been broadly applied, to control flow [14], deliver drugs [15], enhance mixing in micro-/nanofluidic devices [16,17], etc. EDL also occupies a central position in solving many global problems (e.g., water desalination [18], sewage treatment [19], etc.).
The Gouy–Chapman–Stern–Grahame (GCSG) model is a widely accepted EDL model, which holds that EDL is composed of a Stern layer and a diffuse layer [20,21,22,23]. The ions in the Stern layer are considered to be adsorbed and fixed on the solid surface, while the ions in the diffuse layer but near the Stern layer cannot move freely due to the strong electrostatic force locally. As the distance from the Stern layer is beyond a certain level, the electrostatic force is not enough to restrain the thermal motion of the ions, and the ions diffuse randomly in the solution, which also shows that the diffuse layer as a dynamic structure can slide on the charged surface. The position where the ions are ready to move is called the shear plane. The ζ potential is the potential at the shear plane [23].
The ζ potential is a key parameter to characterize the electric structure in EDL. It plays an important role in interface science (e.g., interfacial reaction [24,25], self-assembly [26,27], etc.) and electrokinetic flow dynamics [28,29,30]. However, limited by electrochemical analysis techniques [31,32,33,34] on ζ potential, the hypothesis of chemically uniform interfaces with unique ζ potential was usually considered in past research. Our understanding of local ζ potential and its influence on diverse disciplines is nearly blank.
In recent years, with the in-depth study of nanoscale interface phenomena and the development of new materials with complex surfaces and superstructures, researchers have gradually begun to pay attention to the influence of non-uniform surface charge on the electrical properties of the interface and explore the factors affecting the surface charge distribution and the corresponding potential.
On one hand, relevant research shows that chemically homogeneous slab channel walls exhibit non-uniform ζ potential under the action of an outflow field. For example, Lis et al. [35] using vibrational sum frequency generation (v-SFG) spectroscopy preliminarily reveal that surface potentials vary with external flows temporally. Subsequently, Werkhoven et al. [36] showed theoretically that a pressure-driven flow can induce a strong heterogeneous surface charge and ζ potential on a chemically homogeneous channel wall. Ober et al. [37] show that liquid flow along the surface of CaF2 creates a reversible space charge gradient by v-SFG spectroscopy as well. In a current investigation, relying on a novel fluorescence photobleaching electrochemistry analyzer (FLEA), Meng et al. [38] reported a two-dimensional ζ potential induced by external flow at the solid–liquid interface.
On the other hand, non-uniform ζ potential can lead to unpredicted outcomes in diverse fields. Paratore et al. [39] demonstrated that the flow that surrounds a surface with non-uniform surface charge density (ζ potential as well) can induce complex flow patterns without a physical constraint. Yang et al. [40] revealed that gold superlattices with nanoscale variations in ζ potential distribution can significantly augment electrochemical reactions, underscoring the profound impact of ζ potential heterogeneity on electrochemical processes.
In the investigation, according to the previous experimental and theoretical reports, we hope to show how the non-uniform ζ potential distribution influences the electric volume force and its influence on the flow, through the numerical simulation method.

2. Theory

The momentum transport process of EK flow can be described by the Navier–Stokes equation with the electric volume force term, which is expressed as:
ρ f u t + u · u = P + μ 2 u + F e
where ρ f is the fluid density, P is the pressure on the fluid, μ is the dynamic viscosity. u = u 1 x 1 ^ + u 2 x 2 ^ + u 3 x 3 ^ is the velocity vector, with u 1 , u 2 and u 3 being the velocity components in x 1 , x 2 , and x 3 directions, respectively, and x 1 ^ , x 2 ^ and x 3 ^ the corresponding unit vectors in Cartesian coordinates (see Figure 1). F e is the electric volume force (EVF), which can be expressed as [41,42,43].
F e = ρ e E 1 2 E · E ε + 1 2 ρ E E ε ρ f T
where ρ e is electric charge density, E is the vector of electric field strength. In homogeneous dielectric media with incompressible approximation, the variation of ε is ignored, thus F e = ρ e E , with ρ e and E being:
ρ e = ε 2 φ
E = φ
where φ is the electric potential in the EDL, ε = ε 0 ε r is electric permittivity with ε 0 being the vacuum permittivity, and ε r is the relative permittivity of fluid. When the ζ = ζ x 1 , x 2 potential changes slowly, the potential distribution in the EDL can be approximately described by Poisson–Boltzmann theory as:
φ x 1 , x 2 , x 3 = 4 k B T e z i tanh 1 tanh e z i ζ 4 k B T e x p x 3 λ
where k B is the Boltzmann constant, T is temperature, e is the elementary charge, and z i is the relative valence state of ions in electrolytes. λ = ε k B T / 2 N A e 2 i c i z i 2 is the Debye length, which is a characteristic length of the EDL. N A is Avogadro constant, and c i is the concentration of ions in electrolytes.
Under steady state, divide both sides of Equation (1) by ρ f , and then calculate the curl on both sides with incompressible fluid relation. Thus,, we have:
u · ω ω · u = ν 2 ω + T
where ω = × u is the vorticity, ν = μ / ρ f is the kinematic viscosity coefficient, and T = 1 ρ f × F e = T 1 x 1 ^ + T 2 x 2 ^ + T 3 x 3 ^ is the driving term of vorticity equation. According to Maxwell’s equations and assuming the magnetic field is steady, by substituting Equations (3) and (4) into T , we find
T = 1 ρ f ρ e × E = ε ρ f 2 φ × φ
For convenience, let ζ i = ζ / x i , ζ i j = 2 ζ / x i x j , ζ i j k = 3 ζ / x i x j x k (with i , j , k are either 1 or 2), from Equation (5), it is obtained
φ = B N ζ 1 M x 1 ^ + B N ζ 2 M x 2 ^ + ψ A B λ M x 3 ^
2 φ = · φ = B N ζ 11 M 2 A B N ζ 1 2 ψ M + 2 A B 3 N 2 ζ 1 2 ψ M 2 + B N ζ 22 M 2 A B N ζ 2 2 ψ M + 2 A B 3 N 2 ζ 2 2 ψ M 2 + 2 ψ A 3 B 3 λ 2 M 2 ψ A B λ 2 M
and
2 φ = · ( 2 φ ) = { B N ζ 111 M + 2 B N 2 ζ 1 3 ψ 2 M + 4 A 2 B N ζ 1 3 ψ 2 M 6 A B N ζ 11 ζ 1 ψ M + 6 A B 3 N 2 ζ 11 ζ 1 ψ M 2 B 3 N 3 ζ 1 3 M 2 12 A 2 B 3 N 2 ζ 1 3 ψ 2 M 2 + 4 A B N ζ 2 ζ 12 ψ M + 4 A 2 B N ζ 1 ζ 2 2 ψ 2 M + 2 A B 3 N 2 ζ 1 ζ 22 ψ M 2 + 4 A B 3 N 2 ζ 2 ζ 12 ψ M 2 + B N ζ 1 λ 2 M 12 A 2 B 3 N 2 ζ 1 ζ 2 2 ψ 2 M 2 + 8 A 2 B 5 N 3 ζ 1 ζ 2 2 ψ 2 M 3 8 A 2 B 3 N ζ 1 λ 2 M 2 + 8 A 4 B 5 N ζ 1 λ 2 M 3 } x 1 ^ + { B N ζ 222 M + 2 B N 2 ζ 2 3 ψ 2 M + 4 A 2 B N ζ 2 3 ψ 2 M 6 A B N ζ 22 ζ 2 ψ M + 6 A B 3 N 2 ζ 22 ζ 2 ψ M 2 B 3 N 3 ζ 2 3 M 2 12 A 2 B 3 N 2 ζ 2 3 ψ 2 M 2 + 8 A 2 B 5 N 3 ζ 2 3 ψ 2 M 3 + 2 B N 2 ζ 2 ζ 1 2 ψ 2 M 2 B 3 N 3 ζ 2 ζ 1 2 ψ 2 M 2 2 A B N ζ 2 ζ 11 ψ M 4 A B N ζ 1 ζ 12 ψ M + 4 A 2 B N ζ 2 ζ 1 2 ψ 2 M + 2 A B 3 N 2 ζ 2 ζ 11 ψ M 2 + 4 A B 3 N 2 ζ 1 ζ 12 ψ M 2 + B N ζ 2 λ 2 M 12 A 2 B 3 N 2 ζ 2 ζ 1 2 ψ 2 M 2 + 8 A 2 B 5 N 3 ζ 2 ζ 1 2 ψ 2 M 3 8 A 2 B 3 N ζ 2 λ 2 M 2 + 8 A 4 B 5 N ζ 2 λ 2 M 3 } x 2 ^ + { 2 A 2 B 3 N ζ 11 λ M 2 B N ζ 11 λ M 6 A B 3 N 2 ζ 1 2 ψ λ M 2 4 A 3 B 3 N ζ 1 2 ψ λ M 2 + 8 A 3 B 5 N 2 ζ 1 2 ψ λ M 3 + 2 A B N ζ 1 2 ψ λ M + 2 A 2 B 3 N ζ 22 λ M 2 B N ζ 22 λ M 6 A B 3 N 2 ζ 2 2 ψ λ M 2 4 A 3 B 3 N ζ 2 2 ψ λ M 2 + 8 A 3 B 5 N 2 ζ 2 2 ψ λ M 3 + 2 A B N ζ 2 2 ψ λ M + ψ A B λ 3 M 8 ψ A 3 B 3 λ 3 M 2 + 8 ψ A 5 B 5 λ 3 M 3 } x 3 ^
where A x 1 , x 2 = tanh ζ x 1 , x 2 / ψ , B x 3 = exp x 3 / λ , M = A 2 B 2 1 , N = A 2 1 , and ψ = 4 k B T / e z i is thermal potential. Furthermore, after substituting Equations (8)–(10) into F e and T , we have
F e = ε 2 φ φ = ( ε B 2 M 2 N ζ 11 2 A N ζ 1 2 ψ + 2 A B 2 N 2 ζ 1 2 ψ M + N ζ 22 2 A N ζ 2 2 ψ + 2 A B 2 N 2 ζ 2 2 ψ M + 2 ψ A 3 B 2 λ 2 M ψ A λ 2 ) N ζ 1 x 1 ^ + N ζ 2 x 2 ^ + ψ A λ x 3 ^
and
T = ε ρ f 2 φ ×   φ = ε ρ f B 2 N λ M 2 { A 2 + B 2 N 2 M 2 A 2 B 2 N M ζ 2 ζ 1 2 + ζ 2 2 + [ N + 4 A 2 B 2 N M 6 A 2 ζ 2 ζ 22 + N 2 A 2 ζ 2 ζ 11 + 4 A 2 B 2 N M 1 ζ 1 ζ 12 ] + ψ A ζ 112 + ζ 222 } x 1 ^ ε ρ f B 2 N λ M 2 { 4 A ψ A 2 + B 2 N 2 M 2 A 2 B 2 N M ζ 1 ζ 2 2 + ζ 1 2 + [ N + 4 A 2 B 2 N M 6 A 2 ζ 1 ζ 11 + N 2 A 2 ζ 1 ζ 22 + 4 A 2 B 2 N M 1 ζ 2 ζ 12 ] + ψ A ζ 122 + ζ 111 } x 2 ^ + ε ρ f B 2 N 2 M 2 { 4 A ψ B 2 N M 1 ζ 12 ζ 2 2 ζ 1 2 + ζ 1 ζ 2 ζ 11 ζ 22 + ζ 2 ζ 111 + ζ 122 ζ 1 ζ 112 + ζ 222 } x 3 ^
It can be seen when ζ is non-uniform in 2D, ζ i , ζ i j , and ζ i j k can be non-zero. All the three components in Equation (12) can be nonzero, and thus, the driving term of vorticity is:
T 0
Therefore, based on Equation (6), the electric volume force F e with curl can drive the fluid to produce a swirling flow. It is worth noting that, since there is no external electric field applied, we only consider the Coulomb force in Equation (1). The influence of other electric volume forces, e.g., according to electrothermal or dielectric effect, are believed to be absent in this investigation.

3. Numerical Simulation

Based on the theory above, we can make a reasonable prediction that in the steady-state flow field, such as pressure-driven flow, there may be a non-uniform two-dimensional ζ potential distribution at the solid–liquid interface, which will induce the generation of local electrostatic charge and EVF, resulting in a swirling flow near the interface. To validate this conjecture, we established a mathematical model of non-uniform ζ potential distribution and proceeded with a numerical simulation by Comsol Multiphysics on Equations (1)–(5). The influence of non-uniform ζ potential distribution at the solid–liquid interface on the flow is investigated under steady-state approximation with the value of relative tolerance is 1 × 10 6 .

3.1. Pressure-Driven Laminar Flow

In this simulation model, the material is liquid water, with electric conductivity and pH values are σ = 100 μS/cm and p H = 9 , respectively. The thickness of the EDL on the bottom of the microchannel, i.e., Debye length, is estimated to be 11 nm. The specific physical parameters are shown in Table 1.
The basic flow is a pressure-driven steady flow with incompressible conditions. The bulk flow Reynolds number R e = U d / ν < 1 , indicating the basic flow is a typical laminar flow. Here, d = A / 2 w + h is the hydraulic diameter of the flow, U = Q / A is the bulk flow velocity, Q is flow rate, and A is the cross-sectional area of the computational regime.
In the laminar flow interface, the two sides of the rectangular model that are perpendicular to the x 1 direction are set to the inlet and outlet, respectively. At the entrance, the fully developed flow condition has been applied, with the average velocity U a v = 6.78 × 10 5 m/s. At the exit, the boundary condition is suppressing backflow (which adjusts the outlet pressure to reduce the amount of fluid entering the domain through the exit), with static pressure P 0 = 0 Pa. According to the symmetry of the geometric model and the rationality of the physical model, a no slip boundary condition has been imposed on the other four walls (including the upper wall, bottom wall, and two side walls), where the fluid velocity is zero. Electric volume force F e has been applied to the entire domain of the model (see Equation (5), similar to potential φ , and F e is dominant in the EDL at the bottom of microchannel), as shown in Figure 1.

3.2. Electrostatics Module

In this investigation, there is no external electric field applied, and the electric volume force F e is entirely generated by the 2D distribution of ζ potential. To study the effect of the electrical volume force F e due to the non-uniform ζ potential on the flow field, the electrostatic interface has been added into the simulation model. Furthermore, since the electric field is steady, the effect of the magnetic field is negligible.
The charge is conserved throughout the domain, and the electroosmotic potential φ calculated by Equation (5) was applied in the entire physical model. To further reduce the computational cost, we calculate the electric volume force by defining custom variables and analytic functions according to Equations (2)–(5).
As mentioned above, in order to investigate the effect of non-uniform ζ potential distribution on interfacial flow, a mathematical model of ζ potential distribution (as plotted in Figure 2) was established to characterize the ζ potential distribution in the solid–liquid interface at the bottom of the microchannel only (as shown in Figure 1) based on the recent report of Meng et al. [38], as
ζ x 1 , x 2 = k x 1 + A * cos 2 π k 2 x 2 + 1 ζ 0
where ζ 0 is the classical ζ potential value at the interface of water-fused silica [44]. The slope k represents how quickly the ζ potential increases in the streamwise direction (with ζ < 0 ). The dimensionless amplitude A * represents the variation of ζ potential along the spanwise direction. The wavenumber k 2 determines the possible periodicity of ζ potential in the spanwise direction.
By substituting Equation (14) into Equation (5), electroosmotic potential φ could be calculated. In turn, the electric charge density ρ e , electric field strength E , EVF, and T could be calculated accordingly by Equations (3), (4), (11) and (12).
Finally, it is worth noting that due to the choice of steady-state study method, initial values of velocity, pressure, and potential are not considered when the boundary conditions are set.

3.3. Modeling and Meshing

The computational regime is a cuboid with dimensions in the order of micrometers. In the x 1 - x 2 plane, the geometric model is built around the origin O ( 0 , 0 ) . In the investigation, the potential φ (see Equation (5)) decreases exponentially with x 3 and is dominant only in the EDL at the solid–liquid interface. To capture the influence on flow dynamics, grid cells with nanoscale resolution are required. However, to be able to capture the evolution of flow, a large computational domain is required. To compromise these two requirements, the computation regime is L = 5 μm long, w = 4 μm wide, and h = 1 μm high. Thus, we can capture the influence of micrometer variation of ζ potential on the flow. In the modeling process, the geometric model is divided into three regions in the x 3 direction, including a large bulk region ( h b = 0.9 μm) in the center of the channel and two symmetrical thin-layer regions ( h t = 50 nm) near the bottom and upper wall.
In the meshing process, mapping operation is used in the x 1 - x 2 plane, and sweeping operation is carried out upward from the bottom to the upper wall in the x 3 direction and downward from the ceiling to the under wall in the opposite of the x 3 direction, respectively, in the two thin-layer regions. The grid can be denser in the target area by the arrangement of geometric sequence. Detailed grid distribution of each region is shown in Figure 3.
During the sweeping operation in the x 3 direction, the grid in the thin-layer region is divided down to the nanometer scale to ensure resolution. The minimum and maximum grid sizes are 1.33 nm and 2.66 nm, respectively. While in the bulk region, the grid is much coarser with a fixed size of 180 nm. In the x 1 - x 2 plane, both bulk region and thin-layer region have the same minimum grid size and the same maximum grid size. In the x 1 direction, the minimum and maximum grid sizes are 6.68 nm and 20 nm, respectively. In the x 2 direction, the minimum and maximum grid sizes are 4.42 nm and 13.26 nm, respectively. Hexahedral grid elements are used for all regions in the model. Detailed parameters for complete grid statistics of the model are shown in Table 2.

3.4. Grid Independence Analysis

To ensure computational validity and stability, we conducted a grid independence analysis by systematically varying the grid size in all dimensions while keeping other parameters such as boundary conditions and solvers unchanged. The relative errors ( δ ) of the simulated velocity field between the models of different grid quantities ( N ) are plotted in Figure 4. It can be seen that δ decreases asymptotically towards zero, with the increase in grid quantity. The minimum value is 0.24%. The results show that when the grid quantity is equal to or beyond N = 7,106,000 , δ is within 1%, indicating a negligible dependency of computational results on grid density [45]. Consequently, to strike a balance between computational efficiency and precision, we adopted a grid resolution of N = 7,106,000 for subsequent analyses in this study.

4. Results

4.1. Pressure-Driven Basic Flow

In the previous investigations on interfacial flow, the influence of interfacial charge distribution and ζ potential distribution was not taken into account. The basic flow (e.g., common Poiseuille flow) is driven only by pressure.
In this case, the velocity field distribution of the basic flow is shown in Figure 5a, where x 1 * = x 1 / l , x 2 * = x 2 / l , x 3 * = x 3 / λ with characteristic length l = h / 2 = 0.5 μm and Debye length λ = 11 nm. We applied a small volume rate ( Q = 0.001   μ L / h ), where the basic flow velocity is 1.2 × 10 4 m/s at 0.5 µm from the wall, to simulate a relatively large wall strain rate (240 s−1). Slice and quiver diagrams are used to show the distribution of the intensity (where the magnitude of intensity is represented by the color) and direction of the fluid velocity u in the x 1 - x 3 planes at different streamwise positions. The direction of the fluid velocity after normalization at each position is represented by the vector arrow. Due to the absence of the external force, the velocity of the basic flow shows an approximately parabolic distribution in the spanwise direction near the wall, as can be found in Figure 5b with x 1 * = 4 .
In contrast, as shown in Figure 5c, since the EDL is extremely thin, the streamwise component of velocity, i.e., u 1 , shows a linear variation with the dimensionless vertical position x 3 * . This is an approximation of the parabolic velocity distribution of Poiseuille flow in the thin EDL.

4.2. Uniform ζ Distribution

When considering the presence of uniformly distributed ζ potential (i.e., ζ = ζ 0 ) at the solid–liquid interface, based on Equation (5), it can be seen that there is an electroosmotic potential distribution in Figure 6a. It shows a uniform distribution in the x 1 - x 2 plane and exponential decay in the vertical direction. Accordingly, F e are all in the negative x 3 direction with decreasing magnitudes from bottom to top, as shown in Figure 6b. T is exactly zero everywhere. In other words, the uniform distribution of the ζ potential induces irrotational F e . Under the influence of F e , the velocity distribution of the flow field is shown in Figure 6c. Compared with Figure 5a, it is obvious that the direction of the velocity u has not changed, which still dominated by the x 1 -directional component u 1 . Furthermore, the parabolic velocity distribution in the x 2 and x 3 directions are still present.
From the perspective of vortex dynamics, the source term T in Equation (6) does not exist in either the basic flow or the flow field considering uniform ζ potential. Therefore, it does not disturb the basic flow, nor drive the fluid to produce vorticity, and cannot cause the vorticity to deform (stretch, tilt, twist, etc.).

4.3. One-Dimensional ζ Distribution

In this subsection, we consider if an 1D distribution of ζ potential with a certain slope in the x direction is formed due to nonuniform interface charge redistribution by the basic flow. Let A * = 0 in Equation (14), we have:
ζ x 1 = k x 1 + 1 ζ 0
Taking k = 222 1/m and ζ 0 = 36 mV as an example, the 1D distribution of ζ potential see Figure 7a, and the distributions of φ and F e induced by this ζ potential are shown in Figure 7b and Figure 7c, respectively. Based on the increase in ζ in the x 1 direction, the strength of F e also shows an increasing trend from upstream to downstream. Compared with Figure 6b, although the x 3 component of F e is still dominant, the x 1 component of the F e is not zero anymore. This subsequently means that the x 2 component of T , i.e., T 2 0 , as shown in Figure 7c.
Although the flow with 1D ζ potential distribution can induce the electric volume force with F e 1 0 and F e 3 0 , in this case, we still have the streamwise and vertical components of T , i.e., T 1 = T 3 = 0 , which could be seen more intuitively in Figure 7d. According to Equation (6), ω 1 = ω 3 = 0 . There is no vorticity generated in streamwise and vertical directions through the electrical volume force, except for T 2 and ω 2 . As a result, the velocity distribution of the flow field in Figure 7e shows a clear difference from Figure 5a and Figure 6c. The dominancy of u 1 is no longer always true; in contrast, u 3 becomes dominant when 0 < x 3 * < 1 . Compared with the basic flow, both u 1 (Figure 5c) and u have exhibited significant increment. The parabolic distribution in the x 2 direction of basic flow (see Figure 5b) is no longer present when considering 1D ζ potential distribution.

4.4. Two-Dimensional ζ Distribution

In this subsection, the influence of the 2D distribution of ζ potential on flow is investigated. In Equation (14), we take k = 222 1/m, A * = 0.1 , ζ 0 = 36 mV, and k 2 = 1 / π as an example. The ζ potential exhibits a 2D distribution on a solid–fluid surface, as plotted in Figure 2.
Based on Equations (2)–(5), a 3D electroosmotic potential φ distribution is generated in the EDL at the bottom of the microchannel (as shown in Figure 8a), which in turn generates local electric charge and induces a non-uniform electric volume force. As shown in Figure 8b, F e exhibits strong x 3 -directional component. The directions of F e are mainly pointing to the bottom wall of the microchannel. The strength of F e generally shows a decaying trend along x 3 , reaching a maximum at one Debye length above the bottom wall. The strength of F e also exhibits variations along the streamwise direction and the spanwise direction. It increases gradually downstream and decreases from the center of the flow field, forming an approximate symmetry distribution along the spanwise direction.
Different from the F e induced by 0D (Figure 6b) and 1D (Figure 7c) distributions of ζ potential, in the model of 2D ζ potential, F e can be rotational in 3D. All three components of T can be non-zero, as shown in Figure 8c,d. T increases in the streamwise direction, which is consistent with that of F e . In the spanwise, T shows an approximate symmetry with respect to the center of the channel, but in opposite directions (see Figure 8d). Two maxima can be observed at x 2 * = −1.5 and 1.5.
The fluid velocity distribution under the driven of rotational F e is shown in Figure 9a. The vertical velocity component u 3 was observed near the bottom within one Debye length (i.e., 0 < x 3 * < 1 ). While u 1 component begins to dominate at x 3 * 1 . This phenomenon suggests that even though most of the electric volume force is offset by pressure gradient, such a large volume force will still have a significant impact on the velocity field. At the height of x 3 * 2.5 and higher, u is concentrated in the center of the microchannel in the spanwise to reach the maximum, which shows significant difference from those in Figure 5a, Figure 6c and Figure 7e.
In order to further explore the influence of 2D ζ potential on the interfacial flow, the velocity distribution in the x 1 - x 2 plane that one Debye length from the bottom ( x 3 * = 1 ) is shown in Figure 9b by a surface diagram. It can be seen that the velocity u of the flow, considering the effect of F e caused by 2D ζ potential distribution, is about 3 folds higher than the basic flow. It exhibits an obvious increasing trend from upstream to downstream in the streamwise direction. In this plane, u 3 , rather than u 1 , is dominant among the three components of the flow velocity, indicating that the electric volume force (after being mostly balanced by the pressure term of Equation (1)) can cause the development of flow at the bottom wall. This implies, when considering the EK flow driven by 2D ζ potential, the momentum transport near the interface can be significantly promoted.
The non-zero T further generates ω in the flow field, as shown in Figure 9c. Particularly, the streamwise component of vorticity is non-zero, i.e., ω 1 0 , and is gradually strengthened from upstream to downstream. In the upstream region, ω mainly aligns along the x 2  direction with approximately uniform distribution on the same x 3 . While in the downstream, this distribution is broken. ω becomes bending towards the streamwise direction, indicating the possible generation of streamwise vortices.
Streamwise vortices commonly exist in flow boundary layers, e.g., the well-known coherent vortex structures in the turbulent boundary layer. However, the origin of such vortices remains unclear. We hope to give an insight into this problem through the vision of electrokinetic mechanism, and provide a new approach to control the generation and evolution of vortices in boundary layers.

5. Discussion

In this section, the influence of the parameters (such as slope k , dimensionless amplitude A * , initial zeta potential ζ 0 , and wavenumber k 2 ) on the flow are systematically investigated by numerical simulations.

5.1. Influence of k

In this subsection, the influence of slope k on the flow is investigated at A * = 0.1 , ζ 0 = 36 mV, and k 2 = 1 / π . In Figure 10a, it is obvious that u 1 shows a linear increment with x 3 * , which is an approximation of the parabolic distribution caused by fluid viscosity in a thin layer. When the slope k gradually increases, a stronger variation of ζ potential along the streamwise direction is induced. The streamwise component of electric volume force gradually increases as well (see Figure 10b), accompanied by the increasing u 1 and velocity gradient u 1 / x 3 . Thus, the wall shear stress is:
τ = μ u 1 x 3
which considering 2D ζ potential distribution can be significantly stronger than that of basic flow, as can be seen more clearly in Figure 10c.

5.2. Influence of A *

The influence of A * is investigated at k = 222 1/m, ζ 0 = 36 mV and k 2 = 1 / π , as shown in Figure 5c. The comparison between the basic flow and the flow under 1D ζ potential (i.e., A * = 0 ) shows that the fluid velocity can be enhanced even if T 1 = T 3 = 0 . When spanwise variation of ζ is taken into account, both u 1 and velocity gradient u 1 / x 3 are strengthened at larger A * , indicating stronger wall shear stress as well.

5.3. Influence of ζ 0

ζ 0 is the ζ potential in quiescent fluid. It is a constant in a chemically uniform interface, determined by material properties and irrelevant to the flow conditions. Similar as k and A * , as ζ 0 is decreased from 10 mV [28], 36 mV [44], 45.8 mV [22], to 100 mV [20], u 1 (Figure 11) and its gradient u 1 / x 3 increase apparently.

5.4. Influence of k 2

At last, we discuss the influence of spanwise distributions of ζ potential on the flow, under k = 222 1/m, A * = 0.1 , and ζ 0 = 36 mV. Larger k 2 indicates a wavier distribution of ζ potential in the x 2 direction, which induces a stronger gradient of ζ potential and more peak value positions of F e as well. At k 2 = 2 / π as an example, the electric volume force F e and the corresponding fluid velocity u both exhibit wavy distribution along the spanwise direction, as shown in Figure 12a,b. Even if the peak values of F e at k 2 = 2 / π are equal to that of k 2 = 1 / π shown in Figure 8b, the peak value of u becomes decreased a little bit relative to Figure 9a. The increasing k 2 also decreases the streamwise velocity u 1 and its gradient u 1 / x 3 (Figure 12c). Therefore, a smaller τ of flow can be concluded.
All in all, it can be seen, when taking the influence of 2D ζ potential into account, the wall shear stress can be significantly enhanced. For example, at k = 222 1/m, A * = 0.1 , ζ 0 = 36 mV, and k 2 = 1 / π , the maximum τ can be up to 1.074 Pa, which is about 2.5 folds larger than the basic flow. A stronger wall shear stress indicates a stronger interaction between wall and flow. This could be important in many shear-sensitive circumstances, especially on biological materials. One example is that the flow shear stress on cells and tissues may be larger than previously expected. Studying the different control parameters of 2D potential distribution can help us better understand the effect of charge distribution on interfacial flow, and according to Equation (16), we can even realize of adjustment of wall shear stress by changing these parameters.

5.5. Influence of Other Parameters

In addition to the parameters of the zeta potential distribution discussed above, there are also several other parameters that affect the EK flow at the solid–liquid interface, for example, the electric conductivity and viscosity of the fluid. On the one hand, the electric conductivity is closely related to Debye length, which affects the distribution of electric potential φ and charge density ρ e in the EDL. This in turn affects electric volume force F e and the flow as delineated in Equations (2)–(5). On the other hand, the influence of viscosity is also important. The EK flow induced by non-uniform ζ potential is intrinsically a balance between EVF and viscosity [46,47], primarily due to the balance of electric volume force and the viscous effect. The changing of viscosity could affect both the velocity profile and wall shear stress τ near the wall. However, a comprehensive exploration on these parameters accompanied with zeta potential variation are beyond the scope of this manuscript and will be studied in future works.

6. Conclusions

In this investigation, the influence of the non-uniform 2D distribution of ζ potential on the flow at the solid–liquid interface has been investigated numerically. Since the electric volume force F e can be rotational, it further drives the fluid to produce a swirling flow. The streamwise vortices, which are important in the turbulent boundary layer, are theoretically predicted in this laminar flow model when considering the 2D distribution of ζ potential. This observation implies it is necessary to reconsider the origin of streamwise vortices in turbulent boundary layer, from the perspective of electrokinetic flow. In addition, when considering the influence of ζ potential distribution, the wall shear stress can be apparently larger relative to no ζ potential considered. This is important for shear-sensitivity circumstances, particularly in biomedical and medical devices, and in vivo. We hope the present work provides a theoretical foundation to understand the complex development and evolution mechanism through various electrical and dynamics parameters, and provides a new way to control the vortices and shear stress in boundary layers.

Author Contributions

Conceptualization, W.Z.; methodology, Y.H.; software, Y.H.; validation, W.Z.; formal analysis, Y.H. and W.Z.; investigation, Y.H.; resources, W.Z.; data curation, Y.H.; writing—original draft preparation, Y.H.; writing—review and editing, W.Z.; visualization, Y.H.; supervision, W.Z.; project administration, W.Z.; funding acquisition, W.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Natural Science Foundation of China, grant number 51927804.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Toney, M.F.; Howard, J.N.; Richer, J.; Borges, G.L.; Gordon, J.G.; Melroy, O.R.; Wiesler, D.G.; Yee, D.; Sorensen, L.B. Voltage-dependent ordering of water molecules at an electrode–electrolyte interface. Nature 1994, 368, 444–446. [Google Scholar] [CrossRef]
  2. Bazant, M.Z.; Squires, T.M. Induced-Charge Electrokinetic Phenomena: Theory and Microfluidic Applications. Phys. Rev. Lett. 2004, 92, 066101. [Google Scholar] [CrossRef]
  3. Davidson, S.M.; Andersen, M.B.; Mani, A. Chaotic Induced-Charge Electro-Osmosis. Phys. Rev. Lett. 2014, 112, 128302. [Google Scholar] [CrossRef]
  4. Fernández-Mateo, R.; Calero, V.; Morgan, H.; Ramos, A.; García-Sánchez, P. Concentration–Polarization Electroosmosis near Insulating Constrictions within Microfluidic Channels. Anal. Chem. 2021, 93, 14667–14674. [Google Scholar] [CrossRef]
  5. Sparreboom, W.; Berg, A.v.d.; Eijkel, J.C.T. Principles and applications of nanofluidic transport. Nat. Nanotechnol. 2009, 4, 713–720. [Google Scholar] [CrossRef]
  6. Whitesides, G.M. The origins and the future of microfluidics. Nature 2006, 442, 368–373. [Google Scholar] [CrossRef] [PubMed]
  7. Wu, Y.; Ma, X.; Li, K.; Yue, Y.; Zhang, Z.; Meng, Y.; Wang, S. Bipolar Electrode-based Sheath-Less Focusing and Continuous Acoustic Sorting of Particles and Cells in an Integrated Microfluidic Device. Anal. Chem. 2024, 96, 3627–3635. [Google Scholar] [CrossRef] [PubMed]
  8. Wu, Y.; Zhang, H.; Yue, Y.; Wang, S.; Meng, Y. Controllable and Scaffold-Free Formation of 3D Multicellular Architectures Using a Bipolar Electrode Array. Adv. Mater. Technol. 2024, 9, 2301626. [Google Scholar] [CrossRef]
  9. Escosura-Muñiz, A.d.l.; Ambrosi, A.; Merkoçi, A. Electrochemical analysis with nanoparticle-based biosystems. TrAC Trends Anal. Chem. 2008, 27, 568–584. [Google Scholar] [CrossRef]
  10. Karden, E.; Ploumen, S.; Fricke, B.; Miller, T.; Snyder, K. Energy storage devices for future hybrid electric vehicles. J. Power Sources 2007, 168, 2–11. [Google Scholar] [CrossRef]
  11. Hu, Y.; Yang, J.; Hu, J.; Wang, J.; Liang, S.; Hou, H.; Wu, X.; Liu, B.; Yu, W.; He, X.; et al. Synthesis of Nanostructured PbO@C Composite Derived from Spent Lead-Acid Battery for Next-Generation Lead-Carbon Battery. Adv. Funct. Mater. 2018, 28, 1705294. [Google Scholar] [CrossRef]
  12. Cui, P.; Wang, S. Application of microfluidic chip technology in pharmaceutical analysis: A review. J. Pharm. Anal. 2019, 9, 238–247. [Google Scholar] [CrossRef]
  13. Wu, Y.; Gu, Q.; Wang, Z.; Tian, Z.; Wang, Z.; Liu, W.; Han, J.; Liu, S. Electrochemiluminescence Analysis of Multiple Glycans on Single Living Cell with a Closed Bipolar Electrode Array Chip. Anal. Chem. 2024, 96, 2165–2172. [Google Scholar] [CrossRef]
  14. Tiflidis, C.; Westerbeek, E.Y.; Jorissen, K.F.A.; Olthuis, W.; Eijkel, J.C.T.; Malsche, W.D. Inducing AC-electroosmotic flow using electric field manipulation with insulators. Lab Chip 2021, 21, 3105–3111. [Google Scholar] [CrossRef]
  15. Pikal, M.J. The role of electroosmotic flow in transdermal iontophoresis. Adv. Drug Deliv. Rev. 2001, 46, 281–305. [Google Scholar] [CrossRef] [PubMed]
  16. Wu, M.; Gao, Y.; Ghaznavi, A.; Zhao, W.; Xu, J. AC electroosmosis micromixing on a lab-on-a-foil electric microfluidic device. Sens. Actuators B 2022, 359, 131611. [Google Scholar] [CrossRef]
  17. Kneller, A.R.; Haywood, D.G.; Jacobson, S.C. AC Electroosmotic Pumping in Nanofluidic Funnels. Anal. Chem. 2016, 88, 6390–6394. [Google Scholar] [CrossRef] [PubMed]
  18. Porada, S.; Weinstein, L.; Dash, R.; Wal, A.v.d.; Bryjak, M.; Gogotsi, Y.; Biesheuvel, P.M. Water Desalination Using Capacitive Deionization with Microporous Carbon Electrodes. ACS Appl. Mater. Interfaces 2012, 4, 1194–1199. [Google Scholar] [CrossRef]
  19. Wang, C.; Gao, Q. 3D Numerical Study of the Electrokinetic Motion of a Microparticle Adsorbed at a Horizontal Oil/Water Interface in an Infinite Domain. ACS Omega 2022, 7, 4062–4070. [Google Scholar] [CrossRef] [PubMed]
  20. Santiago, J.G. Electroosmotic Flows in Microchannels with Finite Inertial and Pressure Forces. Anal. Chem. 2001, 73, 2353–2365. [Google Scholar] [CrossRef]
  21. Schoch, R.B.; Lintel, H.v.; Renaud, P. Effect of the surface charge on ion transport through nanoslits. Phys. Fluids 2005, 17, 100604. [Google Scholar] [CrossRef]
  22. Brown, M.A.; Abbas, Z.; Kleibert, A.; Green, R.G.; Goel, A.; May, S.; Squires, T.M. Determination of Surface Potential and Electrical Double-Layer Structure at the Aqueous Electrolyte-Nanoparticle Interface. Phys. Rev. X 2016, 6, 011007. [Google Scholar] [CrossRef]
  23. Saha, P.; Zenyuk, I.V. Electrokinetic Streaming Current Method to Probe Polycrystalline Gold Electrode-Electrolyte Interface Under Applied Potentials. J. Electrochem. Soc. 2021, 168, 046511. [Google Scholar] [CrossRef]
  24. Putnis, C.V.; Ruiz-Agudo, E. The Mineral–Water Interface: Where Minerals React with the Environment. Elements 2013, 9, 177–182. [Google Scholar] [CrossRef]
  25. Schaefer, J.; Gonella, G.; Bonn, M.; Backus, E.H.G. Surface-specific vibrational spectroscopy of the water/silica interface: Screening and interference. Phys. Chem. Chem. Phys. 2017, 19, 16875–16880. [Google Scholar] [CrossRef]
  26. Zhang, S.; Chen, J.; Liu, J.; Pyles, H.; Baker, D.; Chen, C.-L.; Yoreo, J.J.D. Engineering Biomolecular Self-Assembly at Solid–Liquid Interfaces. Adv. Mater. 2021, 33, 1905784. [Google Scholar] [CrossRef] [PubMed]
  27. Li, X.; Chen, C.; Niu, Q.; Li, N.-W.; Yu, L.; Wang, B. Self-assembly of nanoparticles at solid–liquid interface for electrochemical capacitors. Rare Met. 2022, 41, 3591–3611. [Google Scholar] [CrossRef]
  28. Xue, J.; Zhao, W.; Nie, T.; Zhang, C.; Ma, S.; Wang, G.; Liu, S.; Li, J.; Gu, C.; Bai, J.; et al. Abnormal Rheological Phenomena in Newtonian Fluids in Electroosmotic Flows in a Nanocapillary. Langmuir 2018, 34, 15203–15210. [Google Scholar] [CrossRef] [PubMed]
  29. Zhao, W.; Liu, X.; Yang, F.; Wang, K.; Bai, J.; Qiao, R.; Wang, G. Study of Oscillating Electroosmotic Flows with High Temporal and Spatial Resolution. Anal. Chem. 2018, 90, 1652–1659. [Google Scholar] [CrossRef] [PubMed]
  30. Dutta, P.; Beskok, A. Analytical Solution of Time Periodic Electroosmotic Flows:  Analogies to Stokes’ Second Problem. Anal. Chem. 2001, 73, 5097–5102. [Google Scholar] [CrossRef] [PubMed]
  31. Huang, G.; Xu, B.; Qiu, J.; Peng, L.; Luo, K.; Liu, D.; Han, P. Symmetric electrophoretic light scattering for determination of the zeta potential of colloidal systems. Colloids Surf. A 2020, 587, 124339–124345. [Google Scholar] [CrossRef]
  32. Dukhin, A.S.; Parlia, S. Studying homogeneity and zeta potential of membranes using electroacoustics. J. Membr. Sci. 2012, 415, 587–595. [Google Scholar] [CrossRef]
  33. Bellaistre, M.C.d.; Renaud, L.; Kleimann, P.; Morin, P.; Randon, J.; Rocca, J.-L. Streaming current measurements in zirconia-coated capillaries. Electrophoresis 2004, 25, 3086–3091. [Google Scholar] [CrossRef]
  34. Gallardo-Moreno, A.M.; Vadillo-Rodríguez, V.; Perera-Núñez, J.; Bruqueab, J.M.; González-Martín, M.L. The zeta potential of extended dielectrics and conductors in terms of streaming potential and streaming current measurements. Phys. Chem. Chem. Phys. 2012, 14, 9758–9767. [Google Scholar] [CrossRef] [PubMed]
  35. Lis, D.; Backus, E.H.G.; Hunger, J.; Parekh, S.H.; Bonn, M. Liquid flow along a solid surface reversibly alters interfacial chemistry. Science 2014, 344, 1138–1142. [Google Scholar] [CrossRef]
  36. Werkhoven, B.L.; Everts, J.C.; Samin, S.; Roij, R.v. Flow-Induced Surface Charge Heterogeneity in Electrokinetics due to Stern-Layer Conductance Coupled to Reaction Kinetics. Phys. Rev. Lett. 2018, 120, 264502. [Google Scholar] [CrossRef]
  37. Ober, P.; Boon, W.Q.; Dijkstra, M.; Backus, E.H.G.; Roij, R.v.; Bonn, M. Liquid flow reversibly creates a macroscopic surface charge gradient. Nat. Commun. 2021, 12, 4102. [Google Scholar] [CrossRef] [PubMed]
  38. Meng, S.; Han, Y.; Zhao, W.; Zhu, Y.; Zhang, C.; Feng, X.; Zhang, C.; Zang, D.; Jing, G.; Wang, K. Electrokinetic origin of swirling flow on nanoscale interface. arXiv 2024, arXiv:2402.04279. [Google Scholar]
  39. Paratore, F.; Boyko, E.; Kaigala, G.V.; Bercovici, M. Electroosmotic Flow Dipole: Experimental Observation and Flow Field Patterning. Phys. Rev. Lett. 2019, 122, 224502. [Google Scholar] [CrossRef]
  40. Yang, X.; Rong, C.; Zhang, L.; Ye, Z.; Wei, Z.; Huang, C.; Zhang, Q.; Yuan, Q.; Zhai, Y.; Xuan, F.-Z.; et al. Mechanistic insights into C-C coupling in electrochemical CO reduction using gold superlattices. Nat. Commun. 2024, 15, 720. [Google Scholar] [CrossRef]
  41. Ramos, A.; Morgan, H.; Green, N.G.; Castellanos, A. Ac electrokinetics: A review of forces in microelectrode structures. J. Phys. D Appl. Phys. 1998, 31, 2338–2353. [Google Scholar] [CrossRef]
  42. Zhao, W.; Wang, G. Scaling of velocity and scalar structure functions in ac electrokinetic turbulence. Phys. Rev. E 2017, 95, 023111. [Google Scholar] [CrossRef]
  43. Nan, K.; Shi, Y.; Zhao, T.; Tang, X.; Zhu, Y.; Wang, K.; Bai, J.; Zhao, W. Mixing and Flow Transition in an Optimized Electrokinetic Turbulent Micromixer. Anal. Chem. 2022, 94, 12231–12239. [Google Scholar] [CrossRef]
  44. Hu, Z.; Zhao, W.; Chen, Y.; Han, Y.; Zhang, C.; Feng, X.; Jing, G.; Wang, K.; Bai, J.; Wang, G.; et al. Onset of Nonlinear Electroosmotic Flow under an AC Electric Field. Anal. Chem. 2022, 94, 17913–17921. [Google Scholar] [CrossRef] [PubMed]
  45. Lee, M.; Park, G.; Park, C.; Kim, C. Improvement of Grid Independence Test for Computational Fluid Dynamics Model of Building Based on Grid Resolution. Adv. Civ. Eng. 2020, 2020, 8827936. [Google Scholar] [CrossRef]
  46. Suresh, V.; Homsy, G.M. Stability of time-modulated electroosmotic flow. Phys. Fluids 2004, 16, 2349–2356. [Google Scholar] [CrossRef]
  47. Hu, Z.; Zhao, T.; Zhao, W.; Yang, F.; Wang, H.; Wang, K.; Bai, J.; Wang, G. Transition from periodic to chaotic AC electroosmotic flows near electric double layer. AIChE J. 2021, 67, e17148. [Google Scholar] [CrossRef]
Figure 1. The computational domain, coordinate system, and the boundary conditions of the physical model in the numerical simulation.
Figure 1. The computational domain, coordinate system, and the boundary conditions of the physical model in the numerical simulation.
Micromachines 15 00419 g001
Figure 2. The non-uniform 2D ζ potential distribution with k = 222 1/m, A * = 0.1 , ζ 0 = 36 mV and k 2 = 1 / π .
Figure 2. The non-uniform 2D ζ potential distribution with k = 222 1/m, A * = 0.1 , ζ 0 = 36 mV and k 2 = 1 / π .
Micromachines 15 00419 g002
Figure 3. The meshing of solution domain: (a) the thin-layer region (the schematic is scaled up to 10 times in the x 3 direction) and (b) the large bulk region.
Figure 3. The meshing of solution domain: (a) the thin-layer region (the schematic is scaled up to 10 times in the x 3 direction) and (b) the large bulk region.
Micromachines 15 00419 g003
Figure 4. Grid independence analysis. Plot of relative errors δ on the computation of velocity, versus the grid quantity N , with k = 55.5 1/m, A * = 0.1 , ζ 0 = 36 mV and k 2 = 1 / π .
Figure 4. Grid independence analysis. Plot of relative errors δ on the computation of velocity, versus the grid quantity N , with k = 55.5 1/m, A * = 0.1 , ζ 0 = 36 mV and k 2 = 1 / π .
Micromachines 15 00419 g004
Figure 5. Velocity distributions. (a) 3D distribution of fluid velocity u of basic flow. The magnitudes of u , i.e., u , is demonstrated by color. The direction of u , evaluated by u / u , is demonstrated by the arrows with unit length. (b) Plot of the streamwise velocity component of basic flow, u 1 , versus x 2 * at different x 3 * with x 1 * = 4 . (c) Plot of streamwise velocity component, u 1 , versus x 3 * at different A * , where x 1 * = 4 and x 2 * = 0 with k = 222 1/m, ζ 0 = 36 mV and k 2 = 1 / π .
Figure 5. Velocity distributions. (a) 3D distribution of fluid velocity u of basic flow. The magnitudes of u , i.e., u , is demonstrated by color. The direction of u , evaluated by u / u , is demonstrated by the arrows with unit length. (b) Plot of the streamwise velocity component of basic flow, u 1 , versus x 2 * at different x 3 * with x 1 * = 4 . (c) Plot of streamwise velocity component, u 1 , versus x 3 * at different A * , where x 1 * = 4 and x 2 * = 0 with k = 222 1/m, ζ 0 = 36 mV and k 2 = 1 / π .
Micromachines 15 00419 g005
Figure 6. Potential, electric volume force, and velocity fields generated by uniform ζ potential with ζ = 36 mV. 3D distribution of (a) electroosmotic potential φ , (b) electric volume force F e , and (c) fluid velocity u .
Figure 6. Potential, electric volume force, and velocity fields generated by uniform ζ potential with ζ = 36 mV. 3D distribution of (a) electroosmotic potential φ , (b) electric volume force F e , and (c) fluid velocity u .
Micromachines 15 00419 g006
Figure 7. Electric field and flow field generated by 1D ζ potential with k = 222 1/m and ζ 0 = 36 mV. (a) 1D distribution of ζ potential. 3D distribution of (b) electroosmotic potential φ , (c) electric volume force F e , (d) the curl of electric volume force T , and (e) fluid velocity u under this 1D ζ potential distribution.
Figure 7. Electric field and flow field generated by 1D ζ potential with k = 222 1/m and ζ 0 = 36 mV. (a) 1D distribution of ζ potential. 3D distribution of (b) electroosmotic potential φ , (c) electric volume force F e , (d) the curl of electric volume force T , and (e) fluid velocity u under this 1D ζ potential distribution.
Micromachines 15 00419 g007aMicromachines 15 00419 g007b
Figure 8. Potential, F e and T generated by 2D ζ potential with k = 222 1/m, A * = 0.1 , ζ 0 = 36 mV and k 2 = 1 / π . (a) 3D distribution of φ . (b) 3D distribution of electric volume force F e . (c,d) 3D distribution of T .
Figure 8. Potential, F e and T generated by 2D ζ potential with k = 222 1/m, A * = 0.1 , ζ 0 = 36 mV and k 2 = 1 / π . (a) 3D distribution of φ . (b) 3D distribution of electric volume force F e . (c,d) 3D distribution of T .
Micromachines 15 00419 g008
Figure 9. Velocity and vorticity fields generated by 2D ζ potential with k = 222 1/m, A * = 0.1 , ζ 0 = 36 mV and k 2 = 1 / π . (a) 3D distribution of fluid velocity u . (b) 2D velocity distribution of flow considering and without considering ζ potential distribution at x 3 * = 1 . (c) 3D distribution of the curl of fluid velocity ω , denoted by ω * , with ω * = ω / ω r e f . Here, ω r e f is the maximum ω of basic flow.
Figure 9. Velocity and vorticity fields generated by 2D ζ potential with k = 222 1/m, A * = 0.1 , ζ 0 = 36 mV and k 2 = 1 / π . (a) 3D distribution of fluid velocity u . (b) 2D velocity distribution of flow considering and without considering ζ potential distribution at x 3 * = 1 . (c) 3D distribution of the curl of fluid velocity ω , denoted by ω * , with ω * = ω / ω r e f . Here, ω r e f is the maximum ω of basic flow.
Micromachines 15 00419 g009
Figure 10. Influence of k on the flow, where x 1 * = 4 and x 2 * = 0 , A * = 0.1 , ζ 0 = 36 mV, and k 2 = 1 / π . (a) Plot of streamwise velocity component, u 1 , versus dimensional height x 3 * at different k . (b) Plot of streamwise electric volume force component F e 1 versus slope k at the height of x 3 * = 1 . (c) Plot of wall shear stress, τ , versus dimensional height x 3 * at different k .
Figure 10. Influence of k on the flow, where x 1 * = 4 and x 2 * = 0 , A * = 0.1 , ζ 0 = 36 mV, and k 2 = 1 / π . (a) Plot of streamwise velocity component, u 1 , versus dimensional height x 3 * at different k . (b) Plot of streamwise electric volume force component F e 1 versus slope k at the height of x 3 * = 1 . (c) Plot of wall shear stress, τ , versus dimensional height x 3 * at different k .
Micromachines 15 00419 g010aMicromachines 15 00419 g010b
Figure 11. Influence of ζ 0 on the flow, where x 1 * = 4 and x 2 * = 0 , k = 222 1/m, A * = 0.1 , and k 2 = 1 / π . Plot of streamwise velocity component, u 1 , versus dimensional height x 3 * at different ζ 0 .
Figure 11. Influence of ζ 0 on the flow, where x 1 * = 4 and x 2 * = 0 , k = 222 1/m, A * = 0.1 , and k 2 = 1 / π . Plot of streamwise velocity component, u 1 , versus dimensional height x 3 * at different ζ 0 .
Micromachines 15 00419 g011
Figure 12. Influence of k 2 on the flow, with k = 222 1/m, A * = 0.1 , and ζ 0 = 36 mV. 3D distribution of (a) electric volume force F e  and (b) fluid velocity u with wavenumber k 2 = 2 / π . (c) Plot of streamwise velocity component, u 1 , versus dimensional height x 3 * at different k 2 , where x 1 * = 4 and x 2 * = 0 .
Figure 12. Influence of k 2 on the flow, with k = 222 1/m, A * = 0.1 , and ζ 0 = 36 mV. 3D distribution of (a) electric volume force F e  and (b) fluid velocity u with wavenumber k 2 = 2 / π . (c) Plot of streamwise velocity component, u 1 , versus dimensional height x 3 * at different k 2 , where x 1 * = 4 and x 2 * = 0 .
Micromachines 15 00419 g012
Table 1. Properties of the material.
Table 1. Properties of the material.
Material ρ f
( k g / m 3 )
μ
( × 10 5   P a s )
ε r
Liquid water99789.5780
Table 2. Detailed information of the grid in the model.
Table 2. Detailed information of the grid in the model.
ParameterValue
Number of grid elements7,106,000
Number of grid vertices7,275,576
Minimum grid element quality1.0
Average grid element quality1.0
Grid element volume ratio 3.427 × 10 5
Grid volume 20 μm3
Minimum   grid   size   in   x 1 direction6.68 nm
Maximum   grid   size   in   x 1 direction20 nm
Minimum   grid   size   in   x 2 direction4.42 nm
Maximum   grid   size   in   x 2 direction13.26 nm
Minimum   grid   size   in   x 3 direction1.33 nm
Maximum   grid   size   in   x 3 direction180 nm
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Han, Y.; Zhao, W. Numerical Simulation of the Influence of Non-Uniform ζ Potential on Interfacial Flow. Micromachines 2024, 15, 419. https://doi.org/10.3390/mi15030419

AMA Style

Han Y, Zhao W. Numerical Simulation of the Influence of Non-Uniform ζ Potential on Interfacial Flow. Micromachines. 2024; 15(3):419. https://doi.org/10.3390/mi15030419

Chicago/Turabian Style

Han, Yu, and Wei Zhao. 2024. "Numerical Simulation of the Influence of Non-Uniform ζ Potential on Interfacial Flow" Micromachines 15, no. 3: 419. https://doi.org/10.3390/mi15030419

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