#### Abstract

This paper is focused on deriving an explicit analytical solution for the prediction of the electrostatic potential, commonly used on electrokinetic research and its related applications. Different from all other analytic techniques, this approach provides a simple way to ensure the convergence of series of solution so that one can always get accurate enough approximations. This new approach can be a useful tool in electrical field applications such as the separation of a mixture of macromolecules and the removal of contaminants in soil cleaning processes.

#### 1. Introduction

In every electrokinetic related application, the electrostatic potential is a key function to be identified in order to understand its effect on flow regimes. The methods to calculate this electrostatic potential function depend largely on approximations based on the range of values of the applied electrical field; therefore, it would be extremely useful to identify a procedure that yields a more general solution valid for a wide range of magnitude of the electrostatic potential. Thus, the same solution will be explicit for the electrostatic potential and valid for the small and large values without restrictions. To our best knowledge, such type of solution has not yet been reported in literature [1–5].

It becomes very much desired that an explicit mathematical expression for the electrostatic potential should be available. This is not necessarily easy to accomplish since the conservation of charge yields a nonlinear differential equation. This is based on the fact that the equation that is the common starting point for the description of the electrostatic potential is the Poisson-Boltzmann equation [6]: where is the dimensionless electrical potential and is the dimensionless inverse Debye length. The nonlinear term on the right-hand side of (1.1) is related to the free charge density.

The main focus of the previous studies was concentrated on how to approximate the hyperbolic function, on the right-hand side of (1.1), in order to obtain an effective analytical solution of the Poisson-Boltzmann differential equation. A very common simplification invokes the Debye-Hückel approximation is usually written as

As a result, the nonlinear Poisson-Boltzmann equation reduces to the following linear equation: Such approximation is valid for the case when [7]. Other contributions that extend the range of valid solution to propose splitting electrical potential values in three regions [8] as a way of simplifying the hyperbolic sine function: Based on the Debye-Hückel approximation, a new solution strategy for the differential equations of the electrostatic potential was proposed by Oyanader and Arce [4, 9]. The authors introduced the correction function to the inverse Debye length, . The function improves the Debye-Hückel approximation. It is a recursive function of the electrical potential and has a polynomial form whose coefficients are adjusted to predict the correct values of the hyperbolic sine function.

In this paper, we will suggest an alternative approach to the search for an explicit analytical solution for (1.1) by homotopy analysis method (HAM) [10–21]. In particular, a series solution is obtained and the constant coefficients are given by recursive formulas.

The homotopy analysis method is based on homotopy, a fundamental concept in topology and differential geometry. Briefly speaking, by means of the HAM, one constructs a continuous mapping of an initial guess approximation to the exact solution of considered equations. An auxiliary linear operator is chosen to construct such kind of continuous mapping, and an auxiliary parameter is used to ensure the convergence of solution series. The method enjoys great freedom in choosing initial approximations and auxiliary linear operators. By means of this kind of freedom, a complicated nonlinear problem can be transferred into an infinite number of simpler, linear subproblems, as shown by Liao and Tan [19].

#### 2. Approach Based on the HAM

Consider the nonlinear Poisson-Boltzmann equation: According to (2.1a) and (2.1b), it is obvious that the solution can be expressed by the base functions
as
where is a coefficient. It provides the *Solution Expression* of .

Based on (2.1a), we define a nonlinear operator where denotes the embedding parameter. Let denote an initial guess of the exact solution which satisfies the boundary conditions (2.1b), an auxiliary parameter, and an auxiliary linear operator. All of , , and will be chosen later with great freedom. Then, we construct a one-parameter family of differential equations subject to the boundary conditions Obviously, when , because of the property of any linear operator , (2.5) and (2.6) have the solution and when , since , (2.5) and (2.6) are equivalent to the original ones, (2.1a) and (2.1b), provided that Thus, based on (2.7) and (2.8), as the embedding parameter increases from to , varies continuously from the initial approximation to the exact solution . This kind of deformation is totally determined by the so-called zeroth-order deformation equations (2.5) and (2.6).

By Taylor's theorem, can be expanded in a power series of as follows: where is called the th-order homotopy-derivative of .

Fortunately, the homotopy-series (2.9) contains an auxiliary parameter . Besides, we have great freedom to choose the auxiliary linear operator , as illustrated by Liao [16]. If the auxiliary linear parameter and the nonzero auxiliary parameter are properly chosen so that the power series (2.9) of converges at then, under these assumptions, we have the so-called homotopy-series solution:

According to the fundamental theorems in calculus, each coefficient of the Taylor series of a function is unique. Thus, is unique and is determined by . Therefore, the governing equations and boundary conditions of can be deduced from the zeroth-order deformation equations (2.5) and (2.6). For brevity, define the vectors

Differentiating the zero-order deformation equation (2.5) times with respect to and then dividing by and finally setting , we have the so-called high-order deformation equation: where Using the definitions (2.4) and (2.14), we have the explicit expression Now we have to find It is clear from (2.10) that , and by using the definition of the operator , it holds Besides, one has Thus, Leibnitz's rule for derivatives of product implies that Setting in the above expression, one has It is clear that and it holds for any integer that Similarly, In this line we have that with the following recurrence formulas: From the above recurrence formulas, we have and so on. So, by means of symbolic computation software such as Mathematica and Matlab, it is not difficult to get for large values of .

Note that the high-order deformation equations (2.13) are linear ODEs. So, according to (2.11), the original nonlinear problem is transferred into an infinite number of linear ODEs.

Both the auxiliary linear operator and the initial guess are chosen under the so-called *Rule of Solution Expression*: the auxiliary linear operator and the initial guess must be chosen so that the solutions of the high-order deformation equations (2.13) exist and besides they obey the Solution Expression (2.3). So, for the solutions to obey the Solution Expression (2.3) and the boundary conditions (2.1b), we choose the initial guess of the solution:
Because the original equation (2.1a) is of second order, we simply choose such a second-order auxiliary linear operator:
which has the property
for any real coefficients and . Let denote a particular solution of (2.13), then its general solution is expressed by
where and are determined by the conditions in (2.13), that is,
In this way, we get one by one in the order

At the th-order approximation, the solution can be expressed as follows: where is a coefficient.

#### 3. Illustrative Example

Consider the system consisting of two parallel walls, each subjected to a dimensionless electrostatic potential equal to . In order to take advantage of the symmetry, the coordinates have been placed in the vertical center of the capillary domain. In dimensionless coordinates, the wall located at the right-hand side is at the location while the one on the left-hand side is at (see Figure 2). In terms of the electrostatic forces, the system is defined by (2.1a) jointly with appropriate boundary conditions [3, 4]:

The HAM explicit analytical solution (2.31) depends on the auxiliary parameter Let denote the set of all values of which ensure the convergence of the series (2.31). According to Liao's proof [16], all of these series solutions must converge to the solution of the original problem (3.1a)-(3.1b). For a fixed in the range of plot the -curve, versus . Because the limit of all convergent series solutions (2.31) is the same for a given , there exists a horizontal line segment above the region . So, by plotting the -curve at a high enough order approximation, one can find an approximation , as shown in Figure 1 in the case of and . In Table 1, the value ensures the convergence of the series (2.31).

We can conclude that, for any given value of , we can always find a proper value of to ensure the convergence of the series of the solution

Let us consider the average residual errors of (3.1a): The average residual errors of the series (2.31) at different orders of approximations when are listed in Table 2. Obviously, as the order increases, the average residual errors decrease. This clearly indicates that our series solution (2.31) is indeed the solution of (3.1a)-(3.1b).

In a previous work [22], and for coordinates system as shown in Figure 1, the authors have obtained the approximate analytical solution of the model described by (3.1a)-(3.1b), and it is given by which is only valid in the range .

The proposed technique does not require small parameters in the equations under study. As a result, the technique completely eliminates the difficulty arising in the classical perturbation methods. The homotopy perturbation method (HPM) [23] was used to obtain an explicit analytical solution to the nonlinear Poisson- Boltzmann problem (3.1a)-(3.1b). The HPM solution is not accurate enough since it depends on an approximated form of (3.1a): and it is very difficult to obtain more than the first few terms of the series solution since the constant coefficients are not given by recursive formulas.

Figure 3 shows the electrostatic potential variation, along the transversal coordinate , predicted by the numerical solution of the complete Poisson-Boltzmann equation by Runge-Kutta method (RK4), Debye-Hückel approximation (3.3), 15th-order HAM approximation with , and approximation [4].

In particular, Figure 3 demonstrates that the use of the proposed approach, considerably improves the prediction of the electrostatic potential, , yielding values nearly identical with the numerical solutions. Obviously, we dare say that the best possible approximation with the numerical solution has been achieved by using HAM solution (2.31).

#### 4. Conclusion

Some models show nonlinear sources, such as, electrostatic potential. The general situation is that these models lack analytical or closed form solutions that are quite handy when coupling is required. In this study we developed a simple procedure to obtain an explicit and purely analytic solution for the prediction of the electrostatic potential; that is, the structure of the solution is known and the constant coefficients are given by recursive formulas and it is unnecessary to use any numerical methods to get any coefficients. This new approach embraces excellent results for a wide range of electrostatic potential values and can be highly adaptable to many situations.

#### Acknowledgment

The authors would like to express their great appreciation for the valuable comments and suggestions made by the anonymous reviewers.