A Network Model for Electroosmotic and Pressure-Driven Flow in Porous Microfluidic Channels

Abstract

In this work, the network simulation method is presented as a tool for the numerical resolution of the electroosmotic and pressure-driven flow problem in microchannels with rectangular and cylindrical geometries. Based on the Brinkman equation for steady flow and constant porosity, the network model is designed using spatial discretization. An equivalent electrical circuit is obtained by establishing an analogy between the physical variable fluid velocity and electric potential. The network model is solved quickly and easily employing an electrical circuit resolution code, providing solutions for the velocity profile in the channel cross-section and the total circulating flow. After simulating two practical cases, the suitability of the grid is discussed, relating the relative errors made in the variables of interest with the number of cells used. Finally, two other applications, one for rectangular geometries and the other for cylindrical channels, show the effects the main parameters controlling the flow in these types of channels have on velocities and total flow: the zeta potential of the soil pores, applied potential and pressure gradients, and the boundary condition modified by the zeta potential in the walls of the channel.

1. Introduction

Electroosmosis is an electrokinetic phenomenon that stimulates the flow of pore fluid through a porous medium when electrodes apply a direct current (DC) electrical field [1]. Under the influence of an applied electrical field, positively charged ions are electrically attracted toward the negative electrode [2]. The motion of these ions is explained by the double-layer charge, which is the surface charge on the particle and the corresponding counter-ion charge in the pore fluid [3]. The motion of the cations drags the pore fluid within the porous medium. As the particles move, momentum is transferred to the surrounding fluid molecules, generating an electroosmotic flow between the electrodes, generally from anode to cathode [4].

The electroosmotic flow caused by an electric field is much greater than that generated by a hydraulic field. Therefore, electroosmotic dewatering offers advantages over conventional treatments due to its accelerated dewatering processes [5,6]. Electroosmotic dewatering is preferable to traditional dewatering techniques when dealing with low-permeability materials, such as fine soils, sediments, and sludge [7]. The pioneering work conducted by Casagrande to accelerate dewatering in consolidated clay soils with low hydraulic conductivity [8] is a technique that has been successfully applied to industrial applications, including soil consolidation [9,10], remediation of contaminated soils [11,12,13], food engineering [14], the dewatering of sludge [15], residues from metallurgical processes [16,17,18,19], and the nuclear industry [20], among other applications.

Dewatered flow is generated by two gradients: (i) the hydraulic gradient explained by Darcy′s law and (ii) the electrical gradient generated by the application of the DC electrical field. An electrical gradient is more effective than a hydraulic gradient for low-permeability soils. Thus, electroosmosis is effective for fine-particle soils [21].

The dewatered rate is a function of the double-layer charge, the porous medium, the pore fluid, and the electrical field [1]. Generally, electroosmotic flow is expressed through the electroosmotic permeability of the medium, which can be defined as the fluid flux per area of porous medium and the unit of electrical gradient. This coefficient depends on the zeta potential, fluid viscosity, soil porosity, and the electrical permeability of the soil [22,23]. According to the Smoluchowski theory, the zeta potential is the most important variable affecting electroosmotic flow. The zeta potential is a complex function of the interfacial chemistry between fluid and solid particles. When a solid particle moves in a fluid because of electroosmosis, a shear plane surrounding the solid particle is formed. The difference between the electrical potentials of the formed plane and the fluid is the zeta potential [24,25].

In this research, the network simulation method [26,27] is proposed as a numerical tool for solving these types of problems. The model analyzed is that of the flow induced in porous microchannels by applying both an electric field and a hydraulic potential gradient. It is defined by a type of Brinkman equation [28,29], in which steady flow and constant porosity are assumed. The network model and the process to be modeled are equivalent in that both are governed by the same equations, with a correspondence between the dependent variables in the mathematical model and the electrical variables in the equivalent circuit [30]. Once the network model design is completed, the equivalent circuit can be solved simply using an electrical circuit simulation code, such as NgSpice or PSpice [31,32,33,34]. Finally, after verifying the precision and reliability of our method with existing analytical solutions [35], two applications are presented (one for rectangular domains and the other for cylindrical ones). In these applications, the influence of the main physical parameters governing the problem on the variable fluid velocity and total circulating flow is analyzed.

The network simulation method has been successfully employed in numerous fields of applied engineering, such as ceramic coatings [36], dispersion of atmospheric pollutants [37], reinforced concrete corrosion [38], soil consolidation [39,40], seepage [41], heat transport [42,43,44], and many physical problems in engineering [30]. This technique makes use of the powerful algorithms of the circuit resolution codes [32,33] that are able to successfully cope with coupled and strong non-linear mathematical models, including Gear´s fixed time methods [45], trapezoidal integration [31], and iterative methods, such as Runge–Kutta. Thanks to these algorithms, stability in the convergence of the numerical solution is considerably improved, and a significant reduction in the local truncation error is achieved, providing high efficiency and accuracy to the network simulation method. In transient problems, time discretization is ultimately implemented automatically by the circuit resolution code [33]. In other words, the user can set a value for the time step, but the code will make the divisions necessary to reach convergence. Therefore, this feature can be considered both an advantage (the stability and convergence of the transient calculation is guaranteed) and a disadvantage (we do not have absolute control over the time step). Finally, the network simulation method requires in-depth knowledge of the problem, including the governing equations and the boundary conditions, when designing and choosing the mesh. Too-fine meshes can hinder stability and convergence in the calculations or lead to excessively long computation times, while too-coarse meshes can lead to large errors in the solution or even make convergence unattainable.

In this paper, a new numerical tool for solving problems of electroosmotic and pressure-driven flow in microporous channels is presented, based on the analogy with electrical circuits. The relative simplicity of programming this powerful and accurate tool (which requires knowledge of a few rules in electrical circuit implementation) has allowed its application to microfluidic channels with rectangular and cylindrical geometries formed by materials with different properties. This has simplified the way to identify and analyze the effects that the different variables of interest have on water flow through porous fine-particle soils. These variables are: the zeta potential of the pores, applied potential gradient, applied pressure gradient, and modified zeta potential in the channel walls.

2. Mathematical Models

The physical–mathematical model addressed here is that of electroosmotic and pressure-driven flow in porous material microchannels. The model is defined by a type of Brinkman equation [28,29], in which the flow is assumed to be stationary, and the soil has constant porosity, Equation (1). This expression, which governs the average volume flow (with units of velocity, m/s) through a porous medium, derives directly from the Navier–Stokes equations. It includes viscous forces, which allow the effects of the porosity of the soil and variations in the zeta potential near the channel walls to be considered. The last variable, which describes the intensity of the double-layer static electrical field in the limit between the soil grains and the fluid, is defined by the linearized Poisson–Boltzmann [46], Equation (7). This governs the potential distribution in the porous medium due to a zeta potential in the channel walls.

The equations presented below have been explained in detail in the work of Scales and Tait [35], but, as they form the base of our network models design, the equations are summarized here for the convenience of the reader.

2.1. Governing Equation for a Two-Parallel-Plate Channel

For this first case, which approximates the flow in channels of rectangular geometry, the Brinkman equation for volume-averaged velocity (u) due to potential ($\nabla {\mathsf{\varphi }}_z$) and pressure ($\nabla {\mathrm{P}}_{\mathrm{z}}$) gradients is expressed by:

$$\frac{{\partial }^2\mathrm{u}}{\partial {\mathrm{y}}^2}-{\mathsf{\lambda }}^2\mathrm{u}=\frac{\mathrm{n}}{\mathsf{\mu }}\left(\nabla {\mathrm{P}}_{\mathrm{z}}+\left({\mathsf{\rho }}_{\mathrm{eff}}+{\mathsf{\rho }}_{\mathrm{k}}\right)\nabla {\mathsf{\varphi }}_z\right),$$

(1)

where u (or also u_z) represents the velocity profile with which the fluid passes through the porous medium and which is a function of position y within the channel (distance perpendicular to the walls). This is: $\mathrm{u}={\mathrm{u}}_{\mathrm{z}}\left(\mathrm{y}\right)$.

In Equation (1), coefficient $\mathsf{\lambda },$ known as inverse Brinkman screening length (representing a kind of double-layer thickness measure), is defined as:

$$\mathsf{\lambda }=\sqrt{\frac{\mathrm{n}\sqrt{\mathsf{\tau }}}{\mathrm{K}}\frac{{\mathsf{\eta }}_{\mathrm{e}}}{\mathsf{\eta }}}$$

(2)

where $\mathsf{\tau }$ is the tortuosity (a pore sinuosity factor), defined as:

$$\mathsf{\tau }={\left(\frac{\mathrm{l}}{\mathrm{L}}\right)}^2$$

(3)

$\mathrm{K}$ is the intrinsic permeability of the medium, which has the expression:

$$\mathrm{K}=\frac{{\mathrm{nm}}^2}{{\mathrm{k}}_{\mathrm{o}}\sqrt{\mathsf{\tau }}}$$

(4)

where $\mathrm{m}$, the average hydraulic radius of the pores in the case of tortuous cylinders, is expressed as:

$$\mathrm{m}=\frac{{\mathrm{a}}_{\mathrm{p}}}{2}$$

(5)

Continuing with the coefficients of Equation (1), ${\mathsf{\rho }}_{\mathrm{eff}}$ is the effective charge density for the electroosmotic flow due to the zeta potential of the porous material and, in the case of cylindrical pores, takes the following expression:

$${\mathsf{\rho }}_{\mathrm{eff}}=\frac{8{\mathsf{ϵ}}_{\mathrm{w}}{\mathsf{\psi }}_{\mathrm{o}}}{{\mathrm{a}}_{\mathrm{p}}{}^2}\left(\frac{2{\mathsf{\kappa }\mathrm{a}}_{\mathrm{p}}-1}{{\left({\mathsf{\kappa }\mathrm{a}}_{\mathrm{p}}\right)}^2}-1\right)$$

(6)

Finally, ${\mathsf{\rho }}_{\mathrm{k}}$ is the charge density in the porous medium for the electroosmotic flow due to the zeta potential of the walls. Its value is obtained after solving the linearized Poisson–Boltzmann Equation (7), which governs the potential distribution in the channel, $\mathsf{\psi }$, due to the modified zeta potential of channel walls ${\mathsf{\psi }}_{\mathrm{w}}$.

$$\frac{{\partial }^2\mathsf{\psi }}{\partial {\mathrm{y}}^2}={\mathsf{\tau }\mathrm{k}}^2\mathsf{\psi }$$

(7)

With

$${\mathsf{\rho }}_{\mathrm{k}}=-{\mathsf{ϵ}}_{\mathrm{w}}{\mathsf{\kappa }}^2\mathsf{\psi }$$

(8)

where ρ_k is a function of the position y within the channel (since $\mathsf{\psi }$ is as well), ${\mathsf{\rho }}_{\mathrm{k}}\left(\mathrm{y}\right)=-{\mathsf{ϵ}}_{\mathrm{w}}{\mathsf{\kappa }}^2\mathsf{\psi }\left(\mathrm{y}\right)$. The parameter $\mathsf{\kappa }$ is the inverse Debye screening length [47], very similar to that of the inverse Brinkman screening length ($\mathsf{\lambda }$).

Boundary Conditions

In general, we will assume the non-slip condition for the channel walls, which results in zero velocities for these contours.

$${\mathrm{u}}_{\mathrm{w}1}={\mathrm{u}}_{\mathrm{w}2}=0$$

(9)

Regarding the modified zeta potential in the channel walls, the values of ${\mathsf{\psi }}_{\mathrm{w}1}$ and ${\mathsf{\psi }}_{\mathrm{w}2}$ will always take a constant value, which will be zero in the cases of uncharged walls.

Furthermore, in cases that consist of two different porous materials, continuity will be assumed for velocity at the interface that separates the two regions, that is:

$${\mathrm{u}}_{\mathrm{int}1}={\mathrm{u}}_{\mathrm{int}2}$$

(10)

In addition, the viscous shear along this interface must also be constant, with:

$${\mathsf{\eta }}_{\mathrm{e}1}\frac{\partial {\mathrm{u}}_{\mathrm{int}1}}{\partial \mathrm{y}}={\mathsf{\eta }}_{\mathrm{e}2}\frac{\partial {\mathrm{u}}_{\mathrm{int}2}}{\partial \mathrm{y}}$$

(11)

2.2. Porous Cylinder-Governing Equation and Boundary Conditions

For the case of radial geometries, the Brinkman Equation [28,29] for flow in a porous cylinder has the following expression:

$$\frac{{\partial }^2\mathrm{u}}{\partial {\mathrm{r}}^2}+\frac{1}{\mathrm{r}}\frac{\partial \mathrm{u}}{\partial \mathrm{r}}-{\mathsf{\lambda }}^2\mathrm{u}=\frac{\mathrm{n}}{\mathsf{\mu }}\left(\nabla {\mathrm{P}}_{\mathrm{z}}+\left({\mathsf{\rho }}_{\mathrm{eff}}+{\mathsf{\rho }}_{\mathrm{k}}\right)\nabla {\mathsf{\varphi }}_z\right)$$

(12)

where u (or also u_z) represents the velocity profile with which the fluid passes through the porous medium and which is a function of position r within the channel (radial distance to cylinder axis). This is: $\mathrm{u}={\mathrm{u}}_{\mathrm{z}}\left(\mathrm{r}\right)$.

The definition of the parameters and coefficients in Equation (12) is the same as for those in Equation (1). However, to obtain parameter ${\mathsf{\rho }}_{\mathrm{k}}$, Equation (8), it is necessary to solve the equation for the potential distribution in the channel ($\mathsf{\psi }$) due to the modified zeta potential of channel walls ${\mathsf{\psi }}_{\mathrm{w}}$ in radial coordinates. This is:

$$\frac{{\partial }^2\mathsf{\psi }}{\partial {\mathrm{r}}^2}+\frac{1}{\mathrm{r}}\frac{\partial \mathsf{\psi }}{\partial \mathrm{r}}={\mathsf{\tau }\mathrm{k}}^2\mathsf{\psi }$$

(13)

Regarding the boundary conditions, symmetry in the cylinder axis ($\mathrm{r}=0$) will be assumed. In the outer wall, we will consider the non-slip condition (${\mathrm{u}}_{\mathrm{wRa}}=0$) for velocity, while the modified zeta potential will take any constant value (for cases with an uncharged wall, we will have ${\mathsf{\psi }}_{\mathrm{wRa}}=0$).

3. Network Models

The network simulation method [26,27] is a numerical technique to study and simulate physical processes that can be defined using a mathematical model or a complete set of equations. Based on this model, the procedure consists of two well-defined stages: firstly, developing a network model (or electrical circuit) equivalent to the process (by spatial discretization of the governing equations) and, secondly, simulating this model using a program for solving electrical circuits, such as NgSpice or PSpice [31,32,33,34]. The network model and the physical process are formally equivalent in that both are governed by the same differential equations in finite spatial differences in terms of the elementary volume (or cell) and boundary conditions. Thanks to this equivalence and the reliability of existing circuit resolution programs (capable of obtaining their exact solution), errors in the simulation will only be associated with the choice of mesh size [48] (division of the spatial domain of the problem in n cells), as we will see in the Verification and Applications section. To achieve stability and solution convergence, circuit resolution codes, such as NgSpice [33,49], use the Newton–Raphson method to solve the non-linear equations that describe the circuit. This is an interactive algorithm that terminates once the following two conditions are met between the last iteration (k) and the current one (k + 1): (i) the currents in the non-linear branches converge to within a tolerance of 0.1% or 10⁻¹² A, and (ii) the node voltages converge to within a tolerance of 0.1% or 10⁻⁶ V.

To design the network models, we established an analogy between fluid velocity (the physical variable of the problem) and electric potential. The different electrical devices, which represent the governing equation terms (or addends), are placed between the different nodes of the elementary cell, which consists of a central node (position i) and two ends (positions i + Δ and i − Δ) in 1D domains, as can be seen in Figure 1.

It is important to note that the network simulation method would also be suitable to simply deal with the non-stationary state (if necessary) since there is a direct analogy between the addends of the governing equations with time derivatives (typical of the non-steady state) and the capacitors of an electrical circuit. Capacitors are devices that allow storage (charge or discharge) and perfectly reproduce transient phenomena. Circuit resolution codes [32,33] also incorporate different analysis modes, such as the transient mode, that is implemented using the TRAN sentence [49]. Recently, the network simulation method has been successfully used in non-stationary problems, including soil consolidation [39,40], ceramic coatings [36], and catalytic oxidation processes [50].

3.1. Discretization of the Governing Equation and Implementation of the Electric Circuit in Rectangular Domains

For rectangular domains, Equation (1), expressed in spatial finite differences according to the nomenclature in Figure 1, is as follows:

$$\left[\frac{{\mathrm{u}}_{\mathrm{i}+\Delta }-{\mathrm{u}}_{\mathrm{i}}}{\frac{{\left(\Delta \mathrm{y}\right)}^2}{2}}-\frac{{\mathrm{u}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{i}-\Delta }}{\frac{{\left(\Delta \mathrm{y}\right)}^2}{2}}\right]-{\mathsf{\lambda }}^2{\mathrm{u}}_{\mathrm{i}}=\frac{\mathrm{n}}{\mathsf{\mu }}\left(\nabla {\mathrm{P}}_{\mathrm{z}}+\left({\mathsf{\rho }}_{\mathrm{eff}}+{\mathsf{\rho }}_{\mathrm{k}}\right)\nabla {\mathsf{\varphi }}_z\right)$$

(14)

It is important take into account that in the network simulation method, $\frac{{\partial }^2\mathrm{u}}{\partial {\mathrm{y}}^2}=\frac{\partial }{\partial \mathrm{y}}\left(\frac{\partial \mathrm{u}}{\partial \mathrm{y}}\right)$ is approximated in first derivative between the ends of the cell ($\partial \mathrm{y}\cong \Delta \mathrm{y}$), while the second derivative is approximated between the central node and each one of the end nodes, with two addends with half lengths ($\partial \mathrm{y}\cong \frac{\Delta \mathrm{y}}{2}$).

According to the network simulation method, each of the addends in Equation (14) can be considered electric currents

$${\mathrm{J}}_{\mathrm{R}+\Delta }=\frac{{\mathrm{u}}_{\mathrm{i}+\Delta }-{\mathrm{u}}_{\mathrm{i}}}{\frac{{\left(\Delta \mathrm{y}\right)}^2}{2}}, {\mathrm{J}}_{\mathrm{R}-\Delta }=\frac{{\mathrm{u}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{i}-\Delta }}{\frac{{\left(\Delta \mathrm{y}\right)}^2}{2}}, {\mathrm{J}}_{\mathrm{G},\mathsf{\lambda }}={\mathsf{\lambda }}^2{\mathrm{u}}_{\mathrm{i}}, {\mathrm{J}}_{\mathrm{G},\nabla }=\frac{\mathrm{n}}{\mathsf{\mu }}\left(\nabla {\mathrm{P}}_{\mathrm{z}}+\left({\mathsf{\rho }}_{\mathrm{eff}}+{\mathsf{\rho }}_{\mathrm{k}}\right)\nabla {\mathsf{\varphi }}_z\right)$$

(15)

which are balanced at the central node of the elementary cell, so that:

$${\mathrm{J}}_{\mathrm{R}+\Delta }-{\mathrm{J}}_{\mathrm{R}-\Delta }-{\mathrm{J}}_{\mathrm{G},\mathsf{\lambda }}={\mathrm{J}}_{\mathrm{G},\nabla }$$

(16)

The terms ${\mathrm{J}}_{\mathrm{R}+\Delta }$ and ${\mathrm{J}}_{\mathrm{R}-\Delta }$ are linear and can, therefore, be implemented using individual resistors. Since $\mathrm{I}=\mathrm{V}/\mathrm{R}$, according to Ohm’s law, the values of these devices are defined as:

$${\mathrm{R}}_{\mathrm{i}+\Delta }=\frac{{\left(\Delta \mathrm{y}\right)}^2}{2}, {\mathrm{R}}_{\mathrm{i}-\Delta }=\frac{{\left(\Delta \mathrm{y}\right)}^2}{2}$$

(17)

These devices are placed between the nodes where the potential drop occurs (in the physical analogy, velocity), as reflected in Figure 1.

The terms ${\mathrm{J}}_{\mathrm{G},\mathsf{\lambda }}$ and ${\mathrm{J}}_{\mathrm{G},\nabla }$ are non-linear, so they must be implemented in the circuit as voltage-controlled current sources. In these elements, the current is specified directly through the following expressions:

$${\mathrm{G}}_{\mathsf{\lambda }}={\mathsf{\lambda }}^2{\mathrm{u}}_{\mathrm{i}}, {\mathrm{G}}_{\nabla }=\frac{\mathrm{n}}{\mathsf{\mu }}\left(\nabla {\mathrm{P}}_{\mathrm{z}}+\left({\mathsf{\rho }}_{\mathrm{eff}}+{\mathsf{\rho }}_{\mathrm{k}}\right)\nabla {\mathsf{\varphi }}_z\right)$$

(18)

where ${\mathrm{u}}_{\mathrm{i}}$ is directly read in the central node of corresponding cell i. These devices must be implemented between the central node of the cell and the common ground node.

According to Equation (8), parameter ${\mathsf{\rho }}_{\mathrm{k}}$ is a variable that depends on the potential distribution in the channel, $\mathsf{\psi }$, with which it will be necessary to solve this variable in an additional circuit (we will call it secondary), independent of the circuit that simulates the evolution of velocity $\mathrm{u}$ (we will call it main). In other words: since variable $\mathsf{\psi }$ is independent of velocity $\mathrm{u}$, the secondary circuit can be solved separately in every cell, while the main circuit has to read variable $\mathsf{\psi }$ in the secondary circuit. In this way, the expression of generator ${\mathrm{G}}_{\nabla }$ would be:

$${\mathrm{G}}_{\nabla }=\frac{\mathrm{n}}{\mathsf{\mu }}\left(\nabla {\mathrm{P}}_{\mathrm{z}}+\left({\mathsf{\rho }}_{\mathrm{eff}}-{\mathsf{ϵ}}_{\mathrm{w}}{\mathsf{\kappa }}^2{\mathsf{\psi }}_{\mathrm{i}}\right)\nabla {\mathsf{\varphi }}_z\right)$$

(19)

The value of ${\mathsf{\psi }}_{\mathrm{i}}$ will, therefore, be obtained from the reading of this variable at node i (corresponding to the same cell) of the secondary circuit, which is defined from the expression in spatial finite differences in Equations (7) and (20). Figure 2 shows the nomenclature and the network model for this variable.

$$\left[\frac{{\mathsf{\psi }}_{\mathrm{i}+\Delta }-{\mathsf{\psi }}_{\mathrm{i}}}{\frac{{\left(\Delta \mathrm{y}\right)}^2}{2}}-\frac{{\mathsf{\psi }}_{\mathrm{i}}-{\mathsf{\psi }}_{\mathrm{i}-\Delta }}{\frac{{\left(\Delta \mathrm{y}\right)}^2}{2}}\right]={\mathsf{\tau }\mathrm{k}}^2{\mathsf{\psi }}_{\mathrm{i}}$$

(20)

As before, the addends in Equation (20) can be assumed to be different electric currents that balance at the central node of the volume element, so that:

$${\mathrm{J}}_{\mathrm{R}+\Delta ,\mathsf{\psi }}=\frac{{\mathsf{\psi }}_{\mathrm{i}+\Delta }-{\mathsf{\psi }}_{\mathrm{i}}}{\frac{{\left(\Delta \mathrm{y}\right)}^2}{2}}, {\mathrm{J}}_{\mathrm{R}-\Delta ,\mathsf{\psi }}=\frac{{\mathsf{\psi }}_{\mathrm{i}}-{\mathsf{\psi }}_{\mathrm{i}-\Delta }}{\frac{{\left(\Delta \mathrm{y}\right)}^2}{2}}, {\mathrm{J}}_{\mathrm{G},\mathrm{k}}={\mathsf{\tau }\mathrm{k}}^2{\mathsf{\psi }}_{\mathrm{i}}$$

(21)

$${\mathrm{J}}_{\mathrm{R}+\Delta ,\mathsf{\psi }}-{\mathrm{J}}_{\mathrm{R}-\Delta ,\mathsf{\psi }}={\mathrm{J}}_{\mathrm{G},\mathrm{k}}$$

(22)

The terms ${\mathrm{J}}_{\mathrm{R}+\Delta \mathsf{\psi }}$ and ${\mathrm{J}}_{\mathrm{R}-\Delta ,\mathsf{\psi }}$ are linear, so they are implemented as resistors, while ${\mathrm{J}}_{\mathrm{G},\mathsf{\psi }}$ is non-linear and is made through a voltage-controlled current source.

$${\mathrm{R}}_{\mathrm{i}+\Delta ,\mathsf{\psi }}=\frac{{\left(\Delta \mathrm{y}\right)}^2}{2}, {\mathrm{R}}_{\mathrm{i}-\Delta ,\mathsf{\psi }}=\frac{{\left(\Delta \mathrm{y}\right)}^2}{2}, {\mathrm{G}}_{\mathrm{k}}={\mathsf{\tau }\mathrm{k}}^2{\mathsf{\psi }}_{\mathrm{i}}$$

(23)

where ${\mathsf{\psi }}_{\mathrm{i}}$ is directly read in the central node of corresponding cell i.

Boundary Conditions

As mentioned in Section 2.1, the values of the velocity and the modified zeta potential in the channel walls (${\mathrm{u}}_{\mathrm{w}1}$, ${\mathrm{u}}_{\mathrm{w}2}$, ${\mathsf{\psi }}_{\mathrm{w}1}$, and ${\mathsf{\psi }}_{\mathrm{w}2}$) are always constant. This means that these boundary conditions can be implemented in electrical circuits through voltage sources, which are connected between the boundary node and the common ground node. In this way, the following elements are defined for the main circuit:

$${\mathrm{V}}_{\mathrm{uw}1}={\mathrm{u}}_{\mathrm{w}1}=0 \mathrm{Between} \mathrm{the} \mathrm{lower} \mathrm{boundary} \mathrm{node} \mathrm{and} \mathrm{the} \mathrm{ground} \mathrm{node},$$

(24)

$${\mathrm{V}}_{\mathrm{uw}2}={\mathrm{u}}_{\mathrm{w}2}=0 \mathrm{Between} \mathrm{the} \mathrm{upper} \mathrm{boundary} \mathrm{node} \mathrm{and} \mathrm{the} \mathrm{ground} \mathrm{node}.$$

(25)

For the secondary circuit, the specifications of elements ${\mathrm{V}}_{\mathsf{\psi }\mathrm{w}1}$ and ${\mathrm{V}}_{\mathsf{\psi }\mathrm{w}2}$ will be, respectively, values ${\mathsf{\psi }}_{\mathrm{w}1}$ and ${\mathsf{\psi }}_{\mathrm{w}2}$.

The schematic for these boundary conditions is shown in Figure 3.

3.2. Network Model Design in Radial Domains

Given the similarities between Equations (1) and (12), the network model for radial domains is similar to that of rectangular geometry, with the addition of one more term, Figure 4. Thus, the expression in spatial finite differences of governing Equation (12) remains:

$$\left[\frac{{\mathrm{u}}_{\mathrm{i}+\Delta }-{\mathrm{u}}_{\mathrm{i}}}{\frac{{\left(\Delta \mathrm{r}\right)}^2}{2}}-\frac{{\mathrm{u}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{i}-\Delta }}{\frac{{\left(\Delta \mathrm{r}\right)}^2}{2}}\right]+\frac{{\mathrm{u}}_{\mathrm{i}+\Delta }-{\mathrm{u}}_{\mathrm{i}-\Delta }}{{\mathrm{r}}_{\mathrm{i}}\Delta \mathrm{r}}-{\mathsf{\lambda }}^2{\mathrm{u}}_{\mathrm{i}}=\frac{\mathrm{n}}{\mathsf{\mu }}\left(\nabla {\mathrm{P}}_{\mathrm{z}}+\left({\mathsf{\rho }}_{\mathrm{eff}}+{\mathsf{\rho }}_{\mathrm{k}}\right)\nabla {\mathsf{\varphi }}_z\right)$$

(26)

Each of the addends of Equation (26) can be considered an electric current, which balances with the others at the central node of the elementary cell (Equation (28)).

$${\mathrm{J}}_{\mathrm{R}+\Delta }=\frac{{\mathrm{u}}_{\mathrm{i}+\Delta }-{\mathrm{u}}_{\mathrm{i}}}{\frac{{\left(\Delta \mathrm{r}\right)}^2}{2}}, {\mathrm{J}}_{\mathrm{R}-\Delta }=\frac{{\mathrm{u}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{i}-\Delta }}{\frac{{\left(\Delta \mathrm{r}\right)}^2}{2}}, {\mathrm{J}}_{\mathrm{G},\mathrm{ur}}=\frac{{\mathrm{u}}_{\mathrm{i}+\Delta }-{\mathrm{u}}_{\mathrm{i}-\Delta }}{{\mathrm{r}}_{\mathrm{i}}\Delta \mathrm{r}}$$

$${\mathrm{J}}_{\mathrm{G},\mathsf{\lambda }}={\mathsf{\lambda }}^2{\mathrm{u}}_{\mathrm{i}}, {\mathrm{J}}_{\mathrm{G},\nabla }=\frac{\mathrm{n}}{\mathsf{\mu }}\left(\nabla {\mathrm{P}}_{\mathrm{z}}+\left({\mathsf{\rho }}_{\mathrm{eff}}+{\mathsf{\rho }}_{\mathrm{k}}\right)\nabla {\mathsf{\varphi }}_z\right)$$

(27)

$${\mathrm{J}}_{\mathrm{R}+\Delta }-{\mathrm{J}}_{\mathrm{R}-\Delta }+{\mathrm{J}}_{\mathrm{G},\mathrm{ur}}-{\mathrm{J}}_{\mathrm{G},\mathsf{\lambda }}={\mathrm{J}}_{\mathrm{G},\nabla }$$

(28)

The new term that appears with respect to the rectangular domains (${\mathrm{J}}_{\mathrm{r}}$) is non-linear, so it is implemented using a voltage-controlled current source. Thus, the elements that make up the elementary cell network model for a radial domain (Figure 4) have the following specifications:

$${\mathrm{R}}_{\mathrm{i}+\Delta }={\mathrm{R}}_{\mathrm{i}-\Delta }=\frac{{\left(\Delta \mathrm{r}\right)}^2}{2}, {\mathrm{G}}_{\mathrm{ur}}=\frac{{\mathrm{u}}_{\mathrm{i}+\Delta }-{\mathrm{u}}_{\mathrm{i}-\Delta }}{{\mathrm{r}}_{\mathrm{i}}\Delta \mathrm{r}}, {\mathrm{G}}_{\mathsf{\lambda }}={\mathsf{\lambda }}^2{\mathrm{u}}_{\mathrm{i}}, {\mathrm{G}}_{\nabla }=\frac{\mathrm{n}}{\mathsf{\mu }}\left(\nabla {\mathrm{P}}_{\mathrm{z}}+\left({\mathsf{\rho }}_{\mathrm{eff}}+{\mathsf{\rho }}_{\mathrm{k}}\right)\nabla {\mathsf{\varphi }}_z\right)$$

(29)

where ${\mathrm{u}}_{\mathrm{i}}$ is directly read in the central node of the corresponding cell i.

Again, the ${\mathsf{\rho }}_{\mathrm{k}}$ parameter depends on the potential distribution in the channel, $\mathsf{\psi }$, Equation (8), which makes it necessary to solve this variable in a secondary circuit so that the main circuit of $\mathrm{u}$ obtains the value of ${\mathsf{\rho }}_{\mathrm{k}}$ from readings of variable $\mathsf{\psi }$ in this additional circuit. Being radial coordinates, we express Equation (13) in spatial finite differences as follows:

$$\left[\frac{{\mathsf{\psi }}_{\mathrm{i}+\Delta }-{\mathsf{\psi }}_{\mathrm{i}}}{\frac{{\left(\Delta \mathrm{r}\right)}^2}{2}}-\frac{{\mathsf{\psi }}_{\mathrm{i}}-{\mathsf{\psi }}_{\mathrm{i}-\Delta }}{\frac{{\left(\Delta \mathrm{r}\right)}^2}{2}}\right]+\frac{{\mathsf{\psi }}_{\mathrm{i}+\Delta }-{\mathsf{\psi }}_{\mathrm{i}-\Delta }}{{\mathrm{r}}_{\mathrm{i}}\Delta \mathrm{r}}={\mathsf{\tau }\mathrm{k}}^2{\mathsf{\psi }}_{\mathrm{i}}$$

(30)

Figure 5 shows the nomenclature and the network model for this variable.

The electric currents that form Equation (30) are expressed as

$${\mathrm{J}}_{\mathrm{R}+\Delta ,\mathsf{\psi }}=\frac{{\mathsf{\psi }}_{\mathrm{i}+\Delta }-{\mathsf{\psi }}_{\mathrm{i}}}{\frac{{\left(\Delta \mathrm{r}\right)}^2}{2}}, {\mathrm{J}}_{\mathrm{R}-\Delta ,\mathsf{\psi }}=\frac{{\mathsf{\psi }}_{\mathrm{i}}-{\mathsf{\psi }}_{\mathrm{i}-\Delta }}{\frac{{\left(\Delta \mathrm{r}\right)}^2}{2}}, {\mathrm{J}}_{\mathrm{G},\mathsf{\psi }\mathrm{r}}=\frac{{\mathsf{\psi }}_{\mathrm{i}+\Delta }-{\mathsf{\psi }}_{\mathrm{i}-\Delta }}{{\mathrm{r}}_{\mathrm{i}}\Delta \mathrm{r}}, {\mathrm{J}}_{\mathrm{G},\mathrm{k}}={\mathsf{\tau }\mathrm{k}}^2{\mathsf{\psi }}_{\mathrm{i}}$$

(31)

and they are balanced at the central node of the elementary cell as follows:

$${\mathrm{J}}_{\mathrm{R}+\Delta ,\mathsf{\psi }}-{\mathrm{J}}_{\mathrm{R}-\Delta ,\mathsf{\psi }}+{\mathrm{J}}_{\mathrm{G},\mathsf{\psi }\mathrm{r}}={\mathrm{J}}_{\mathrm{G},\mathrm{k}}$$

(32)

The terms ${\mathrm{J}}_{\mathrm{R}+\Delta \mathsf{\psi }}$ and ${\mathrm{J}}_{\mathrm{R}-\Delta ,\mathsf{\psi }}$ are linear, so they are implemented as resistors, while ${\mathrm{J}}_{\mathrm{G},\mathsf{\psi }\mathrm{r}}$ and ${\mathrm{J}}_{\mathrm{G},\mathsf{\psi }}$ are non-linear and are made by both voltage-controlled current sources.

$${\mathrm{R}}_{\mathrm{i}+\Delta ,\mathsf{\psi }}={\mathrm{R}}_{\mathrm{i}-\Delta ,\mathsf{\psi }}=\frac{{\left(\Delta \mathrm{r}\right)}^2}{2}, {\mathrm{G}}_{\mathsf{\psi }\mathrm{r}}=\frac{{\mathsf{\psi }}_{\mathrm{i}+\Delta }-{\mathsf{\psi }}_{\mathrm{i}-\Delta }}{{\mathrm{r}}_{\mathrm{i}}\Delta \mathrm{r}}, {\mathrm{G}}_{\mathrm{k}}={\mathsf{\tau }\mathrm{k}}^2{\mathsf{\psi }}_{\mathrm{i}}$$

(33)

where ${\mathsf{\psi }}_{\mathrm{i}}$ is directly read in the central node of the corresponding cell i.

Boundary Conditions

According to what was expressed in Section 2.2, the values of the velocity and the modified zeta potential in the outer wall of the channel (${\mathrm{u}}_{\mathrm{wRa}}$, ${\mathsf{\psi }}_{\mathrm{wRa}}$) are always constant, so these boundary conditions are implemented through voltage sources, which are connected between the boundary node and the common ground node.

However, for the cylinder axis ($\mathrm{r}=0$), symmetry is assumed for variable $\mathrm{u}$ and variable $\mathsf{\psi }$. This condition can be implemented in the network simulation method by placing an infinite value resistor to prevent the flow of electrical current (adiabatic condition).

In this way, the following elements are defined for the main circuit:

$${\mathrm{V}}_{\mathrm{uwRa}}={\mathrm{u}}_{\mathrm{wRa}}=0 \mathrm{Between} \mathrm{the} \mathrm{outer} \mathrm{boundary} \mathrm{node} \mathrm{and} \mathrm{the} \mathrm{ground} \mathrm{node}$$

(34)

$${\mathrm{R}}_{\mathrm{uax}}=\infty \mathrm{Between} \mathrm{the} \mathrm{axial} \mathrm{node} \mathrm{and} \mathrm{the} \mathrm{ground} \mathrm{node}.$$

(35)

For the secondary circuit, the specifications of elements ${\mathrm{V}}_{\mathsf{\psi }\mathrm{wRa}}$ and ${\mathrm{R}}_{\mathsf{\psi }\mathrm{ax}}$ will be, respectively, the values ${\mathsf{\psi }}_{\mathrm{wRa}}$ and $\infty$.

The schematic for these boundary conditions is shown in Figure 6.

4. Verification and Applications

In this section, we will begin by verifying the solutions obtained using our numerical tool with the analytical solutions proposed by Scales and Tait [35] for rectangular geometries. Subsequently, and once the precision of the network simulation method has been tested, two applications will be presented (one for a rectangular channel and the other for a porous cylinder) for heterogeneous media composed of two materials with different properties.

The Network Simulation method has also been verified with other numerical techniques, such as the Finite Element Method [30], demonstrating that their numerical solutions have the same precision.

4.1. Verification of the Numerical Solution for the Network Simulation Method

The purpose of this verification is to quantify the deviations between our numerical solution and the analytical theoretical one and to determine how much we have to refine our grid (using parameter $n_c$, which is the number of cells into which we divide the width ($c$) of our domain) to correctly reproduce the effects of having the medium and the walls of the channel alternately charged or uncharged.

4.1.1. Charged Medium Case with Uncharged Walls

For the case of a parallel plate channel with a single charged medium (${\mathsf{\psi }}_{\mathrm{o}}\ne 0$ V) and uncharged walls (${\mathsf{\psi }}_{\mathrm{w}}=0$ V), we consider a soil that has the following characteristics:

$\mathrm{n}=0.7$, $\mathsf{\tau }=1.4$, ${\mathrm{a}}_{\mathrm{p}}=8\times {10}^{-7}$ m, ${\mathrm{k}}_{\mathrm{o}}=2$, ${\mathsf{\eta }}_{\mathrm{e}}=\mathsf{\eta }={10}^{-6}$ m²/s, $\mathsf{\mu }={10}^{-3}$ N·s/m², $\nabla {\mathrm{P}}_{\mathrm{z}}=2\times {10}^4$ N/m³, ${\mathsf{ϵ}}_{\mathrm{w}}=4.51\times {10}^{-10}$ F/m, ${\mathsf{\psi }}_{\mathrm{o}}=0.02$ V, $\mathsf{\kappa }={10}^8$ m⁻¹, ${\mathsf{\psi }}_{\mathrm{w}1}={\mathsf{\psi }}_{\mathrm{w}2}=0$ V, $\nabla {\mathsf{\varphi }}_{\mathrm{z}}=-4.92\times {10}^4$ V/m, $\mathrm{c}={10}^{-5}$ m, ${\mathrm{n}}_{\mathrm{c}}=\left[10,100\right]$.

Figure 7a,b show the results for 10, 15, 25, and 100 cells, together with the analytical solution (black line). As can be seen, the network simulation method obtains satisfactory solutions with a very low number of cells. With 100 cells, the maximum relative error made in estimating velocity $\mathrm{u}$ does not exceed 0.37% (

Table 1. Fluid velocity and flow per unit section relative errors as a function of the number of cells. Parallel plate channel. Single material. Charged medium with uncharged walls.

Number of Cells	>100	25	15	10
Max. fluid velocity ($\mathrm{u}$) error (%)	0.37	13.62	24.50	41.50
Flow per unit section ($\mathrm{U}$) error (%)	0.24	5.95	9.56	13.95

). For the quantification of the flow per unit section $\mathrm{U}$ (obtained for this geometry as the mean of velocity $\mathrm{u}$), the relative error falls to 0.24%. For this mesh, the computation time was about only 5 s on an Intel Core i7-6700 CPU 4.00 GHz computer, which demonstrates the suitability and strengths of the network simulation method for these types of problems.

4.1.2. Uncharged Medium Case with Charged Walls

For a parallel plate channel with a single uncharged medium (${\mathsf{\psi }}_{\mathrm{o}}=0$ V, $\nabla {\mathrm{P}}_{\mathrm{z}}=0$ N/m³) and charged walls (${\mathsf{\psi }}_{\mathrm{w}1}\ne 0$ V, ${\mathsf{\psi }}_{\mathrm{w}2}\ne 0$ V), we consider a soil that has the following characteristics:

$\mathrm{n}=0.7$, $\mathsf{\tau }=1.4$, ${\mathrm{a}}_{\mathrm{p}}=1.6\times {10}^{-6}$ m, ${\mathrm{k}}_{\mathrm{o}}=2$, ${\mathsf{\eta }}_{\mathrm{e}}=\mathsf{\eta }={10}^{-6}$ m²/s, $\mathsf{\mu }={10}^{-3}$ N·s/m², $\nabla {\mathrm{P}}_{\mathrm{z}}=0$ N/m³, ${\mathsf{ϵ}}_{\mathrm{w}}=4.51\times {10}^{-10}$ F/m, ${\mathsf{\psi }}_{\mathrm{o}}=0$ V, $\mathsf{\kappa }={10}^7$ m⁻¹, ${\mathsf{\psi }}_{\mathrm{w}1}=0.03$ V, ${\mathsf{\psi }}_{\mathrm{w}2}=0.01$ V, $\nabla {\mathsf{\varphi }}_{\mathrm{z}}=-4.92\times {10}^4$ V/m, $\mathrm{c}={10}^{-5}$ m, ${\mathrm{n}}_{\mathrm{c}}=\left[10,500\right]$.

Figure 8a–c show the results for 10, 15, 25, 100, and 500 cells, together with the analytical solution (black line). In this case, the network simulation method obtains satisfactory solutions with meshes greater than 25 cells. With 500 cells, the maximum relative error made in estimating velocity $\mathrm{u}$ does not exceed 0.16% (

Table 2. Fluid velocity and flow per unit section relative errors as a function of the number of cells. Parallel plate channel. Single material. Uncharged medium with charged walls.

Number of Cells	500	100	25	15	10
Max. fluid velocity (u) error (%)	0.16	2.03	31.20	64.80	90.40
Flow per unit section (U) error (%)	0.11	0.50	3.81	16.59	32.91

). For the quantification of flow per unit section $\mathrm{U}$ (obtained as the mean of velocity $\mathrm{u}$), the relative error falls to 0.11%. In this case, the computation time was approximately 15 s, longer than the simulated cases with 100 cells but still considered very short for these types of problems.

4.2. Applications

4.2.1. Two-Parallel-Plate Channel Composed of Two Porous Materials

In this subsection, eight different cases will be presented corresponding to a rectangular channel composed of two soils with different properties. The geometry, problem parameters, and characteristics of the materials that form the channel in the first case, of uncharged walls and that will serve as a reference, are the following:

${\mathrm{n}}_1=0.8$, ${\mathsf{\tau }}_1=1.4$, ${\mathrm{a}}_{\mathrm{p}1}={10}^{-6}$ m, ${\mathrm{n}}_2=0.4$, ${\mathsf{\tau }}_2=1.8$, ${\mathrm{a}}_{\mathrm{p}2}=4\times {10}^{-7}$ m, ${\mathrm{k}}_{\mathrm{o}}=2$, ${\mathsf{\eta }}_{\mathrm{e}1}={\mathsf{\eta }}_{\mathrm{e}2}=\mathsf{\eta }={10}^{-6}$ m²/s, $\mathsf{\mu }={10}^{-3}$ N·s/m², $\nabla {\mathrm{P}}_{\mathrm{z}}={10}^4$ N/m³, ${\mathsf{ϵ}}_{\mathrm{w}}=4.43\times {10}^{-10}$ F/m, ${\mathsf{\psi }}_{\mathrm{o}}=0.03$ V, ${\mathsf{\kappa }}_1={\mathsf{\kappa }}_2={10}^8$ m⁻¹, ${\mathsf{\psi }}_{\mathrm{w}1}={\mathsf{\psi }}_{\mathrm{w}2}=0$ V, $\nabla {\mathsf{\varphi }}_{\mathrm{z}}=-2\times {10}^4$ V/m, $\mathrm{c}={10}^{-5}$ m, ${\mathrm{c}}_1/\mathrm{c}=0.4$, ${\mathrm{c}}_2/\mathrm{c}=0.6$, ${\mathrm{n}}_{\mathrm{c}}=500$.

Successive cases 2 to 8 have the same geometry and material characteristics making up the channel. While certain parameters of the problem vary from the reference case (

Table 3. Parameters that are different from the reference case. Two-parallel-plate channel composed of two porous materials.

Case	Parameters
2	${\mathsf{\psi }}_{\mathrm{o}}=0.06$ V
3	$\nabla {\mathsf{\varphi }}_{\mathrm{z}}=-6\times {10}^4$ V/m
4	${\mathsf{\psi }}_{\mathrm{w}1}={\mathsf{\psi }}_{\mathrm{w}2}=0.03$ V
5	${\mathsf{\psi }}_{\mathrm{w}1}={\mathsf{\psi }}_{\mathrm{w}2}=0.05$ V
6	${\mathsf{\psi }}_{\mathrm{w}1}={\mathsf{\psi }}_{\mathrm{w}2}=0.07$ V
7	$\nabla {\mathrm{P}}_{\mathrm{z}}=3\times {10}^5$ N/m3
8	$\nabla {\mathrm{P}}_{\mathrm{z}}={10}^6$ N/m3

), the rest remain the same.

Figure 9 shows the results for the eight cases addressed, which were solved using a grid of 500 cells. As can be seen, both in case 2 (where ${\mathsf{\psi }}_{\mathrm{o}}$ has doubled) and case 3 (where $\nabla {\mathsf{\varphi }}_{\mathrm{z}}$ has tripled), the increase in velocity profiles $\mathrm{u}$ and flow per unit section $\mathrm{U}$ (

Table 4. Flow per unit section in a double-material parallel-plate channel.

Case	1 (Ref)	2	3	4	5	6	7	8
Flow per unit section U (m/s) × 10−5	8.62	17.24	25.86	9.13	9.46	9.80	9.43	11.36
% related to reference case 1	100.0	200.0	300.0	105.9	109.7	113.6	109.3	131.7

) is proportional to these changes (same dimensionless velocity profile $\mathrm{u}$/[${\mathsf{\psi }}_{\mathrm{o}}{\mathsf{ϵ}}_{\mathrm{w}}\nabla {\mathsf{\varphi }}_{\mathrm{z}}$/$\mathsf{\mu }$]). We can see how the effect of the charged walls (cases 4 to 6) considerably increases the velocity $\mathrm{u}$ of the fluid in the vicinity of these contours. However, the effect is not as considerable for the flow total $\mathrm{U}$ that circulates throughout the total section. Finally, the increase in the applied pressure gradient $\nabla {\mathrm{P}}_{\mathrm{z}}$ (cases 7 and 8) has the same qualitative effect as increasing ${\mathsf{\psi }}_{\mathrm{o}}$ or $\nabla {\mathsf{\varphi }}_{\mathrm{z}}$, although for quantitative purposes, its effect is much smaller, since to increase the total flow $\mathrm{U}$ 31.7%, it was necessary to raise $\nabla {\mathrm{P}}_{\mathrm{z}}$ 100 times.

4.2.2. Porous Cylindrical Channel Composed of Two Different Materials

As in the previous subsection, eight different cases will be presented corresponding to a porous cylinder composed of two soils with different properties. The geometry, problem parameters, and characteristics of the materials that form the channel in the first case, which will serve as a reference, are the following:

${\mathrm{n}}_1=0.5$, ${\mathsf{\tau }}_1=1.8$, ${\mathrm{a}}_{\mathrm{p}1}=3\times {10}^{-7}$ m, ${\mathrm{n}}_2=0.8$, ${\mathsf{\tau }}_2=1.5$, ${\mathrm{a}}_{\mathrm{p}2}=6·{10}^{-7}$ m, ${\mathrm{k}}_{\mathrm{o}}=2$, ${\mathsf{\eta }}_{\mathrm{e}1}={\mathsf{\eta }}_{\mathrm{e}2}=5\times {10}^{-7}$ m²/s, $\mathsf{\eta }={10}^{-6}$ m²/s, $\mathsf{\mu }={10}^{-3}$ N·s/m², $\nabla {\mathrm{P}}_{\mathrm{z}}=-2\times {10}^4$ N/m³, ${\mathsf{ϵ}}_{\mathrm{w}}=4.51\times {10}^{-10}$ F/m, ${\mathsf{\psi }}_{\mathrm{o}}=-0.02$ V, ${\mathsf{\kappa }}_1={\mathsf{\kappa }}_2={10}^8$ m⁻¹, ${\mathsf{\psi }}_{\mathrm{wRa}}=-0.02$ V, $\nabla {\mathsf{\varphi }}_{\mathrm{z}}=-3\times {10}^4$ V/m, $\mathrm{Ra}={10}^{-5}$ m, ${\mathrm{Ra}}_1/\mathrm{Ra}=0.8$, ${\mathrm{Ra}}_2/\mathrm{Ra}=0.2$, ${\mathrm{n}}_{\mathrm{c}}=500$.

Successive cases 2 to 8 have the same geometry and material characteristics making up the cylinder. While certain parameters of the problem vary from the reference case (

Table 5. Parameters that are different from the reference case. Porous cylindrical channel composed of two different materials.

Case	Parameter
2	${\mathsf{\psi }}_{\mathrm{o}}=-0.04$ V
3	$\nabla {\mathsf{\varphi }}_{\mathrm{z}}=-1\times {10}^4$ V/m
4	${\mathsf{\psi }}_{\mathrm{wRa}}=-0.04$ V
5	${\mathsf{\psi }}_{\mathrm{wRa}}=-0.08$ V
6	${\mathsf{\psi }}_{\mathrm{wRa}}=-0.12$ V
7	$\nabla {\mathrm{P}}_{\mathrm{z}}=-6\times {10}^5$ N/m3
8	$\nabla {\mathrm{P}}_{\mathrm{z}}=-2\times {10}^6$ N/m3

), the rest remain the same.

The results for the eight cases (using a 1000-cell grid) are shown in Figure 10. Analogously to rectangular domains, the variations in ${\mathsf{\psi }}_{\mathrm{o}}$ and $\nabla {\mathsf{\varphi }}_{\mathrm{z}}$ (cases 2 and 3, respectively) result in increases of equal magnitude in the velocity profiles and total circulating flow (

Table 6. Flow per unit section in a double-material porous cylindrical channel.

Case	1 (Ref)	2	3	4	5	6	7	8
Total volume flow Q (m3/s) × 10−14	5.69	11.37	1.90	5.91	6.33	6.76	6.02	6.82
% related to reference case 1	100.0	200.0	33.3	103.8	111.2	118.7	105.8	119.7

). Regarding the effect of the wall potential ${\mathsf{\psi }}_{\mathrm{wRa}}$ (cases 4 to 6), the total circulating flow is not significantly affected either, despite the fact that in the vicinity of this contour, fluid velocity does increase considerably. Finally, the increase in the applied pressure gradient $\nabla {\mathrm{P}}_{\mathrm{z}}$ (cases 7 and 8) has an impact similar to that observed in rectangular cases, both qualitatively and quantitatively.

5. Discussion and Conclusions

In this work, the network simulation method proved to be an effective tool for the numerical resolution of the physical–mathematical model of both electroosmotic and pressure-driven flow in porous microchannels. The technique, based on the analogy between the physical variables of the problem and electrical potential and intensity, allows us to solve an equivalent electrical circuit quickly (calculation times between 5 and 15 s) simply and effectively, obtaining solutions for the velocity profile and circulating flow for rectangular and radial geometries.

The method proved to be highly accurate with undemanding grids. Thus, with meshes of 100 cells, the relative errors are less than 0.4% in channels with a charged medium, although for scenarios with charged walls, the effect of the contours requires a somewhat finer crosslinking, registering relative errors of less than 0.2% with meshes of 500 cells.

In view of the applications presented, it is observed that both the increase in the zeta potential of the pores and the applied potential gradient have a similar effect, regardless of the geometry of the channel, varying, in the same proportion, the velocity profile of the fluid and, therefore, the total flow. When the applied pressure gradient is increased, the variations in velocities and circulating flow are much smaller, making it clear that the electroosmotic flow component can be much more important in porous fine-particle soils than the pressure-driven flow component.

We also found how the variation in the modified zeta potential in the channel walls has a limited effect on the velocity profiles, increasing them ostensibly in the contours but remaining without effect inside the channel. This means that the increase in total flow due to this boundary condition is much less than that achieved when we increase the zeta potential of the pores and applied potential gradient variables.