# Hydrodynamic cell model: general formulation and comparative analysis of different approaches.

INTRODUCTION

For over a century, flowing of fluids through packed beds or porous diaphragms, as well as movement of particle swarms inside fluids, have been intensively studied by many researchers. For all these processes, mathematical modelling is conducted within the framework of a hydrodynamic boundary value problem. Formulation of such a problem includes the Stokes and continuity equations subject to boundary conditions at both the fluid/solid interface and the outer boundaries of the whole system. At the fluid/solid interface, the boundary conditions impose zero velocity of the fluid with reference to the particles. The outer boundary conditions impose the applied pressure difference and the fluid velocity through the system outer boundary.

Cell Model Approach and Spherical Cell Approximation

Due to the complex and, strictly speaking, uncertain geometrical structure of the solid phase, there are many difficulties in formulating and solving the above outlined boundary value problem. The Cell Model approach is a powerful method which enables one to circumvent the geometrical difficulties by analyzing the hydrodynamic field inside a representative cell which contains a group of particles, or a single particle, surrounded by the fluid and having given coordinates inside the cell.

The Cell Model approach is applicable for the system composed by the particles with fixed positions, as in the case of flows through packed beds or porous diaphragms. A more complex situation occurs while describing the systems containing particles that are free to move, as in the case of sedimentation. Considering behaviour of swarms of particles that are free to move, instead of the Cell Model approach, a great number of authors use the hydrodynamic modifications of the Statistical Mechanics methods, which have been developed in the classical studies of Batchelor (1972, 1974, 1982) and Batchelor and Wen (1982). Using these rigorous techniques lead to results which systematically differ from the predictions that can be obtained for the same systems by using different Cell Models. For example, according to the observation of Davis and Acrivos (1985), at low volume fraction of solid particles, [phi], these two types of approaches give different asymptotic forms for sedimentation velocity [u.sub.sed] ([phi]): [u.sub.sed] ([phi])-[u.sub.sed] (0)[[infinity].sup.1/3] (for different Cell Models) and [u.sub.sed] ([phi])-[u.sub.sed] (0) [phi](for Statistical Mechanics methods).

Numerous experimental data, which have been analyzed and summarized by Richardson and Zaki (1954), Barnea and Mizrahi (1973), Garside and A1-Dibouni (1977), Lynch and Herbolzheimer (1983), Davis and Acrivos, (1985), witness that the above mentioned two forms of the asymptotic behaviour describe two limiting cases of the experimentally observed behaviour of sedimenting particles. Thus, the Cell Model approach, which describes one of such limiting cases, [u.sub.sed] ([phi])-[u.sub.sed] (0)[[infinity].sup.1/3], is an important theoretical method for addressing systems with freely moving particles, not only the particles having fixed coordinates.

We will distinguish between two versions of the Cell Model approach. The first of them relates to the cases when porous diaphragms or particle swarms are considered as regular lattices formed by solid particles. For systems of this type, the elementary cell of the lattice serves as a representative cell for the whole system (Uchida, 1949; Sparrow and Loeffler, 1959; Hasimoto, 1959; Masliyah, 1973; Sangani and Acrivos, 1982; Zick and Homsy, 1982; Masliyah and van de Ven, 1986; Cheng and Papanicolaou, 1997; Maxey and Patel, 2001; etc.). Analysis of the hydrodynamic field within such a "regular" cell can give exact description of hydrodynamic parameters attributed to the ordered particulate systems. However, usually, the structure of the system is not known exactly, and, therefore, the results obtained for "regular" cell can only give information about some general properties of a real system that are independent of the specific geometry of the particle arrangement. At the same time, it is the polyhedron geometry of the lattice elementary cells that defines the major mathematical difficulties in solving the hydrodynamic boundary value problems inside the "regular" cells.

The second version of the Cell Model approach is based on the assumption that representative cell for a particle swarm or porous diaphragm is a sphere (for spherical particles) or cylinder (for fibrous particles). In some publications (Epstein and Masliyah, 1972; Dassios et al., 1995; etc.) the authors dealt with spheroid shapes of both the particles and the cells. In spite of its heuristic nature and certain inconsistencies, the second version of the Cell Model approach is widely used in Fluid Mechanics. Hereafter, we will refer to the second version of the Cell Model approach as the Spherical (or Cylindrical) Cell Approximation.

After the pioneering publications of Cunningham (1910), Frazer, (1926), Simha (1952), Happel (1957; 1958) and Kuwabara (1959), both the Spherical and Cylindrical Cell Approximations have been widely employed by different authors for obtaining effective hydrodynamic parameters of disperse systems. In particular, these methods have been used for predicting the hydraulic resistivity of a diaphragm and the sedimentation coefficient for particle swarms (Smith 1966, 1967; Gal-Or, 1970; Yaron and GalOr, 1971; Leclair, and Hamielec, 1971; Epstein and Masliyah, 1972; Masliyah, 1973; Neale and Nader, 1974; Neale and Masliyah, 1975; Mehta and Morse, 1975; Yu and Soong, 1975; Kvashnin, 1979; Dassios et al., 1995; Ziemer and Foerester, 1996; Chen et al., 1999; Datta and Deo, 2002; Filippov et al., 2006; etc.). In other publications, the Spherical Cell Approximation has been employed for determining the effective viscosity of suspensions (Ruiz-Reina et al., 2003; Rubio-Hernandez et al., 2004; Zholkovskiy et al., 2006). Over the three past decades, the Spherical Cell Approximation has been successfully used for addressing electrokinetic phenomena (Levine and Neale, 1974; Levine et al., 1976; Shilov et al., 1981, 2004; Zharkikh and Shilov, 1981; Ohshima, 1997a, b, 1998, 2000a, b; Dukhin et al., 1999a, b, c, 2002; Ding, and Keh, 2001; Carrique et al., 2001a, b, 2002; Lee et al., 2001, 2002; Dukhin and Goetz, 2002; Ruiz-Reina et al., 2003; Rubio-Hernandez et al., 2004; Masliyah and Bhattacharjee, 2006; etc.)

Within the framework of the Spherical Cell Approximation, a representative part for a particle swarm or diaphragm is assumed to be a sphere containing a fluid which surrounds a particle at the sphere centre (Figure 1). The ratio of the particle to cell volume is assumed to be equal to the particle volume fraction, [phi]. Consequently, the particle and cell radii, a and b, satisfy the equality

[phi] = [(a/b).sup.3] (1)

While addressing a hydrodynamic flow within the frameworks of the Spherical Cell Approximation, the Stokes and continuity equation are written for the fluid in the cell bulk. These equations are subject to the boundary conditions at the particle surface and at the cell outer boundary. At the particle surface, the boundary conditions have the same form, as those set at the surface of the "real" particles comprising the actual system. As for the boundary conditions at the cell outer boundary, their form is the focus of intensive discussion in the literature.

Different Versions of the Outer Boundary Conditions

In the literature dealing with the Spherical Cell Approximation, there are different versions of the outer boundary conditions (see Kvashnin, 1979; Ziemer and Foerester, 1996, for example). However, the most frequently used conditions have been proposed by Cunningham (1910), Happel (1958) and Kuwabara (1959). Each of these three authors imposed two outer boundary conditions. The first of the imposed conditions is common for all the three theories. Considering a system of particles placed in a uniform flow having velocity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], with reference to the particles (Figure 1), the common outer boundary condition can be written for the local velocity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], in the form:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] at the cell outer boundary (2)

Boundary condition (2) is insufficient for closing the above discussed hydrodynamic boundary value the problem. For closing, the problem formulation, Cunningham (1910), Happel (1958) and Kuwabara (1959) proposed three different versions of the additional outer boundary condition.

Strictly speaking, in the paper of Cunningham (1910), the proposed condition at the cell outer boundary was to equate the local velocity at the outer boundary, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], and the external flow velocity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Clearly, Equation (2) coincides with the Cunningham condition written for the normal velocity. Hence, to complete the Cunningham formulation, one should present the corresponding condition for tangent velocity. In the invariant form, such a condition can be written as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] Cunningham (1910) (3)

In the reference system linked to the particles, the additional boundary conditions proposed by Happel (1958) and Kuwabara (1959) can be represented as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

0 Kuwabara (1959) (5)

In Equation (4), [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is a viscous stress tensor given by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

where p is the local pressure, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the unit tensor, and (*) signifies a transposed tensor.

Using outer boundary condition (2) and one of conditions (3) to (5), one can solve the above hydrodynamic problem to address the hydrodynamic field within the spherical cell. Finally, using of the obtained distributions, one can determine the viscous force exerted on the particle, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Clearly, using different boundary conditions leads to different predictions of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The Hydrodynamic Interaction Factor

For convenience of the discussion, we will compare the predictions from different models in terms of the dimensionless coefficient L which is introduced by normalizing the hydrodynamic force exerted on the particle, F, with the Stokes factor, 6[pi][eta][u.sub.[infinity]]:

L = F/6[pi][eta]a[u.sub.[infinity]] (7)

Hereafter, we will refer to the coefficient L introduced by Equation (7) as the Hydrodynamic Interaction Factor (HIF). The HIF describes increase in the hydrodynamic force exerted on a given particle as compared with the case of a single particle placed in the flow with the same velocity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

[FIGURE 1 OMITTED]

An alternative, but equivalent, definition of the HIF can be given while considering migration of particles under influence of a given force. For such a case, the HIF describes the decrease in the migration velocity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] as compared with that for of a single particle, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]:

L = [u.sub.0]/[u.sub.d] (8)

The HIF given by Equation (8) can be interpreted as the ratio of the sedimentation velocities attributed to the case of a single particle and suspension, respectively.

Using Equation (7), one can obtain a convenient relationship which expresses the HIF, L, through the effective hydraulic resistivity of a diaphragm composed by the solid particles, R. We will consider a plane parallel diaphragm having linear dimension, [h.sub.k], along the kth axis of the Cartesian coordinate system shown in Figure 2. When the pressure difference, [DELTA]p, is applied between the diaphragm sides that are parallel to the coordinate plane [x.sub.2] - [x.sub.3], the effective hydraulic resistivity, R, can be defined, as:

r = [DELTA]p/[h.sub.1][u.sub.[infinity]] (9)

Making use of Equations (7) and (9) and realizing that [DELTA]P[h.sub.2][h.sub.3] = NF (N is the number constituent particles), one can represent the HIF as:

L = NF/6[pi][eta]a[u.sub.[infinity]]N = R/[R.sub.0] (10)

where [R.sub.0] is the asymptotic value of the hydraulic resistivity corresponding to the limiting case of zero volume fractions:

[R.sub.0] = 6[pi][eta]a[u.sub.[infinity]]N/[h.sub.1][h.sub.2] [h.sub.3][u.sub.[infinity]] = 9[eta]/[2a.sup.2][phi] (11)

where [phi] = 4[phi][[alpha].sup.3]/3[h.sub.1][h.sub.2][h.sub.3] is the volume fraction.

[FIGURE 2 OMITTED]

The dimensionless HIF given by any of equivalent definitions (7), (8) or (10), depends on the volume fraction, only, i.e., L = L([phi]). Due to the hydrodynamic interaction between the particles, L([phi]) is an increasing function of the volume fraction, i.e., L([phi]) [greater than or equal to] L(0) = 1.

Predictions from Different Models

Now, in terms of the HIF, L([phi]), we present the results, which have been obtained by Cunningham (1910), Happel (1958) and Kuwabara (1959) who used different boundary conditions given by Equations (3) to (5):

L ([phi]) = 4 (1 - [[phi].sup.5/3]/4 - 9[[phi].sup.1/3] + 10[phi] - 9 [[phi].sup.5/3] + 4[[phi].sup.2] Cunningham (1910) (12)

L ([phi]) = 2/3 x 3 + 2[[phi].sup.5/3]/2 - 3[[phi].sup.1/3] +3 [[phi].sup.5/3] - 2[[phi].sup.2] Happel (1958) (13)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] Kuwabara (1959) (14)

It should be noted that many authors, who used Kuwabara boundary condition (5) (Levine et al., 1976; Lee et al., 1978; Carrique et al., 2001b; Ohshima, 2000a; Datta and Deo, 2002; Dukhin and Goetz, 2002; etc.), report about another result than that given by Equation (14). In terms of the HIF, such an alternative result is represented as:

L ([phi]) = 1/ 1- 9/5 [[phi].sup.1/3] + [phi] - [[phi].sup.2]/5 (15)

Later in the present paper, such a surprising difference between the predictions from the same model, Equations (14) and (15), will be discussed in details. Now, we only mention that, while dealing with the same hydrodynamic field inside the cell, Kuwabara and the authors of result (15), use different method of obtaining the viscous force, F, represented in the definition of the HIF, Equation (7). In the original paper of Kuwabara (1959), the viscous force magnitude is determined by equating the power developed by the viscous force, [Fu.sub.[infinity]], and the overall dissipation rate that is obtained by integrating the local dissipation rate over the cell volume. Oppositely, while obtaining Equation (14), the viscous force is determined by integrating the tensor of viscous stresses over the particle surface or, equivalently, the cell outer boundary.

Thus, by using different outer boundary conditions (3) to (5), one obtains different expressions for the HIF given by Equations (12) to (15). The curves plotted in Figure 3 according to Equations (12) to (15) demonstrate the differences existing between the predictions from the Cunningham and Happel models. In particular, at [phi] = 0.3, Cunningham's value of the HIF substantially, by factor of about 2.8, exceeds the Happel result. Remarkably, the curves corresponding to the Kuwabara and Happel predictions given by Equations (13) and (14), respectively, nearly merge with each other. Note that such close predictions from the Kuwabara and Happel models are obtained by using different methods of determining the viscous force exerted on the particle. The above discussed alternative version of the Kuwabara model, Equation (15), gives a noticeably higher HIF than that obtained while using the Happel model: at [phi] = 0.3, the disagreement is about 30%.

Thus, the existing differences in the predictions demonstrate the necessity to choose the outer boundary condition giving a better approximation. However, while trying to make such a choice, we face substantial difficulties because all the outer boundary conditions given by Equations (2) to (6) have been proposed by their authors on the basis of qualitative arguments that look as equally reasonable. Thus, it is required to develop quantitative criteria which would be helpful in choosing a better version of the Spherical Cell Approximation.

Objective and Structure of the Paper

The objective of the present paper is to reconsider the basic principles of the Cell Model approach and the Spherical Cell Approximation by primarily focusing on the boundary conditions employed at the cell outer boundary. To this end, in the second section, General Formulation, we will start our analysis with the general problem formulation, which describes a flow through a diaphragm shown in Figure 2. In the third section, Representative Cell of the General Type, we will consider representative cell of the general type and obtain relationships interrelating the hydrodynamic quantities at the cell outer boundary with applied pressure difference, [DELTA]p and the external flow velocity, [u.sub.[infinity]]. The relationships to be derived will be based on the mass, momentum and dissipation balance conditions which should be satisfied when a given volume is a representative part of diaphragm considered in the General Formulation section.

In the fourth section, using the Spherical Cell Approximation, we will reformulate the general problem considered in the second section. In particular, we will obtain the cell outer boundary conditions by specifying the general relationships derived in the third section. While reformulating the problem, we will use two basic principles of the Spherical Cell Approximation:

1. The Spherical Cell, whose radius satisfies Equation (1), is a representative part of the diaphragm;

2. inside the Spherical Cell, the angular symmetry of the hydrodynamic field is the same as in the case of a single spherical particle placed in a uniform flow.

On this basis, in the fifth section, Analysis of Predictions and Errors from Different Models, we will re-derive the results reported in the literature and analyze the errors associated with using different boundary conditions. In the sixth section, Discussion, we will consider the common and different features of the hydrodynamic and electrodynamic versions of the Spherical Cell Approximation and discuss which version of the Cell Model gives the best approximation in describing the pressure and electrically driven flows through a diaphragm and sedimentation of particles.

[FIGURE 3 OMITTED]

GENERAL FORMULATION

Our analysis is aimed to obtain the HIF, L, by using Equation (10) which expresses the function L([phi]) through the hydraulic resistivity of a diaphragm. To this end, we will consider a diaphragm of the thickness [h.sub.1] (Figure 4). Other linear dimensions of the diaphragm, [h.sub.2] and [h.sub.3], are assumed to be much bigger than [h.sub.1]. Accordingly, we will neglect all the edge effects whose contribution is of order of O([h.sub.1]/[h.sub.2,3]). To justify external boundaries of the system shown in Figure 4, we introduce two hypothetical parallel planes, A and B, separated by equal distances, d, from the diaphragm.

We will analyze hydrodynamic flow through a macroscopically homogeneous diaphragm which, at the microscopic level, is built up by solid particles. The macroscopic homogeneity means that, at any place inside the diaphragm, one can single out a representative cell whose largest dimension, b, is much smaller than the diaphragm thickness, [h.sub.1], and the distance, d, between the diaphragm and each of the hypothetical panes, A and B. (Figure 4). Hereafter, a given part of the diaphragm is considered as a representative cell for the diaphragm when averaging the distributions of hydrodynamic quantities over the volume of such a part leads to the same results as the averaging over the whole diaphragm volume. At this stage of the analysis, we deal with arbitrary shaped representative cells containing any number of constituent particles. In summary, the hierarchy of the above introduced length scale parameters can be expressed with the help of the following sequence of inequalities:

b [much less than] d [much less than] [h.sub.1] [much less than] [h.sub.2,3] (16)

Now, we consider a pressure driven creeping flow of a fluid through the diaphragm. The flow is described by using the Stokes and continuity equations attributed to the fluid volume which is bound by the planes A and B and the surfaces of all the particles forming the diaphragm (Figure 4):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (18)

[FIGURE 4 OMITTED]

Equations (16) and (17) are subject to the boundary conditions at the particle surfaces and at the planes A and B. Consequently, we set zero local velocity:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] at the particle surfaces (19)

At each of the planes A and B we impose constant pressure

p = [DELTA]p at the plane A (20)

p = 0 at the plane B (21)

For the general case, governing Equations (17) and (18) subject to boundary conditions (19) to (21), do not form a closed problem formulation: it is necessary to impose additional boundary conditions at the planes A and B (Figure 4). For the present problem, such additional conditions should reflect the macroscopic homogeneity of the membrane. We will consider a limiting version of the first of inequalities (16): a vanishingly small representative cell of the diaphragm b/d [right arrow] 0. It can be shown that, at b/d [right arrow] 0, the expression for the local velocity field, which satisfies Equations (17) to (21), contains the terms diverging at the planes A and B. To eliminate these terms, one should impose the following boundary condition:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] converges at the planes A and B (22)

Thus, governing Equations (17) and (18) subject to boundary conditions (19) to (22) form a closed problem formulation for obtaining the local velocity and pressure, distributions, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], within the fluid which is bound by the planes A and B and the surfaces of the particles constituting the diaphragm.

When the first and second of inequalities (16) are valid, the solution of the problem given by Equations (17) to (22), has an important property which is dictated by the macroscopic homogeneity of the diaphragm. At sufficiently large distances from the external surfaces of the diaphragm, i.e., at the much larger distances than the dimension of the representative cell, b, all the distributions become dependent on the coordinate [x.sub.1], only (Figure 4). For such a case, the liquid velocities in the vicinities of the planes A and B, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], are the constant vectors normal to the planes. Due to the continuity of the local velocity, Equation (17), the velocities [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] coincide with the external flow velocity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (Figure 4). Consequently, one can determine [u.sub.[infinity]], as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] at the planes AA or BB, (23)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the unit vector attributed to the axis [x.sub.1] of the Cartesian coordinate system shown in Figure 4. Thus, after obtaining the function [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], as a solution of the hydrodynamic boundary value problem given by Equations (17) to (22), one can determine [u.sub.[infinity]] with the help of Equation (23). The obtained velocity magnitude, [u.sub.[infinity]], is substituted into Equation (9) for obtaining the hydraulic resistivity R. Finally, making use of the normalization given by Equations (10) and (11) one obtains the HIF, L, for the diaphragm.

As stated above, the formulated hydrodynamic boundary value problem describes the hydrodynamic flow through a diaphragm that can be considered as a system of the representative cells. Provided that the cell dimension, b, satisfies inequalities (16), the above formulation is valid for arbitrary shaped representative cells which can contain any numbers of the particles. Analysis of hydrodynamic field inside a single representative cell instead of the whole system enables substantial simplifications in consideration. In the next section we will derive relationships that are important for such an analysis.

REPRESENTATIVE CELL OF THE GENERAL TYPE

While considering hydrodynamic field inside the representative cell, the governing equations and the boundary conditions at the particle surfaces are given by Equations (17) to (19). For closing the problem formulation, it is necessary to impose boundary conditions at the cell outer boundary. The outer cell boundary conditions should be Cell Model equivalents of boundary conditions (20) to (23) that are imposed at the outer boundaries of the actual system, planes A and B (Figure 4). In the present section, we will derive some general integral relationships which form a basis for obtaining the outer cell boundary conditions.

Mass Balance

Now we will derive a relationship enabling one to connect the external flow velocity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], and the velocity distribution inside the cell. To this end, we will consider two equalities expressing the mass balance for an incompressible fluid flowing through the diaphragm

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (24)

where [S.sub.C] ([x.sub.1]), is a hypothetic plane which is drawn inside the diaphragm to be parallel to the planes A and B. Position of such a parallel plane is defined by its Cartesian coordinate, [x.sub.1]. The plane [S.sub.C] ([x.sub.1]) has the same area, S = [h.sub.2] [h.sub.3], as the planes A and B. The unit normal vector for such a plane, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], is chosen such that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Equation (24) can be obtained by considering a part of the diaphragm enveloped by a surface which includes the planes A and C and a side surface. Consequently, making use of continuity Equation (17) and the Gauss theorem leads to Equation (24) provided that the latter of inequalities (16) is satisfied.

Let us consider the following identical transformations of Equation (24):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (25)

where the operator

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (26)

signifies averaging over a volume V.

Clearly, the averaging over the diaphragm and representative cell volumes should give the same results.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (27)

Consequently, combining Equations (25) and (27) yields

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (28)

Rewriting the latter expression according to Equation (27) one obtains

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (29)

Finally, we transform the above integral over the cell volume into the integral over the cell surface. To this end, we will use the equality

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (30)

which is valid since [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], Equation (18). Combining Equations (29) and (30) and using the Gauss theorem we arrive at the final equation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (31)

Thus, Equation (31) enables one to express the external flow velocity, [u.sub.[infinity]], through a given velocity distribution, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], at the cell outer boundary.

Momentum Balance

It should be noted that measurable quantities do not depend on the reference point for the pressure. Taking this into account, boundary conditions (20) and (21) can be employed for interrelating between the hydrodynamic field inside the cell and the pressure difference applied to the diaphragm. Consequently, we consider the total force exerted on the fluid volume which is bound by the closed complex surface, S, which includes planes A and B, the particle surfaces and the side surface of the system shown in Figure 4. For a stationary creeping flow, such a force should be zero, i.e.:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (32)

In Equation (32), [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the outward normal vector for the surface S. While considering the separate contributions of the integrals over the surfaces constituting the complex surface S, one can neglect by the contribution of the integral over the side surface of the plane-parallel system shown in Figure 4. Such a contribution can be ignored under the condition given by the latter of inequities (16). As for the integrals over the planes A and B, their forms are defined by the uniformity of the hydrodynamic flow in the vicinity of the planes A and B. Recall that this uniformity exists when inequalities (16) are valid. For a uniform flow, the viscous stress tensor given by Equation (6) takes a diagonal form. Consequently, Equation (32) can be rewritten as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (33)

In the second term of Equation (33), the summation is conducted over all the particles forming the diaphragm. Accordingly, the integration under the sign of sum is conducted over the surface of the qth particle. In Equation (33), [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are the outward normal unit vectors attributed to the planes A and B and the qth particle surface, [S.sub.q], respectively.

Combining Equations (20), (21) and (33), and realizing that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (Figure 4), one can rewrite Equation (33) in the following form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (34)

where V = [h.sub.1] [h.sub.2] [h.sub.3] is the diaphragm volume. According to inequalities (16), we neglected, by the thickness of the gap between the diaphragm and the planes A and B, d, as compared with the diaphragm thickness, [h.sub.1]. Thus, Equation (34) interrelates the mean pressure gradient, [DELTA]P/[h.sub.1], with the tensor of viscous stresses, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], at the particle surfaces.

Since the problem formulated in the second section, General Formulation, describes a diaphragm being a set of the representative cells, one can conduct the summation in Equation (34) in two steps: (i) to summate the surface integrals for the particles which belong to a single representative cell; and (ii) to multiply the sum obtained for a single cell by the number of the constituent cells, [N.sub.cell]. While dealing with integrals of the viscous stress tensor, without losing generality, one can replace the sum of the integrals over the particle surfaces by the similar integral over any surface enveloping the volume containing the particles. Considering the external boundary of the representative cell, [S.sub.cell], as such an enveloping surface and following the above described scheme of summation, Equation (34) is transformed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (35)

In Equation (35), [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the unit outward normal vector for the cell boundary and [V.sub.cell] = V/[N.sub.cell] is volume of the individual representative cell. Equation (35) is of great generality since it describes a macroscopically homogeneous diaphragm composed by representative cells having arbitrary shapes and containing any number of the particles. The importance of Equation (35) is that it enables one to compute the mean pressure gradient, [DELTA]P/[h.sub.1], while considering the hydrodynamic field inside a single representative cell, only.

Now, by combining Equations (6) and (35), we will derive an explicit expression for the mean pressure gradient, [DELTA]P/[h.sub.1], through the velocity and pressure distributions, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], inside the cell. Accordingly, making use of Equation (6), we will transform the expression under integral in Equation (35), as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (36)

Let us now consider the second term on the right-hand side of Equation (36). Using the convention of summation, after minor transformation, one obtains

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (37)

where [u.sub.k] are the coordinates of the vector [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] in the Cartesian coordinate system shown in Figure 4. Consequently, integrating the latter expression in Equation (37) over the cell surface yields

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (38)

Thus, the second term on the right-hand side of Equation (36) does not contribute into the integral in Equation (35).

After some straightforward transformations, the third term on the right-hand side of Equation (35) takes the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (39)

Finally, combining Equations (34) to (38) we arrive at the following result

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (40)

Thus, Equation (40) interrelates the mean pressure gradient, [DELTA]P/[h.sub.1], and the velocity and pressure distributions, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], inside the representative cell of the general type.

Dissipation Balance

Let us consider a hypothetic closed surface S* which envelops a volume, [V.sub.*], of a fluid containing motionless particles. The mean volume viscous dissipation rate inside such a volume, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], is determined by using the equality

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (41)

The expression for the overall dissipation rate, which is given by the integral over the surface S* in Equation (41), is obtained by integrating the local dissipation rate over the volume [V.sub.*] with further transformation of the integral over a volume into the surface integral. Details of this derivation can be found elsewhere (for example, Happel and Brenner, 1965).

Now, we will use Equation (41) twice by identifying [V.sub.*] as: (i) the volume, V, between the planes A an B (Figure 4); and (ii) the representative cell volume. Combining Equations (6), (20), (21), (23) and (41), the first of the above identifications leads to the following condition

[<w>.sub.V] = [DELTA]p/[h.sub.1] [u.sub.[infinity]] (42)

While obtaining Equation (42), we took into account that, under conditions given by inequalities (16), one deals with a uniform flow velocity in the vicinity of the planes A and B. As well, under the same conditions, we neglected by volumes of the gaps between the planes and the diaphragm and by the contribution of the integral over the side surface of the system shown in Figure 4.

Now, considering V* as the volume of the representative cell, [V.sub.cell], we rewrite Equation (41) as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (43)

Realizing that the dissipation rates averaged over the whole system and representative cell volumes coincide, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], one obtains from Equations (42) and (43) that:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (44)

Thus, Equation (44) establishes equality between the mean dissipation rate, [u.sub.[infinity]] [DELTA]p/[h.sub.1], and the same quantity obtained by averaging the local dissipation rate over the representative cell volume and transforming the result into the integral over the cell surface.

Finally, we represent Equation (43) in an equivalent form by combining the expression under integral with Equation (6):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (45)

To represent the second and the third terms on the right-hand side of Equation (45) in a convenient form, we consider the following identical transformations

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (46)

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (47)

Combining Equations (44) to (47) we arrive at the final expression

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (48)

Thus, Equation (48) yields a relationship between the hydrodynamic field inside the cell and mean dissipation rate which is obtained as a product of the mean pressure gradient and the external flow velocity, [u.sub.[infinity]] [DELTA]p/[h.sub.1]

Generalized Kuwabara Condition

In the previous subsection, Dissipation Balance, the dissipation balance equation was derived by considering a volume, [V.sub.*], enveloped by a surface [S.sub.*] and particles surrounded by fluid inside the volume. Now, by considering the same volume, we will derive another important equality.

Within the fluid, the pressure [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is a continuous function of coordinates. Therefore, for the pressure distribution within the fluidic part of the volume, V*, making use of the gradient theorem, the following identity can be proved

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (49)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are the outward normal unit vectors attributed to the surface [S.sub.*] and the qth particle surface, respectively. In Equation (49), the summation is conducted over all the particles inside the enveloped volume. Each of the integrals under the sign of sum is taken over the qth particle surface.

Within the frameworks of the present analysis, without losing generality, one can consider the solid particles as a limiting case of liquid particles having infinitely high viscosity. Accordingly, we will consider the pressure inside the qth particle, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Within the qth particle volume, the pressure, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], is a continuous function of the coordinates. The latter enables one to use the gradient theorem with respect to [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (50)

Combining Equations (49) and (50), after minor transformations, we arrive at the following identity

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (51)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the pressure gradient averaged over the whole volume [V.sub.*]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (52)

Thus, Equation (51) gives an identity which is valid for the pressure distribution inside an arbitrary volume, V*, which contains particles (or drops) surrounded by a fluid and enveloped by the surface S* drawn inside the fluid.

Now, we will use Equation (51) twice. First, we identify S* as the surface consisting of the planes A and B and the side surface of the system shown in Figure 4. According to the second and third of inequalities (16), while applying Equation (51), one can neglect by the volume of the gaps between the planes and the diaphragm and by the contribution of the integrals over the side surface. Consequently, one obtains

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (53)

where V = [h.sub.1] [h.sub.2] [h.sub.3] is the diaphragm volume. Recalling that the diaphragm is assumed to be a set of representative cells, the summation in Equation (53) can be conducted in a similar manner to that described above Equation (35). Accordingly, we summate the integrals over the surfaces of the particles which belong to a single representative cell. After that, we multiply the obtained sum by the numbers of the representative cells, [N.sub.cell]. Consequently, taking into account that V/[N.sub.cell] = [V.sub.cell], Equation (53) can be rewritten as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (54)

where the summation is conducted over the particles included into the representative cell.

The second use of Equation (51) is to identify [S.sub.*] as the outer boundary of the representative cell. Such an approach allows Equation (51) to take the form:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (55)

Let us now compare Equations (54) and (55). The right-hand sides of these equations contain the common second term. As well, according to Equation (27) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and, hence, the first terms on the right-hand sides of Equations (54) and (55) are equal, as well. Consequently, by equating the left-hand sides of Equations (54) and (55), we obtain the following relationship:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (56)

Thus, Equation (56) gives an expression for the mean pressure gradient, [DELTA]P/[h.sub.1], through the pressure distribution, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], at the representative cell outer boundary. It should be noted that Equation (56) is of great generality and describes the diaphragm constituted by arbitrary shaped representative cells which contain any number of particles.

In the Momentum Balance subsection, using the force balance, we derived Equation (40) which is of the same generality as Equation (56). Remarkably, Equations (40) and (56) differ from each other due to the presence of the second term in the square brackets on the right-hand side of Equation (40). Hence, Equations (40) and (56) can simultaneously be satisfied provided that the second term in the square brackets in Equation (40) becomes zero:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (57)

Thus, Equation (57), imposes a condition for the vorticity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], at the cell outer boundary. Note that Equation (57) is valid for a representative cell of the general type. As it will be shown in the next section, for the spherical cell containing a single particle, Equation (3) is transformed into the Kuwabara (1959) condition given by Equation (5). Therefore, Equation (57) can be considered as a generalization of the Kuwabara condition for the case of arbitrary shaped cell containing any number of particles.

In summary, in these subsections, we derived four relationships given by Equations (31), (40), (48) and (57) that impose certain conditions for the integrals over the outer cell boundary. Utilizing symmetry properties of a representative cell under consideration, one can transform the derived four equations into the cell outer boundary conditions. For the case of Spherical Cell Approximation, such a transformation will be conducted in the next section where we reformulate the general boundary value problem discussed in the second section, General Formulation.

SPHERICAL CELL APPROXIMATION: REFORMULATION OF GENERAL PROBLEM

The present section will describe the hydrodynamic field inside the Spherical Cell containing a particle surrounded by a fluid (Figure 1). Since the Spherical Cell is assumed to be a representative part for a diaphragm composed by spherical particles, the cell and particle radii, b and a, should satisfy Equation (1).

For addressing the hydrodynamic field inside the Spherical Cell, we will consider the boundary value problem including Stokes and continuity governing Equations (17) and (18). At the particle surface, the governing equations are subject to the boundary condition (19). For closing the problem formulation, we will impose conditions at the cell outer boundary by using four relationships (31), (40), (48) and (57) derived in the third section, Representative Cell of the General Type.

While using the Spherical Cell approximation, an important simplification in the hydrodynamic analysis is associated with the assumption that the angular symmetry of hydrodynamic field inside the cell is the same as in the case of a single particle placed in a uniform flow. Using such an assumption, we will specify Equations (17) to (19) and general relationships (31), (40), (48) and (57) which yield the outer cell boundary conditions.

Hydrodynamic Field inside the Cell for Uncertain Outer Boundary Conditions

Within the Spherical Cell, we will consider the velocity and pressure distributions using the spherical coordinate system shown in Figure 5. Accordingly, the spherical coordinates, r, [theta] and [phi] are interrelated with the Cartesian coordinates, [x.sub.k], employed in the third section, Representative Cell of the General Type, as:

[x.sub.1] r cos ([theta]) (a)

[x.sub.1] r sin ([theta]) cos ([phi]) (b)

[x.sub.1] r sin ([theta]) sin ([phi]) (c) (58)

The unit vectors attributed to each of the coordinate systems are interrelated as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (59)

Assume that the hydrodynamic field inside the cell has the same angular symmetry as that in the case of a single particle placed in a uniform flow whose velocity is directed in parallel to the axis [x.sub.1]. Such an assumption dictates the following angular dependencies for the pressure and velocity distributions

p = [PI] (r) cos ([theta]) (60)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (61)

where, using continuity Equation (18), [U.sub.r](r) and [U.sub.[theta]](r) are expressed as:

[U.sub.r](r) = - 2/[r.sup.2] Y(r) (62)

[U.sub.[theta]](r) = 1/r dY(r)/dr (63)

The functions Y(r) and [PI](r) are determined as solutions of an ordinary differential equation set which is obtained while substituting the distributions given by Equations (60) to (63) into governing Equations (17) and (18). Consequently, following a standard scheme, which can be found elsewhere (see, for example, Happel and Brenner, 1965), one obtains the functions [PI](r) and Y(r) that are dependent on four integration constants. Two of such constants are eliminated by combining Equations (60) to (63) with boundary condition (19) that imposes zero velocity at the particle surface, i.e., at r = a. The functions Y(r) and [PI](r) obtained in such a manner depend on two uncertain constants, A and B. These functions can be represented as:

Y(r) = B [A [(r/a).sup.4] + [(r/a).sup.2] - (5A + 3) 4/2a + (3A + 1) a/2r] (64)

[PI](r) = - [eta] d/dr [[d.sup.2]Y/[dr.sup.2] - 2Y/[r.sup.2] (65)

Thus, Equations (60) to (65) describe hydrodynamic field inside the spherical cell shown in Figure 5 for uncertain outer boundary conditions.

Outer Boundary Conditions

Now, we consider how the general relationships derived in the third section, Representative Cell of the General Type, and given by Equations (31), (40), (48) and (57) are specified while using the Spherical Cell approximation. Inspecting Figure 5 we will use the following equalities for transformation of Equations (31), (40), (48) and (57)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (66)

[FIGURE 5 OMITTED]

Mass balance condition

Combining the mass balance condition given by Equation (31) with Equations (58a) and (66) yields

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (67)

Substituting into the above integral the velocity [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] given by Equation (61) leads to a simple result

[U.sub.r](b) = [u.[infinity]] (68)

Thus, Equation (68) gives the Spherical Cell version of general condition (31).

It should be noted that Equation (68) is equivalent to Equation (2) which, as stated in the Introduction section, is a common boundary condition employed in the models of Cunningham (1910), Happel (1958) and Kuwabara (1959). Such an equivalency becomes obvious when we multiply both sides of Equation (68) by cos([theta]) and consider the following chain of equalities based on Equation (61):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (69)

Hence, the above consideration, which includes the analysis given in the Mass Balance section and Equations (62) to (64), obtained in the present section, can be considered as a derivation of boundary condition (2) which thus has the meaning of the mass balance condition.

For convenience, boundary condition (68) can be rewritten in terms of the function Y(r). Accordingly, making use of Equation (62) one obtains

Y(b) = - 1/2 [u.sub.[infinity]] [b.sup.2] (70)

Thus, by using mass balance condition (70), the value of the function Y(r) at the cell outer boundary is interrelated with the filtration velocity, [u.sub.[infinity]].

Momentum balance condition

Prior to transforming Equation (40), which yields the momentum balance condition, we represent the vorticity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], in the spherical coordinate system shown in Figure 5. Using Equations (61) to (63), the vorticity is expressed in the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (71)

Using Equation (59a) the vector products under integrals in Equation (40) are transformed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (72)

Consequently, combining Equations (40), (60), (67), (71) and (72) one obtains:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (73)

After integration, we arrive at the expression

[PI](b)/b - 2/[b.sup.2] [eta] [[d.sup.2]Y/[dr.sup.2](b) - 2/[b.sup.2] Y(b)] = [DELTA]P/[h.sub.1] (74)

Thus, Equation (74) gives the Spherical Cell version of momentum balance condition.

Dissipation balance condition

The Spherical Cell versions of the dissipation balance condition are obtained by combining Equations (48), (60) to (63) and (66). After straightforward transformations, one obtains:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (75)

Thus, Equation (75) yields the dissipation balance condition for the Spherical Cell approximation

Kuwabara condition

This derivation repeats all the steps which leads to obtaining the second term on the left-hand side of Equation (74). Accordingly, combining Equations (50), (60), (76), (71) and (72), the generalized Kuwabara condition (57) is transformed as:

[d.sup.2]Y/[dr.sup.2](b) - 2/[b.sup.2] Y(b) = 0 (76)

Comparing the above equation with the expression for vorticity given by Equation (71), we obtain a confirmation that Equation (76) is equivalent to the Kuwabara condition given by Equation (5).

Different Sets of the Outer Boundary Conditions

In the Outer Boundary Conditions subsection, we specified the general relationships given by Equations (31), (40), (48) and (57) and obtained four boundary conditions to be imposed at the Spherical Cell outer boundary, Equations (70) and (74) to (76). Note that the velocity and pressure distributions obtained in Hydrodynamic Field inside the Cell for Uncertain Outer Boundary Conditions subsection for uncertain outer boundary conditions depend on two unknown constants, A and B, represented in Equation (64). Hence, two outer boundary conditions are required for eliminating these constant. The third outer condition is necessary to obtain the hydraulic resistivity. Thus, for addressing the hydrodynamic field inside the cell and determining the hydraulic resistivity, it is sufficient to set three outer boundary conditions, only.

The above discussion shows that, for obtaining the hydraulic resistivity, R, and the HIF, L, one should choose only three from four conditions (70) and (74) to (76). While choosing a set of three conditions, the fourth, "spare", condition is not employed in the problem formulation. For diaphragms that are constituted by "regular" representative cells, as in the case of regular lattices, the unutilized condition would be satisfied, automatically, provided that the three employed conditions are satisfied. Since the Spherical Cell is an approximate version of the representative cell, Equations (70) and (74) to (76), should be considered as approximate conditions. Therefore, a solution obtained by using any set of three from four boundary conditions (70) and (74) to (76) can be contradictory with respect to the fourth, omitted, condition. In such a situation, one can expect discrepancies between the predictions from the models that use different sets of three from four outer boundary conditions (70) and (74) to (76). Below, we give brief descriptions for each of such models.

Model 1: Mass-Momentum Balance

This model includes the mass and momentum balance conditions given by Equations (70) and (74), respectively, and Kuwabara condition (76). Accordingly, the dissipation balance condition, Equation (75), is omitted. Hence, while using the Mass Momentum Balance Model, one can expect a violation of dissipation balance condition (75).

Model 2: Mass-Dissipation Balance

Instead of momentum balance condition (74) included in the above discussed model, the Mass-Dissipation Balance Model employs dissipation balance (75). Thus, Model 2 deals with the mass and dissipation balance condition given by Equations (70) and (75) and with Kuwabara condition (76). Using this model, one can expect a violation of momentum balance (74).

Model 3: Momentum-Dissipation Balance

Within the frameworks of this model, one deals with boundary conditions (74) to (76) whereas the mass balance condition given by Equation (70) is omitted. Accordingly, using the Momentum-Dissipation Balance Model may lead to a violation mass balance condition (70).

Models 4 and 5: Mass-Momentum-Dissipation Balance

Now, we consider simultaneous use of all the balance conditions given by Equations (70), (74) and (75). It should be noted that simultaneous use of Equations (70), (74) and (75) leads to two different models. The latter becomes clear while analyzing the structure of the expression which is obtained by substituting the velocity, [u.sub.[infinity]], given by Equation (70) and the mean pressure gradient, [DELTA]p/[h.sub.1], given by Equation (74) into the right-hand side of dissipation balance (75):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (77)

Substituting Y(r) given by Equation (64) into Equation (77) leads to the quadratic equation with respect to the unknown constant A represented in Equation (64). As for the multiplicative constant B from Equation (64), it is cancelled out while using the above substitution. Two roots of the quadratic equation obtained in such a manner define two different values of the coefficient A, and thus two different hydrodynamic fields and hydraulic resistivities.

Remarkably, two different values of A, which are the roots of the quadratic equation obtained by combining Equations (64) and (77), can be determined in an alternative manner by using either the Cunningham (1910) or Happel (1958) boundary conditions given by Equations (3) and (4), respectively. To demonstrate this for the case of Cunningham condition (3), we consider the general form of the dissipation balance which can be obtained by combining Equations (35) and (44)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (78)

As stated in Different Versions of the Outer Boundary Conditions subsection, Equations (2) and (3) are equivalent to the condition proposed in the original paper of Cunningham (1910) who equated the local velocity at the cell outer boundary and the external flow velocity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Thus, using Figure 4, two conditions (2) and (3) can be rewritten as a single equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] at the cell outer boundary (79)

While substituting the velocity [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] given by Equation (79) into the integral on the left-hand side of Equation (78), one obtains an identity which is valid for an arbitrary shaped cell. Being a particular case of Equation (78), Equation (77) should be valid for any hydrodynamic field satisfying Equation (79) at the cell outer boundary. Hence, determining the unknown integration constant, A, from condition (79) yields one of two roots of the quadratic equation which is obtained by combining Equations (64) and (77).

It is important that dissipation balance (78) is also valid when the Happel condition given by Equation (4) is satisfied. To show this, we consider the following transformation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (80)

The first of the above equalities is an identity which, using Happel condition (4), leads to the final expression in (80). Now, combing Equations (2) and (80) and taking into account that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], the left-hand side of Equation (78) is transformed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (81)

Above, starting with the left-hand side of Equation (78), we arrived at its right-hand side given by the latter expression in (81). Thus, simultaneous use of Happel condition (4) and the mass and momentum balance conditions given by Equations (70) and (74) provides the validity of dissipation balance (78).

In summary, combining the mass, momentum and dissipation balance conditions leads to quadratic equation for obtaining the unknown constant A. The roots of this quadratic equation can be determined as solutions of two linear equations that are obtained by specifying the Cunningham (1910) and Happel (1958) conditions given by Equations (3) and (4) for the case of the Spherical Cell Approximation.

[FIGURE 6 OMITTED]

It should be noted that the above observation regarding the satisfaction of the dissipation balance while using the Happel outer condition is specific for the case of a uniform external flow which is considered in the present paper. As shown by Zholkovskiy et al. (2006), in the case of a purely shear external flow, the dissipation balance is violated when the Happel condition is satisfied. At the same time Zholkovskiy et al. demonstrated that, for the shear external flow, the dissipation balance is satisfied while using the Simha (1952) outer boundary condition. Remarkably, the Simha outer boundary condition, which imposes the undisturbed velocity of the external shear flow at the outer cell boundary, is somewhat similar to the Cunningham (1910) condition.

To obtain the Spherical Cell version of Cunningham condition (3) we combine Equations (2), (3), and (61) to (63). Consequently, realizing that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], after some straightforward transformations, Cunningham condition (3) takes the form

dY/dr (b) - 2/b Y(b) = 0 (82)

For obtaining the the Spherical Cell version of Happel condition (4), we consider the following transformations

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (83)

At the outer cell boundary, where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], the first term in the quadratic brackets in the latter expression in (83) is transformed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (84)

Combining Equations (61), (63), (71), (83) and (84) one obtains

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (85)

Using Equation (85), Happel condition (4) takes the form

[d.sup.2]/[dr.sup.2] (b) - 2/b dY/dr (b) + 2/[b.sup.2] Y(b) = 0 (86)

Thus, Equations (82) and (86) give the Spherical Cell versions of Cunningham (1910) and Happel (1958) outer boundary conditions, respectively.

The above analysis demonstrates the existence of two models, 4 and 5, for which the mass, momentum and dissipation balance conditions are satisfied, simultaneously. Each of these two models uses Equations (70) and (74) expressing the mass and momentum balance, respectively. Additionally, Models 4 and 5 use the Cunningham and Happel conditions, respectively. The Spherical Cell versions of the Cunningham and Happel conditions are given by Equations (82) and (86), respectively.

ANALYSIS OF PREDICTIONS AND ERRORS FROM DIFFERENT MODELS

In the present Section, we will use different sets of three outer boundary conditions attributed to each of the five models itemized in the Different Sets of the Outer Boundary Conditions subsection. By using these conditions, we will eliminate two unknown constants, A and B represented in the expression for the function Y(r) given by Equation (64) and determine the hydraulic resistivity, R. For each of the model we will predict the HIF, L, which is obtained by normalizing the hydraulic resistivity, R, according to Equations (10) and (11). Since each of the models violate one of the four conditions given by Equations (31), (40), (48) and (57), we will determine the respective errors originating from such violations.

Model 1: Mass-Momentum Balance

According to the Different Sets of the Outer Boundary Conditions subsection, this model includes the mass balance, momentum balance and Kuwabara conditions given by Equations (70), (74) and (76), respectively. Combining Equation (64) with the Kuwabara boundary condition (76) we obtain one of two unknown constants represented in Equation (64):

A= - 3/5 [phi]/[phi] + 2 (87)

Instead of determining the remaining unknown constant in Equation (64), B, we will combine mass and momentum balance conditions (70) and (74) with the expression (9) which yields the hydraulic resistivity given. Consequently, taking into account the Kuwabara boundary condition (76), one obtains

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (88)

were [R.sub.1] signifies the hydraulic resistivity determined within the frameworks of Model 1. Finally, combining Equations (10), (11), (64), (65), (87) and (88) we arrive at the expression for the HIF, attributed to Model 1, [L.sub.1]

[L.sub.1] = - [a.sup.3]/9[[phi].sup.4/3] Y(b) [{d/dr [[d.sup.2]Y/[dr.sup.2] - 2Y/[r.sup.2]]}.sub.r=b] = 5/(5 + 6[[phi].sup.1/3] + 3[[phi].sup.2/3] + [phi]) [(1 - [[phi].sup.1/3].sup.3](89)

Note that the HIF becomes independent of the multiplicative constant, B, which is cancelled out while substituting (64) into (89).

As it could be expected, the final result in (89) is equivalent to the function L([phi]) given by Equation (15). Thus, in terms of the classification introduced in the Different Sets of the Outer Boundary Conditions subsection, the approach, which leads to result (15) (Levine et al., 1976; Lee et al., 1978; Carrique et al., 2001b; Ohshima, 2000a; Datta and Deo, 2002; Dukhin and Goetz, 2002; etc.), employs the Mass-Momentum Balance Model. According to discussion of the Different Sets of the Outer Boundary Conditions subsection, using such a model can result in a violation of the dissipation balance condition whose Spherical Cell version is given by Equation (75).

In order to evaluate the error originating from the violation of the dissipation balance, we introduce a dimensionless parameter, [[delta].sub.1], as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (90)

The numerator of Equation (90) gives a difference between the right- and left-hand sides of dissipation balance condition (44) where the external flow velocity, [u.sub.[infinity]], and the mean pressure gradient, [DELTA]p/[h.sub.1], are replaced by relevant expressions given by Equations (31) and (35). In the denominator of Equation (90), we placed one of two dissipation rates which turn out to be of lower magnitude.

Combining Equation (75), which is the Spherical Cell version of Equation (44), with Equations (70), (74) and (76), Equation (90) is rewritten as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (91)

Substituting [PI](r) and Y(r) given by Equations (64) and (65) into the above expression and taking into account Equation (87), the normalized error, [[delta].sub.1], is expressed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (92)

Thus, Equation (92) represents the parameter [[delta].sub.1] describing the violation of the dissipation balance while using the Mass-Momentum Balance Model. The function [[delta].sub.1] ([phi]) described by Equation (92) is always positive and approaches zero when [phi] [right arrow]

0. The curve with a maximum plotted in Figure (6) according to Equation (92) displays the behaviour of the error expressed in percents, [[delta].sub.1] x 100%, as a function of the volume fraction.

Model 2: Mass-Dissipation Balance

A stated in Different Sets of the Outer Boundary Conditions subsection, Model 2 is based on the Spherical Cell versions of the Mass Balance, Dissipation Balance and Kuwabara conditions given by Equations (70), (75) and (76), respectively. Because the model under consideration includes the Kuwabara condition (76), the expression for the unknown constant A remains the same as for the above discussed model, i.e., the constant A is determined from Equation (87).

The hydraulic resistivity attributed to Model 2, [R.sup.2], is obtained by combining Equations (9) and (44), as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (93)

Note that all the distributions under integrals in (93) differ from those represented under integrals in Equation (90) attributed to Model 1 due to the multiplicative constant B, only (see Equation (64)). Clearly, the final expression in (93) is independent of such a constant. Since Equations (90) and (93) contain the same functions under integrals, one can use Equation (90) for obtaining the hydraulic resistivity, [R.sub.2], and the HIF, [L.sub.2]. To this end, let us rewrite Equation (90) as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (94)

Realizing that [R.sub.2]/[R.sub.1] = [L.sub.2]/[L.sub.1], the HIF attributed to Model 2 is determined by making use of Equations (88), (93) and (94), as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (95)

The final expression in (95) is the same as that obtained in the original paper of Kuwabara (1959), Equation (14).

According to Different Sets of the Outer Boundary Conditions subsection, using Model 2 leads to the violation of the momentum balance condition given by Equation (74). The error due to the violation of Equation (74) can be easily evaluated as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (96)

where the first and second terms in the numerator give the pressure differences determined from the dissipation and momentum balance conditions, Equations (44) and (35), respectively. The second equality in (92) becomes obvious while comparing the first equalities in (90) and (92) and taking into account that Equations (90) and (92) contain the same distributions under respective integrals.

Model 3: Momentum-Dissipation Balance

Model 3 deals with outer boundary conditions given by Equations (74) to (76). Similarly to the above considered two models, Model 3 includes Kuwabara boundary condition (76). Accordingly, the constant A for Model 3 is given by Equation (87).

Combining Equations (9) and (44), the hydraulic resistivity attributed to Model 3 is represented as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (97)

To obtain the HIF, [L.sub.3], we will use the method similar to that employed above for obtaining [L.sub.2]. Accordingly, Equation (90) is rewritten, as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (98)

Using Equations (88), (89), (92), (97) and (98) and realizing that [R.sub.3]/[R.sub.1] = [L.sub.3]/[L.sub.1], one obtains

[L.sub.3] = [L.sub.1] (1 + [[delta].sub.1]) = 50/(50 + 51[[phi].sup.1/3] + 3[[phi].sup.2/3] - 14[phi] + 20[[phi].sup.5/3] + 24[phi] + 12[[phi].sup.7/3] + 4[[phi].sup.8/3]) [(1 - [[phi].sup.1/3]).sup.3] (99)

Thus, Equation (99) gives the HIF, for the Momentum-Dissipation Balance Model.

Among the above considered three models employing the Kuwabara (1959) condition, Model 3 yields the highest value of the HIF. One can observe that the above obtained values of the HIF satisfy the equality

[L.sub.1] = [square root of [L.sub.2][L.sub.3] (100)

Thus, the HIF predicted for Model 1 is the mean geometrical value of those obtained using Models 2 and 3.

While using Model 3, the mass balance condition is violated. The corresponding error can be evaluated as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (101)

Note that the error, [[delta].sub.3], is negative, i.e., the Momentum-Dissipation Balance Model under predicts the external flow velocity, [u.sub.[infinity]].

Model 4: The Cunningham Version of the Mass-Momentum-Dissipation Balance

As stated in Different Sets of the Outer Boundary Conditions subsection, combining conditions (70), (74) and (75), which provide the mass, momentum and dissipation balances, leads to a quadratic equation with respect to the unknown constant A. According to Different Sets of the Outer Boundary Conditions subsection, one of the roots of this quadratic equation can be found using condition (82). Such a reformulated problem is equivalent to the Cunningham (1910) model.

Combining Equations (1), (64) and (82), the unknown constant, A, in Equation (64) is determined in the form

A = - 3[phi] ([[phi].sup.1/3] + 1/9[[phi].sup.4/3] + 9[phi] + 4[[phi].sup.2/3] + 4[[phi].sup.1/3] + 4 (102)

The hydraulic resistivity is obtained by combining Equation (9) with the mass and momentum balance conditions given by Equations (70) and (74), respectively:

[R.sub.4] = b[PI](b) + 2[eta][[d.sup.2]/[dr.sup.2](b) - 2/[b.sup.2]Y(b)/ 2Y(b) (103)

Substituting into the above equation the functions [PI](r) and Y(r) given by Equations (64) and (65) and using Equation (102) one obtains the final expression for [R.sub.4]. Normalizing R4 according to Equations (10) and (11), the HIF attributed to the Cunningham model takes the form

[L.sub.4] = 4 (1 - [[phi].sup.5/3]/4 + 7[[phi].sup.1/3] + 4 [[phi].sup.2/3]) [(1 = [[phi].sup.1/3]).sup.4] (104)

Equation (104) is equivalent to the Cunningham (1910) result given by Equation (12).

Clearly, using the Cunningham model leads to the violation of the Kuwabara condition. For evaluating the corresponding error, recall that, in the Generalized Kuwabara Condition subsection, the Kuwabara condition was derived in a generalized form given by Equation (57). The derivation was based on the requirement that two obtained expressions for the external pressure difference, Equations (40) and (56), should simultaneously be valid. Accordingly, the Kuwabara condition given by Equation (57), or its Spherical Cell equivalent Equation (76), imposes zero value of the vorticity term in Equation (40). Due to the presence of such a term, Equation (40) differs from Equation (56). While using the Mass-Momentum-Dissipation Balance Models, such a term takes a non-zero value, and thus Equations (40) and (56) contradict each other. Hence, the relative contribution of the vorticity term into the mean pressure gradient given by Equation (40), can be employed as a parameter characterizing the violation of the Kuwabara condition. Consequently, using Equations (40) and (57), the parameter characterizing the error due to the violation of the Kuwabara condition is defined as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (105)

Combining Equations (1), (64), (65), and (105), after some transformations, the parameter [[delta].sub.4] is represented as a function of the volume fraction, [phi]:

[[delta].sub.4] ([phi]) = - 1/3 (1 = [[phi].sup.1/3]) 3 + 6[[phi].sup.1/3] + 4[[phi].sup.2/3] + 2[phi]/1 + [[phi].sup.1/3] + 2[[phi].sup.2/3] + [phi] + [[phi].sup.4/3] (106)

Thus, Equation (106) describes the error, [[delta].sub.4] ([phi]) due to the violation of the Kuwabara condition while using the Cunningham (1910) model. In Figure (7) we plotted the curve which displays the behaviour of the function [[delta].sub.4] ([phi]) according to Equation (106).

Remarkably, the function [[delta].sub.4] ([phi]) is negative and does not approach to zero when [phi] [right arrow] 0. Moreover, at [phi] [right arrow] 0, the error has the highest magnitude which reaches 100%. Paradoxically, the strongest disagreement between the results of using Equations (40) and (56) occurs for the case of [phi] [right arrow] 0 when the Cunningham model yields an exact prediction, [L.sub.4] = 1. In Comparison of Different Hydrodynamic Cell Models subsection, such a surprising behaviour will be discussed in details.

Model 5: The Happel Version of the Mass-Momentum-Dissipation Balance

Now, we consider the second root of the quadratic equation which is obtained by combining the mass, momentum and dissipation balance conditions given by Equations (70), (74) and (75). Recall that, in the Model 4: The Cunningham Version of the Mass-Momentum-Dissipation Balance subsection, the first of the roots was obtained from Cunningham condition (82). As shown in the Different Sets of the Outer Boundary Conditions subsection, the second root can be determined from the Happel boundary condition given by Equation (86). Consequently, combining Equations (64) and (86) one obtains the integration constant, A, in the form:

A = - [[phi].sup.5/3]/ 3[[phi].sup.5/3] + 2 (107)

For obtaining the hydraulic resistivity, we combine Equations (64), (65), (103) and (107). Making use of the normalization given by Equations (10) and (11) we arrive at the Happel result expressed in terms of the HIF, as:

[L.sub.5] = 2/3 x 3 + 2[[phi].sup.5/3]/(2 + [[phi].sup.2/3] + 2[[phi].sup.2/3]) (1 + [[phi].sup.1/3]) (1 + [[phi].sup.1/3]) [(1 - [[phi].sup.1/3]).sup.3] (108)

Equation (108) is equivalent to the Happel result given by Equation (13).

As stated in the Different Sets of the Outer Boundary Conditions subsection, using the Happel model leads to the violation of Kuwabara condition (76). We will evaluate the corresponding error, [[delta].sub.5], in a similar manner as we evaluated the error [[delta].sub.4] in the Model 4: The Cunningham Version of the Mass-Momentum-Dissipation Balance subsection. Accordingly, we will consider the contribution of the vorticity term into the pressures difference given by Equation (74).

[[delta].sub.5] = -2[eta] [[d.sup.2]Y/[dr.sup.2](b)] - 2/[b.sup.2]Y(b)/-b [PI](b)(109)

Equation (109) differs from Equation (105), which was employed for obtaining [[delta].sub.4], due to the presence of the purely pressure term in the denominator instead of the overall pressure term represented in Equation (105). Such a modification is required to follow the approach used in the previous sections where, while determining the corresponding errors, we placed in the denominator the smaller of two quantities whose difference gives the absolute error.

Combining Equations (64), (65), (107) and (109) one obtains:

[[delta].sub.5] = 2 [(1 - [[phi].sup.1/3]).sup.2] 3 + 6[[phi].sup.1/3] + 4[[phi].sup.2/3] + 2[phi]/3 + 10[[phi].sup.2/3] + 2[[phi].sup.5/3] (110)

Thus, Equation (110) gives the error, due to the violation the Kuwabara condition within the frameworks of the Happel model. The obtained error is positive, i.e., the overall pressure is more than that given by the purely pressure term on the left-hand side of (74). The curve displaying the behaviour of the error [[delta].sub.5]([phi]) expressed in percents is given in Figure 7. In the present case, similarly to the error [[delta].sub.4]([phi]), obtained in the Model 4: The Cunningham Version of the Mass-Momentum-Dissipation Balance subsection, the largest error, 200%, is observed for the case [phi] [right arrow] 0 when the Happel models gives exact prediction of the HIF, [L.sub.5] = 1. Detailed discussion of this behaviour will be given in the Comparison of Different Hydrodynamic Cell Models subsection.

[FIGURE 7 OMITTED]

DISCUSSION

In the previous two Section, Spherical Cell Approximation: Reformulation of General Problem, and Analysis of Predictions and Errors from Different Models, we considered five models giving noticeably different predictions for the HIF defined in the Hydrodynamic Interaction Factor subsection. Each of the models employs three from four conditions which were derived in the third Section, Representative Cell of the General Type, to express the mass, momentum and dissipation balance within the system and the continuity of the local pressure within the fluid. Accordingly, one of four conditions is unavoidably violated within the framework of any of the five models (Analysis of Predictions and Errors from Different Models, fifth section).

The impossibility to satisfy all the conditions, which would be valid at the outer surface of the "regular" cell, witnesses certain imperfectness of the hydrodynamic version of the Spherical Cell Approximation. Due to such imperfectness, the hydrodynamic version of the Spherical Cell Approximation differs from its electrodynamic version for which all the necessary conditions are satisfied, simultaneously (Shilov et al., 1981).

Below, using the electrodynamic version of the Spherical Cell Approximation, we briefly describe the general scheme of obtaining the effective electric resistivity of diaphragm formed by perfect dielectric particles inside a conducting medium. Such a problem is the closest electrodynamic analogue of the hydrodynamic problem considered in the present paper.

Hydrodynamic and Electrodynamic Versions of the Spherical Cell Approximation

Let us consider how the Spherical Cell approximation can be employed for determining the effective electric resistivity, [R.sup.*.sub.el], of the diaphragm which was considered in previous section. Assume that the solid particles forming the diaphragm are perfect dielectrics, and the fluid has a finite specific electric resistivity, [R.sup.[infinity].sub.el]. Using Figure 2, the effective electric resistivity can be defined similarly to the hydraulic resistivity, R, given by Equation (9):

[R.sup.*.sub.el] = [DELTA][PHI]/[h.sub.1] [J.sub.[infinity]] (111)

where [DELTA][PHI] is the electric potential difference between the external sides of the diaphragm; and [J.sub.[infinity]] is the magnitude of the external electric current density.

To determine the effective electric resistivity one should obtain the distributions of the electric potential and current density inside the conducting fluid. The distributions are determined by solving a boundary value problem which is formulated on the basis of the principle of electric charge conservation. The latter principle leads to two relationships

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] in the fluid bulk (112)

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] at the particle surfaces (113)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the local density of the electric current.

By using the relationship between the electric current density, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and electric potential gradient (the differential version of Ohm's law),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (114)

Equation (112) is transformed into the Laplace governing equations describing the electric potential distribution inside the fluid

[[nabla].sup.2][PHI] = 0 (115)

Consequently, Equation (113) is transformed into the boundary condition,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], at the particle surfaces (116)

To close the problem formulation, one should impose boundary conditions at the outer boundary of the diaphragm. By introducing planes A and B shown in Figure 4, the required outer boundary conditions are represented as:

[PHI] = [DELTA][PHI] at the plane A (117)

[PHI] = 0 at the plane B (118)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] at the planes A or B (119)

Boundary conditions (117) to (119) are formulated on the basis of the same arguments and assumptions as those used in the General Formulation section for obtaining boundary conditions (20), (21) and (23). In particular, it is assumed that the diaphragm is constituted by representative cells whose maximum dimension b and the length scale parameters of the system shown in Figure 4 satisfy inequality (16).

Solving the boundary value problem given by Equations (115) and (118), and combining the obtained solution with Equations (111) and (119), one obtains the effective electric resistivity for the diaphragm which is constituted by representative cells.

Considering representative cell of the general type, one can derive two conditions that are complete analogues of the relevant conditions obtained in the third section, Representative Cell of the General Type. The first of such conditions expresses the electric charge balance in the system. Using Equation (112), which is an analogue of hydrodynamic continuity Equation (18), and conducting the derivations similar to those presented in the Mass Balance subsection for describing mass balance condition (31), we arrive at the electric charge balance condition:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (120)

Now, we obtain the second of the above mentioned conditions. Similarly to the pressure, the electric potential is a continuous function of the coordinates within the conductive medium. Consequently, one can repeat all the steps which were presented for obtaining Equation (56). Such an approach leads to an equation which is complete analogue of Equation (56):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (121)

Additionally to Equations (120) and (121), one can suggest a relationship expressing the dissipation balance in the system. Similarly to the condition derived in the Dissipation Balance subsection, we equate the local dissipation rates averaged over the diaphragm and representative cell volumes. The corresponding equality takes the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (122)

The left-hand side of Equation (122) yields the power developed by external source per unit diaphragm volume and time. Expression under integral on the right-hand side of Equation (122) is the local Joule dissipation rate which, using continuity Equation (112), can be rewritten as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (123)

Consequently, making use of Equations (113), (123) and the Gauss theorem, dissipation balance condition (122) can be represented with the help of a surface integral, as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (124)

While using the Spherical Cell Approximation, the angular symmetry of the electric potential distribution is assumed to be the same as in the case of an individual particle placed in a uniform electric field. The latter defines the form of the electric potential which is similar to that given by Equation (60) for the local pressure. Accordingly, using the spherical coordinate system shown in Figure 5 one obtains

[PHI] = [gamma](r) cos ([theta]) (125)

For uncertain outer boundary conditions, the function [gamma](r) is given by the solution of Laplace Equation (115) subject to boundary conditions (116):

[gamma](r) = C (r + [a.sup.3]/2[r.sup.2]) (126)

where the integration constant C is determined using the cell outer boundary conditions.

Similarly to the analysis given in the Outer Boundary Conditions subsection, we will obtain the Spherical Cell versions of conditions (120) to (122) derived for a representative cell of the general type. Consequently, making use of Equations (114) and (125), boundary condition (120) takes the form

1/[R.sup.[[infinity].sub.el] d[gamma]/dr (b) = - [J.sub.[infinity]] (127)

As well, using substitution (125), boundary condition (121) is specified as:

[DELTA][PHI]/[h.sub.1] = - [gamma](b)/b (128)

Two conditions (127) and (128) are sufficient to determine the electric resistivity, [R.sup.*.sub.el]. Accordingly, combining Equations (111) and (126) to (128) yields

[R.sup.*.sub.el] = [R.sup.[infinity]]/sub.el] [PHI](b)/b d[PHI]/dr (b) = 1/2 [R.sup.[infinity]]/sub.el] 2 + [phi]/1 - [phi] (129)

The latter expression in (129) is equivalent to the classical Maxwell (1881) expression specified for the case of perfect dielectric particles inside a conducting medium.

Above, we derived the third condition given by Equation (124) that was not employed for obtaining Equation (129). It might appear that we face the situation similar to that considered in the previous two sections, Spherical Cell Approximation: Reformulation of General Problem, and Analysis of Predictions and Errors from Different Models, where all the derived outer boundary conditions could not be satisfied, simultaneously. However, for the present case of the electrodynamic Spherical Cell Approximation, the omitted dissipation balance condition given by Equation (124) turns out to be valid when two other conditions (127) and (128) are satisfied. To demonstrate this, we transform the right-hand side of Equation (124) using Equations (125), (127) and (128):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (130)

Thus, starting with the left-hand side of dissipation balance Equation (124) we arrived at its right-hand side given by the latter expression in (130).

The above analysis shows that the contradiction, which was observed in the previous two sections, Spherical Cell Approximation: Reformulation of General Problem, and Analysis of Predictions and Errors from Different Models, for the hydrodynamic version of the Spherical Cell Approximation, does not exist in the case of the electrodynamic version of the Spherical Cell Approximation: the result obtained using boundary conditions (128) and (127) does not contradict the "spare" condition given by Equation (122).

In the next section, we will demonstrate that, among the five hydrodynamic Spherical Cell models, better approximations are given by the models which result in smaller error due to the violation of the "spare" condition.

Comparison of Different Hydrodynamic Cell Models

Now, we will conduct a comparative analysis of the errors associated with using five different models considered in the fifth section, Analysis of Predictions and Errors from Different Models. In Figures 7a and b, for convenience of consideration, we present the curves displaying the behaviour of each of the five errors as a function of the volume fraction. Recall that each of the errors describes violation of one of the conditions, which should be valid at the outer boundary of the "regular" cell and is unavoidably omitted while using the Spherical Cell Approximation.

As discussed in the Cunningham, and Happel version of the Mass-Momentum-Dissipation Balance subsections, for infinitely diluted system, [phi] [right arrow] 0, the errors attributed to the models based on the Kuwabara conditions approach zero. Simultaneously, at [phi] [right arrow] 0, the errors due to the violation of the Kuwabara condition within the frameworks of the Cunningham and Happel models have the highest magnitudes [absolute value of [[delta].sub.4]] = 100% and [absolute value of [[delta].sub.4]] = 200%.

As we noted in the Cunningham, and Happel version of the Mass-Momentum-Dissipation Balance subsections, the largest errors associated with the Cunningham and Happel models are observed when both these models yield exact prediction of the HIF, which, at [phi] [right arrow] 0, approaches to unity for all the models. To understand such a surprising behaviour of the functions [[delta].sub.4] ([phi]) and [[delta].sub.5]([phi]), let us combine Equations (1), (64), (65) and (74) to represent the mean pressure gradient given by Equation (74) in the following form

[DELTA]P/[h.sub.1] = - 3[eta] B/[a.sup.4] (5A + 3) [phi] (131)

As well, we present separate expressions for the vorticity and purely pressure parts of the mean pressure gradients, which are represented on the left-hand side of Equation (74).

- 2/[b.sup.2] [eta] [[d.sup.2]Y/[dr.sup.2] (b) - 2/[b.sup.2] Y(b)] = [eta] B/[a.sup.4] [-20A - 2 (5A + 3) [phi]] (132)

- [PI](b)/b = [eta] B/[a.sup.4] [20A - (5A + 3) [phi]] (133)

Additionally, combining Equations (1), (64) and (70), we express the external flow velocity [u.sub.[infinity]], as

[u.sub.[infinitiy]] = - 2 B/[a.sup.2] [1/2 (A + 3) [phi] - 1/2 (5A + 3) [[phi].sup.1/3] + 1 + A/[[phi].sup.2/3]] (134)

While inspecting Equations (131) and (134) at [phi] [right arrow] 0, it becomes clear that, irrespectively of the specific form of the function A([phi]), the hydraulic resistivity given by Equation (9), R = [DELTA][rho]/[h.sub.1] [u.sub.[infinity]], approaches its asymptotic form given by Equation (11) for any set of the outer boundary conditions provided that A([phi]) = o([[phi].sup.2/3]). In particular, at low volume fractions, using Equations (131) and (134) one obtains a common limiting value of the hydraulic resistivity given by Equation (11) for A [equivalent to] 0 (zero hydrodynamic interaction), A([phi]) = O([phi]) (Kuwabara and Cunningham conditions, Equation (87) and (102)) and A([phi]) = O([[phi].sup.5/3]) (Happel condition, Equation (107)). Thus, within a wide range of the errors in determining the constant A, the limiting value of the HIF, remains unaffected by the errors.

As shown above, at [phi] [right arrow] 0, the mean pressure gradient, [DELTA]P/[h.sub.1], and the external flow velocity, [u.sub.[infinity]], become independent of the constant A. For the same limiting case, the ratio of the vorticity to purely pressure term in Equation (74) depends on A (see Equations (132) and (133)). Accordingly, for different outer boundary conditions, which define different values of A, the ratio of the vorticity to purely pressure term takes different values: 0 (for Kuwabara condition), -1/2 (for Cunningham condition) and 2 (for both the Happel condition and the case of zero hydrodynamic interaction, A = 0). Such differences in the behaviour occur due to the presence of the first terms in the quadratic brackets on the right-hand sides of (132) and (133), [+ or -] 20A. Note that these terms are cancelled out while obtaining the overall pressure difference given by Equation (131).

The above discussion gives an explanation why, at low volume fractions, the high errors, which are associated with violation of the Kuwabara condition, do not influence on the limiting value of the HIF L. Such an unaltered limiting value, L = 1, which is independent of the unknown integration constant A, gives the zeroth-order term in the series expansion of the HIF by powers of [[phi].sup.1/3]. Note that the next terms in such an expansion depend on the constant A, and take different forms for different models. To illustrate this, we consider the Taylor series expansion of L by powers of [[phi].sup.1/3] and retain three leading terms, only. Such expansions take a common general form for Equations (12), (13) and (15), which give the HIF attributed to the Cunningham, Happel and the Mass-Momentum Balance, models, respectively:

L = 1 + k[[phi].sup.1/3] + [k.sup.2] [[phi].sup.2/3] (135)

In Equation (135), the coefficient k depends on the constant A and thus on the imposed outer condition. Combining Equations (9) to (11), (131) and (134) the coefficient k can be represented as:

k = 3/2 - A/[phi] (136)

While neglecting by the hydrodynamic interaction between the particles, i.e., assuming that A = 0, one obtains from Equation (136) that k = 3/2. Exactly the same value, k = 3/2, corresponds to the Happel model for which, according to Equation (107), A([phi]) = O([[phi].sup.5/3]). Thus, within the frameworks of the Happel model, the asymptotic behaviour of the HIF, [L.sub.5], at low volume fractions is defined by the hydrodynamic field around a single particle placed in a uniform flow. For such a case, the dependency of [L.sub.5] on the volume fraction occurs due to the averaging of such an unaltered local velocity over the Spherical Cell volume.

For Model 1 (Mass-Momentum Balance), the Kuwabara condition leads to Equation (87) which yields A = -3[phi]/10 + 0([[phi].sup.2]). The Cunningham model leads to Equation (102) which can be represented as A = -3[phi]/4 + 0([[phi].sup.4/3]). Consequently, making use of Equation (136), one obtains k=9/5, for Model l, which is based on the Kuwabara condition, and k = 9/4, for the Cunningham model.

Note that the integration constants, A, obtained using the Kuwabara and Cunningham outer boundary conditions, differ from each other by factor of 2.5. However, due to the presence of the term 3/2 in Equation (136), the difference between the coefficients k following from the Kuwabara and Cunningham boundary conditions is not so substantial.

Let us now consider the behaviour of the errors originated from the violation of the omitted boundary conditions. As it becomes clear from Figure 7a, the result of the error comparison depends on the volume fractions. At sufficiently low volume fractions, [phi] < 0.1, the Cunningham and Happel models, 4 and 5, lead to higher errors, [absolute value of [[delta].sub.4]] 65% and [absolute value of [[delta].sub.4]] > 75%, than Models 1, 2 and 3 which use the Kuwabara condition, [absolute value of [[delta].sub.1]] = [absolute value of [[delta].sub.2]] = [absolute value of [[delta].sub.3]] < 16%. Thus, according to the comparison of the errors represented in Figure 7a, at sufficiently low volume fractions, the results based on Models 1, 2 and 3, which use the Kuwabara condition, should give a better approximation than Models 4 and 5.

When the volume fraction approaches its value corresponding to the close packing, [[phi].sub.*], the error of the Happel model becomes the smallest. In Figure 7b, which displays the range of high volume fractions, when [phi] [right arrow] [phi] = [pi]/6 (simple cubic lattice), [[delta].sub.5] [approximately equal to] 8.5% whereas [[delta].sub.1] [approximately equal to] 20% and [[delta].sub.4] [approximately equal to] 21%. For [phi] [right arrow] [[phi].sub.*] = [square root of 2[pi]/6 (face-centred cubic lattice), [[delta].sub.5] [approximately equal to] 2% whereas [[delta].sub.1] [approximately equal to] 12% and [[delta].sub.4] [approximately equal to] 10%.

For Models 1, 2 and 3, the errors produced due to the violation of the omitted condition have the same magnitudes for all the three models. Therefore, analysis of the errors, whose behaviour displayed in Figures 7a and b, does not give solid arguments for choosing the best version from Models 1, 2 or 3. The negative sign of the errors [[delta].sub.2] and [[delta].sub.3] indicates that the hydraulic resistivity is under predicted using Model 2 and over predicted using Model 3. Thus, at for low volume fractions, the curve displaying a correct dependency, L([phi]), is expected to belong to the field between the curves [L.sub.2]([phi]) and [L.sub.3]([phi]) shown in Figure 8. The curves plotted for Cunningham and Happel cases, run outside such field. However, it should be noted that the curve obtained from the Happel model nearly merges with the lower boundary of the filed.

According to Equation (100), the curve displaying the behaviour of the function L1([phi]) is a mean geometrical value of [L.sub.2]([phi]) and [L.sub.3]([phi]). Hence, the curve [L.sub.1]([phi]) belongs to the field bound by the curves [L.sub.2]([phi]) and [L.sub.3]([phi]) in Figure 8. It is interesting that, for sufficiently low volume fractions, the HIF obtained for the Mass-Momentum Balance Model, [L.sub.1]([phi]), is in an excellent agreement with the results obtained by analyzing flow through the simple cubic lattice formed by particles. To check this, let us compare the predictions given by Equation (89) or, equivalently, Equation (15), with the results obtained by Hasimoto (1959) for the a simple cubic lattice. In terms of the HIF, the Hasimoto result, can be represented in the following form

[L.sub.Hasimoto] = 1/ 1 - 1.7601 [[phi].sup.1/3] + [phi] - 1.5593 [[phi].sup.2] + O ([[phi].sup.8/3]) (137)

Remarkably, the right-hand sides of Equations (15) and (137) differ from each other due to the difference between the denominators that contain the same powers of [[phi].sup.1/3] but different coefficients at [[phi].sup.1/3] and [[phi].sup.2]. Note that, in Equations (15) and (137), the coefficients at [[phi].sup.1/3] are very close to each other. In Equation (15) this coefficient is (-9/5=-1.8) whereas in Hasimoto's Equation (137) it is (-1.7601).

The curves shown in Figure 9, demonstrate an excellent agreement between the HIFs given by Equations (15) and (137) for [phi] < 0.15. The curve in Figure 10, which displays the deviations of [L.sub.1] from Hasimoto result (137), has a maximum that does not exceed 4%. While Equation (15) is in a good agreement with the Hasimoto prediction, the Cunningham and Happel results noticeably deviate from the curve plotted according to Equation (137). When the volume fraction, [phi], becomes higher than 0.15, the magnitude of the discrepancy between the predictions given by Equations (15) and (137) rapidly increases with increasing the volume fraction. Such a behaviour could be expected since, for [phi] > 0.15, the Hasimoto expression result substantially over predicts the HIF.

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

According to the above observation, while increasing the volume fraction, one can expect a decrease in the accuracy of using Equation (15). As well, at sufficiently high volume fractions, the Happel result given by Equation (13) is expected to give the best description. In order to examine whether the expected behaviour does take place, let us compare the predictions given by Equations (12), (13) and (15) with an exact solution obtained by Sangani and Acrivos (1982) for the simple cubic cell.

The curves shown in Figure 11a, demonstrate that Equation (15) gives an excellent approximation of the exact solution for [phi] < 0.2. The estimated discrepancy between the Kuwabara and exact results is lower than 3%. When 0.2 < [phi] < 0.42, Equation (15) over predicts the HIF compared to the exact Sangani and Acrivos solution but the error is less than 10%. Finally, for 0.43 < [phi] < [phi] = [pi]/6, Happel's prediction (13) gives a better approximation of the exact result than Equation (15). However, this approximation is not very accurate: the estimated discrepancy between the Happel and exact results varies from -11.5% to 12.5% (Figure 11b). Within the whole range of [phi], the Cunningham result, substantially over predicts the HIF.

Thus, for [phi] < 0.42, Equation (15), which follows from Model 1 (Mass-Momentum Balance), gives the best approximation of the HIF, L, obtained for the simple cubic lattice by Sangani and Acrivos (1982). For higher volume fractions, 0.43 < [phi] < [pi]/6, the better approximation is given by Equation (13) which follows from Model 5, (Happel's model). For other type of lattices, one can expect that Equation (15) gives even better approximation. For example, according to Hasimoto (1959), while dealing with the body and side centred cubic lattices, the coefficient at [[pi].sup.1/3] in Equation (137) should be (-1.791) instead of (-1.7601), i.e., this coefficient becomes closer to (-1.8) represented in Equation (15).

[FIGURE 10 OMITTED]

The above analysis is valid for the case of the purely pressure driven flow through a diaphragm. However, all the conclusions retain their validity in the presence of a body force acting the liquid, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], provided that such a force can be expressed through a gradient of a harmonic function, [OMEGA], as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (138)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (139)

For addressing hydrodynamic flow in the presence of the body force having the structure given by Equations (138) and (139), the boundary value problem discussed in the second section, General Formulation, should be reformulated while replacing the pressure, p, by the effective pressure, p*, defined as

[p.sub.*] = p + [OMEGA] (140)

The reformulated problem contains the same governing equations and boundary conditions as the general problem given by Equations (17) to (23) and is intended for obtaining the distributions of the local velocity [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and effective pressure [p.sup.*]. In the Sedimentation, and Electrically Driven (Electroosmotic) Flow subsections, we will consider two examples for which the body force can be represented using Equation (138) and (139), and a situation when conditions (138) and (139) are not applicable.

[FIGURE 11 OMITTED]

Sedimentation

Let us assume that the plane parallel system shown in Figures 2 and 4 is a layer formed by particles that are sedimenting with velocity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], in the gravity field, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. We will analyze the sedimentation process by linking the reference system to the particles. Considering balance of the hydrodynamic and gravity forces exerted on such a sedimenting layer, the pressure difference between the external sides of the layer, [DELTA]p, takes the form

[DELTA]p = [p.sub.A] [p.sub.B] = - [[rho].sub.F] (1 - [phi]) + [[rho].sub.p] [phi]] [gh.sub.1] (141)

where [[rho].sub.F] and [[rho].sub.P] are the fluid and particle mass densities, respectively. While deriving Equation (141), we considered planes are A and B shown in Figure 4 and assumed that geometrical parameters of the system satisfy inequalities (16).

In the reference system linked to the particles, one observes a flow through the motionless particle ensemble. Such a flow is driven by both the external pressure difference given by Equation (141) and the gravity body force, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] which can be represented in the form given by Equation (138) where the function [OMEGA] is given as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (142)

In the reference system linked to the particles, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Consequently, replacing the pressure, p, by the effective pressure, [p.sub.*], and combining Equations (9), (11) and (140) to (142) the sedimentation velocity is expressed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (143)

Thus, for addressing sedimentation of concentrated disperse system, one should determine the HIF, L, represented in the denominator of the final expression in (143). The HIF is obtained by analyzing the purely pressure driven flow discussed in the earlier sections. Accordingly, all the above conclusions regarding the choice of a better version of the Spherical Cell Approximation are valid for addressing the sedimentation process.

It should be stressed that Equation (143) describes sedimentation of the particles moving with equal velocities and thus having unchanged positions with reference to each other. According to the discussion presented in Cell Model Approach and Spherical Cell Approximation subsection, such a model has a restricted applicability.

Electrically Driven (Electroosmotic) Flow

Now, we consider electroosmotic flow which is generated when a diaphragm is placed in the external electric field, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. For describing the electroosmotic flow, we will modify the formulation of the general boundary value problem discussed in the second section, General Formulation, by introducing the body force, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], on the right-hand side of Stokes Equation (17):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (144)

where [PHI] is the perturbation of the electric potential due to the presence of the external electric field, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the density of the electric space charge, which exists in thermodynamic equilibrium state, due to the interfacial charge separation. An opposite charge is accumulated at the particle surface due to the different physico-chemical mechanisms (Masliyah and Bhattacharjee, 2006). These two types of charges having opposite signs form the electrical Double Layers surrounding each of the constituent particles.

We will consider purely electroosmotic flow by assuming that the pressure difference applied to the diaphragm is zero. Consequently, boundary conditions (20) and (21) are rewritten as:

p = 0 at the plane A (145)

p = 0 at the plane B (146)

Analysis of the electroosmotic flow using the Spherical Cell Approximation is a complex problem considered by many authors (Levine and Neale, 1974; Shilov et al., 1981, 2004; Zharkikh and Shilov, 1981; Ohshima, 1997a, b, 1998, 2000a, b; Dukhin et al., 1999a, b, c, 2000; Ding and Keh, 2001; Carrique et al., 2001a, b, 2002; Lee et al., 2001, 2002; Dukhin and Goetz, 2002; Masliyah and Bhattacharjee, 2006; etc.). We confine our analysis to consideration of two indicative limiting cases: (a) strong overlap of the Double Layers and (b) vanishingly thin Double Layers (the Smoluchowski limiting case).

Strong Overlap of the Double Layers

In the thermodynamic equilibrium state, for sufficiently short inter-particle distances, the thermal motion of ions is sufficient to provide a uniform concentration distribution for each of the individual ions within the porous space in spite of the electrostatic interactions between the ions and the surface charges. The quantitative criterion for such a uniformity is formulated while comparing the half of a mean distance between the particle surfaces, that can be evaluated as (b - a), with the Debye radius, [r.sub.D], which gives the Double Layer thickness. Consequently, the uniform distribution of ions within the liquid phase is observed when:

b - a/[r.sub.D] [much less than] 1 (147)

Inequality (147) gives a condition of the strong Double Layer overlap.

For a uniform distribution of the ions within the liquid phase, both the local electric resistivity, [R.sub.el], and charge density, [[rho].sub.el], do not depend on coordinates. Note that, for the general case, [R.sub.el] [not equal to] [R.sup.[infinity]].sub.el], where [R.sup.[infinity]].sub.el] was introduced in the Hydrodynamic and Electrodynamic Versions of the Spherical Cell Approximation subsection to describe the specific resistivity of the liquid sufficiently far from the diaphragm. Making use of Equations (121) and (114) (where [R.sup.[infinity]].sub.el] is replaced by [R.sub.el]), one can show that the electric force given by Equation (144) can be represented in the form given by Equations (138) and (139) where the function [OMEGA] takes the form

[OMEGA] = [[rho].sub.el] [PHI] (148)

As stated in the end of the Comparison of Different Hydrodynamic Cell Models subsection, by introducing the effective pressure according to Equation (140), the mathematical problem intended for obtaining electroosmotic velocity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], is reduced to the hydrodynamic problem formulated in the second section, General Formulation. Accordingly, combining Equations (9), (140), (145), (146) and (148), one obtains

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (149)

The final expression in (149) represents the electroosmotic velocity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], through the applied electric potential difference, [DELTA][PHI], and the hydraulic resistivity, R, which is obtained by solving the hydrodynamic boundary value problem formulated in the second section, General Formulation.

Taking into account the electroneutrality of the diaphragm, the local charge density, [[rho].sub.el], can be interrelated with the particle charge q, as:

V [NV.sub.p]) [[rho].sub.el] = - qN (150)

where N total number of particles, [V.sub.p] = 4[pi][a.sup.3]/3 is the particle volume. Combining Equations (10), (11), (149) and (150), we arrive at the equation obtained by Shilov et al. (2004):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (151)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the electric field strength averaged over the diaphragm volume.

Thus, for the case of strong overlap of the Double Layers, the electroosmotic velocity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], is given by Equation (151) where the hydrodynamic interactions manifest themselves through the HIF, L, represented in the denominator. Under conditions of strong Double Layer overlap, similarly to the sedimentation process, all the conclusion, which were made in the Comparison of Different Hydrodynamic Cell Models subsection regarding the better choice of the model, are valid for the electrically driven flow.

Vanishingly thin Double Layers (the Smoluchowski limiting case)

Now we consider the limiting case of a vanishingly thin space charge layer for which

[r.sub.D]/a [right arrow] 0 (152)

As shown by Smoluchowski (1903), under condition (152), the electroosmotic flow through diaphragm can be addressed on the basis of the boundary value problem formulated with the only difference: the tangent velocity at the particle surfaces is not zero, as it is imposed by boundary condition (19), and satisfies the Smoluchowski slipping condition:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] at the particle surfaces (153)

where the distribution of the electric potential inside the liquid, [PHI], is determined as a solution of Equation (115) subject to boundary conditions (116) and (119). In Equation (153), [mu] is the electroosmotic mobility. As for the normal velocity it takes zero value as it is defined by Equation (19):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (154)

Strictly speaking, Smoluchowski boundary condition (153) is set at the outer boundary of the vanishingly thin Double Layer. According to Smoluchowski (1903), in spite of the zero thickness of the space charge layer, the electric body force acting within this layer produces finite increase of the tangent velocity which takes zero value at the actual particle surface.

Thus, considering electroosmotic flow under the limiting case given by Equation (152), the hydrodynamic field is determined as a solution of the boundary value problem given by hydrodynamic Equations (17) and (18) subject to boundary conditions (153) and (154), at the particle surfaces, and (22), (145) and (146), at the planes A and B (Figure 4). The distribution of the electric potential, [PHI], which is represented in boundary condition (153), is obtained as a solution of another boundary value problem given by Equations (115), (116) and (119).

Solution for the above outlined boundary value problems yields the hydrodynamic and electric fields in the system. Consequently, using Equation (23), one obtains the electroosmotic velocity, as [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. As shown by Smoluchowski (1903), irrespectively of the system geometry, following the above described scheme one obtains the Smoluchowski expression electroosmotic velocity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], in the following form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (155)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (156)

is the external electric field strength. Thus, Equation (155) gives an exact expression for the Smoluchowski limiting case of the electroosmotic velocity. Below, we compare this known result with the predictions from different versions of the Spherical Cell Approach.

The distribution of the electric potential within the Spherical Cell, [PHI], is given by Equation (126). Accordingly, eliminating the unknown constant, C, with the help of Equation (127) and using Equation (156) one obtains

[PHI] = - [E.sub.[infinity]]/1 - [phi] (r + [a.sup.3]/2[r.sup.2]) cos ([theta]) (157)

Using Equations (61), (63) and (157), boundary condition (153) is rewritten as:

[dY.sub.eo]/dr (a) = 3/ 2 (1 - [phi]) [mu] [E.sub.[infinity]] a (158)

Another boundary condition at the particle surface is derived by combining (61), (62) and (154)

Y (a) = 0 (159)

The boundary condition imposing a given pressure difference is set by Equation (74). Accordingly, taking into account conditions of zero pressure difference, Equations (145) and (146), one obtains

[PI] (b) + 2/b [eta] [[dsup.2]Y/[dr.sup.2] (b) - 2/[b.sup.2] Y(b) (160)

Using Equations (60) to (65), the solution of hydrodynamic Equations (17) and (18) subject to boundary conditions (158) to (160) is expressed through function Y(r) which depends of one unknown constant C, as:

Y (r) = 3/2 (1 - [phi]) (5C + 3) [mu] [E.sub.[infinity]] [a.sup.2] [C [(r/a).sup.4] + [(r/a).sup.2] - a/r (C + 1)] (161)

Let us now determine the unknown constant C using the Kuwabara, Cunningham and Happel boundary conditions which, in terms of the function Y(r) are given by Equations (76), (82) and (86), respectively:

C = 0 Kuwabara (162)

C = - 3 [[psi].sup.5/3]/2 + 3[[phi].sup.5/3] Cunningham (163)

C = [[phi].sup.5/3]/1 - [[phi].sup.5/3] Happel (164)

Combining Equations (70) and (161) yields an expression for the electroosmotic velocity [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Substituting one of Equations (162) to (163) into the obtained expression we arrive at the final expressions for the electroosmotic velocity that corresponds to the Kuwabara, Happel and Cunningham boundary conditions, respectively:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] Kuwabara (165)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] Cunningham (166)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] Happel (167)

[FIGURE 12 OMITTED]

Thus, prediction from the Kuwabara model given by Equation (165) coincides with the exact Smoluchowski (1903) result, Equation (155), for which the electroosmotic velocity is independent of the volume fraction. The electroosmotic velocities obtained using both the Cunningham and Happel boundary conditions depend on the volume fraction and take the Smoluchowski value at [phi] [right arrow] 0. The curves plotted in Figure 12 according to Equations (165) to (167) demonstrate that using the Happel conditions leads to an over prediction, of the exact result that is less than 14%. The Cunningham model appreciably under predicts the exact result.

CONCLUSIONS

1. We derived four integral relationships, which are given by Equations (31), (40), (48) and (57), for determining the hydraulic resistivity of a diaphragm that can be considered as a system of representative cells formed by spherical particles. With the help of integrals over the outer boundary of a representative cell, the derived relationships interrelate the hydrodynamic field inside a representative cell with the external flow velocity, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] the applied pressure difference, [DELTA]p, and the consumed dissipation rate, [u.sub.[infinity]] [DELTA]p.

Three of the derived expressions, Equations (31), (40) and (48), describe the balances of, respectively, mass, momentum, and dissipation in a diaphragm constituted by identical, representative, cells. The fourth condition (the generalized Kuwabara condition) given by Equations (57) reflects the fact that the local pressure [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is a continuous function of the coordinates within the fluid. According to the derivation given, the generalized Kuwabara condition is valid for the systems where the momentum balance given by Equation (40) is satisfied.

2. Within the frameworks of the Spherical Cell Approximation, derived relationships (31), (40), (48) and (57) are transformed into four boundary conditions imposed at the Spherical Cell outer boundary, Equations (70) and (74) to (76). At the same time, only three outer boundary conditions are required for predicting the diaphragm hydraulic resistivity. For the case of an actual representative cell (elementary cell of a lattice, for example), all the four conditions given by Equations (31), (40), (48) and (57) are satisfied, simultaneously. Another situation takes place for the case the Spherical Cell Approximation: while using any group of three from four boundary conditions (31), (40), (48) and (57), the omitted condition becomes unavoidably violated.

3. It was shown that the mass, momentum and dissipation balance conditions given by Equations (70), (74) and (75) are satisfied, simultaneously, if, and only if, either the Cunningham or Happel outer boundary conditions are satisfied.

4. Choosing different groups of three outer boundary conditions, we formulated five models which differ from each other by the omitted, "spare" condition: (1) Mass-Momentum Balance Model; (2) Mass-Dissipation Balance Model; (3) Momentum-Dissipation Balance Model; and (4 and 5) the Cunningham, and Happel versions of the Mass-Momentum Dissipation.

5. On the basis of each of the five proposed models, we determined the Hydrodynamic Interaction Factor (HIF), i.e., the hydraulic resistivity normalized by its asymptotic value which corresponds to low volume fractions. The obtained results coincide with the predictions obtained earlier by others. For each of the models we evaluated the error due to the violation of the respective omitted, "spare" condition.

6. We compared the hydrodynamic version of the Spherical Cell Approximation with its electrodynamic version which leads to the classical Maxwell expression for the diaphragm effective electric resistivity. It was demonstrated that the contradictions observed for the hydrodynamic version do not exist in the case of the electrodynamic version. Although, for the electrodynamic version, the number of the obtained boundary conditions also exceeds the number defined by the order of the employed differential equations, all the obtained conditions turn out to be satisfied, simultaneously.

7. Comparative analysis of the errors due to the violation of the "spare" conditions, which are omitted within the frameworks of each of the five hydrodynamic models, shows that, at sufficiently low volume fractions, Models 1-3 that are based on the Kuwabara condition give a better approximation than the Cunningham and Happel models, 4 and 5. Comparison of the prediction from Model 1, Equation (15), with the results of Hasimoto (1959) and Sangani and Acrivos (1982) obtained for the simple cubic lattice confirm this conclusion. For [phi] < 0.2, the estimated discrepancy between the Kuwabara and exact Sangani and Acrivos result is smaller than 3% and when 0.2 < [phi] < 0.42, the discrepancy is smaller than 10%.

When the volume fraction approaches the value attributed to the close packing, the Happel model leads to the smallest error among the five models. This fact, which follows from the error analysis, correlates with the results of comparison with the Sangani and Acrivos (1982) solution. For 0.43 < [phi] < [[phi].sub.*] = [pi]/6, Happel's prediction (13) gives a better approximation of the exact result than Equation (15). However, this approximation is not very accurate: the estimated discrepancy between the Happel and exact results varies from -11.5% to 12.5% (Figure 11b).

8. The above conclusions regarding the best versions of the Spherical Cell Approximation are valid for purely pressure driven flow. As well, these conclusions are valid in the presence of a body force of arbitrary origin provided that such a body force can be represented as a gradient of a scalar function of coordinates. The latter was illustrated by two examples: the sedimentation and the electroosmotic flow for the case of strong overlap of the Double Layers.

9. We considered electroosmotic flow for the case of vanishingly thin Double Layer (Smoluchowski limit) as an example of the system for which the body force can not be represented as a gradient. For the Smoluchowski limit of the electroosmotic velocity, using the Kuwabara condition leads to the classical Smoluchowski exact result which is independent of the volume fraction. The predictions from the Cunningham and Happel models depend on volume fraction and coincide with the Smoluchowski result for [phi] [arrow right] 0.

10. The conducted analysis of two limiting cases for electroosmotic velocity shows that the Kuwabara boundary condition should be recommended for addressing electroosmosis. For the strong overlap of the Double Layers, similarly to the case of pressure driven flow, the Kuwabara condition gives the best approximation compared to other models when [phi] < 0.42. For vanishingly thin Double Layers, using the Kuwabara condition yields exact prediction for all volume fractions. All other situations are intermediate with respect to the above considered limiting cases.

ACKNOWLEDGEMENT

The authors gratefully acknowledge the financial support from the National Academy of Sciences of Ukraine within the frameworks of the program "Nano-Structured Systems, Nano-Materials and Nano-Technologies."
```NOMENCLATURE

a radius of a spherical particle
b radius of the Spherical Cell
d thickness of the gap between the
diaphragm and hypothetic planes
A and B (Figure 4)
[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] external electric field strength
[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] F force exerted on particle
[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] gravity acceleration

[h.sub.1], [h.sub.2],
[h.sub.3] linear dimension of diaphragm
[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] unit tensor
[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] unit vectors attributed to Cartesian
coordinate system
[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] unit vectors attributed to spherical
coordinate system
[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] electric current density

[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] external electric current density

L Hydrodynamic Interaction Factor (HIF)
N number of particles
[N.sub.cell] number of constituent representative
cells
[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] outward normal vector

p local pressure
[DELTA]p pressure difference between external
sides of a diaphragm
R effective hydraulic resistivity of
diaphragm
[R.sub.0] effective hydraulic resistivity for
low volume fractions
[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] outward normal vector
[R.sup.[infinity]].sub.el] specific electric resistivity of liquid
[R.sup.*.sub.el] effective electric resistivity of
diaphragm
S area
[U.sub.r](r), [U.sub.
of the velocity components within the
Spherical Cell
[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] local velocity of fluid

[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] external flow velocity

[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] electroosmotic velocity

[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] sedimentation velocity

[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] Smoluchowski limit of electroosmotic
velocity
[u.sub.r], [u.sub.[theta]] components of the local velocity vector
V diaphragm volume
[V.sub.cell] cell volume
[V.sub.p] particle volume
[x.sub.1], [x.sub.2],
[x.sub.3] Cartesian coordinates
Y(r) radial part of the streaming function

Greek Symbols

[[delta].sub.k] error due to the violation of the
omitted outer boundary condition
within the frameworks of the kth model
[eta] viscosity
[mu] electroosmotic mobility
of pressure within the Spherical Cell
[rho] electric charge density
[[rho].sub.F] liquid mass density
[[rho].sub.P] particle mass density
[MATHEMATICAL EXPRESSION NOT
REPRODUCIBLE IN ASCII] viscous stress tensor
[phi] volume fraction of particles
[PHI] electric potential
of electric potential within the
Spherical Cell
```

Manuscript received May 4, 2007; revised manuscript received June 20, 2007; accepted for publication June 20, 2007.

REFERENCES

Barnea, E. and J. Mizrahi, "A Generalized Approach to the Fluid Dynamics of Particulate Systems: Part 1. General Correlation for Fluidization and Sedimentation in Solid Multiparticle Systems," Chem. Eng. J. 5, 171-189 (1973).

Batchelor, G. K., "Sedimentation in a Dilute Dispersion of Spheres," J. Fluid Mech. 52, 245-268 (1972).

Batchelor, G. K., "Transport Properties of Two-Phase Materials with Random Structure," Ann. Rev. Fluid Mech. 6, 227-255 (1974).

Batchelor, G. K., "Sedimentation in a Dilute Polydisperse System of Interacting Spheres Part 1. General Theory," J. Fluid Mech. 119, 379-408 (1982).

Batchelor, G. K. and C. S. Wen, "Sedimentation in a Dilute Polydisperse System of Interacting Spheres. Part 2. Numerical Results," J. Fluid Mech., 124, 495-528 (1982).

Carrique, F., F. J. Arroyo and A. V. Delgado, "Electrokinetics of Concentrated Suspensions of Spherical Colloidal Particles: Effect of a Dynamic Stern Layer on Electrophoresis and DC Conductivity," J. Colloid Interface Sci. 243, 351-361 (2001a).

Carrique, F., F. J. Arroyo and A. V. Delgado, "Sedimentation Velocity and Potential in a Concentrated Colloidal Suspension. Effect of a Dynamic Stern Layer," Colloids Surf. A: Physicochem. Eng. Aspects 195 157-169 (2001b).

Carrique, F., F. J. Arroyo and A. V. Delgado, "Electrokinetics of Concentrated Suspensions of Spherical Colloidal Particles with Surface Conductance, Arbitrary Zeta Potential, and Double-Layer Thickness in Static Electric Fields," J .Colloid Interface Sci. 252, 126-137, (2002).

Chen, S. H, S. L. Lin, C. Kung and Y. M. Nian, "Concentration Effects on Sedimentation Velocity of Monodispersed Aerosol Particles," J. Chem. Eng. Japan 32, 635-644 (1999).

Cheng, H. and G. Papanicolaou, "Flow Past Periodic Arrays of Spheres at Low Reynolds Number," J. Fluid Mech. 335, 189-212 (1997).

Cunningham, E., "On the Steady State Fall of Spherical Particles through Fluid Medium," Proc. Roy. Soc. (London) A83, 357-365 (1910).

Dassios, G., M. Hadjinicolaou, F. A. Coutelieris and A. C. Payatakes, "Stokes Flow in Spheroidal Particle-in-Cell Models with Happel and Kuwabara Boundary Conditions," Int. J. Eng. Sci. 33, 1465-1490 (1995).

Datta, S. and S. Deo, "Stokes Flow with Slip and Kuwabara Boundary Conditions," Proc. Indian Acad. Sci. (Math. Sci.), 112, 463-475 (2002).

Davis, R. H. and A. Acrivos, "Sedimentation of Noncolloidal Particles at Low Reynolds Numbers," Ann. Rev. Fluid Mech. 17, 91-118 (1985).

Ding, J. M. and H. J. Keh, "The Electrophoretic Mobility and Electric Conductivity of a Concentrated Suspension of Colloidal Spheres with Arbitrary Double-Layer Thickness," J. Colloid Interface Sci. 236, 180-193 (2001).

Dukhin A. S., H. Ohshima, V. N. Shilov and P. J. Goetz, "Electroacoustics for Concentrated Dispersions," Langmuir 15, 3445-3451 (1999a).

Dukhin A. S., V. N. Shilov and Y. B. Borkovskaya, "Dynamic Electrophoretic Mobility in Concentrated Dispersed Systems. Cell Model," Langmuir 15, 3452-3457 (1999b).

Dukhin A. S., V. N. Shilov, H. J. Ohshima and P. J. Goetz, "Electroacoustic Phenomena in Concentrated Dispersions: New Theory and CVI Experiment," Langmuir 15, 6692-6706 (1999c).

Dukhin, A. S. and P. J. Goetz, "Ultrasound for Characterizing Colloids," in "Studies in Interfacial Science," D. Mobius and R. Miller, Eds., Elsevier (2002), p. 372

Epstein, N. and J. H. Masliyah, "Creeping Flow through Clusters of Spheroids and Elliptical Cylinders," Loughborough Univ. Technol. Chem. Eng. J. 3, 169-175 (1972).

Filippov, A. N., S. I. Vasin and V. M. Starov, "Mathematical Modeling of the Hydrodynamic Permeability of a Membrane Built Up from Porous Particles with a Permeable Shell," Colloids Surf. A 282-283, 272-278 (2006).

Frazer, R. A., "On the Motion of Circular Cylinders in a Viscous Fluid," Philosophical Transactions of the Royal Society of London, A 225, 93-130 (1926).

Gal-Or, B., "On Motion of Bubbles and Drops," Can. J. Chem. Eng. 48, 526-530, (1970).

Garside, J. and M. R. A1-Dibouni, "Velocity-Voidage Relationship for Fluidization and Sedimentation in Solid-Liquid Systems," Ind. Eng. Chem. Process Des. Dev. 16, 206-214 (1977).

Happel, J., "Viscosity of Suspensions of Uniform Spheres," J. Appl. Phys 28, 1288-1292 (1957).

Happel, J., "Viscous Flow in Multiparticle Systems: Slow Motion of Fluids Relative to Beds of Spherical Particles," AICHE J. 4, 197-201 (1958).

Happel, J. and H. Brenner, "Low Reynolds Number Hydrodynamics, with Special Applications to Particulate Media," Prentice-Hall: Englewood Cliffs, NJ (1965).

Hasimoto, H., "On the Fundamental Solutions of the Stokes Equations and Their Application to Viscous Flow Past a Cubic Array of Spheres," J. Fluid Mech. 5, 317-328 (1959).

Kvashnin, A. G., "Cell Model of Suspension of Spherical Particles," Fluid Dynamics 14, 598-602 (1979).

Kuwabara, S., "The Forces Experienced by Randomly Distributed Parallel Circular Cylinders or Spheres in a Viscous Flow at Small Reynolds Numbers," J. Phys. Soc., Japan 14, 527-532 (1959).

Leclair, B. P. and A. E. Hamielec, "Viscous Flow through Particle Assemblages at Intermediate Reynolds Numbers--A Cell Model for Transport in Bubble Swarms," Can. J. Chem. Eng. 49, 713-720 (1971).

Lee, K. W., L. D. Reed and J. A. Gieseke, "Pressure Drop across Packed Beds in the Low Knudsen Number Regime," J. Aerosol Sci. 9, 557-565 (1978).

Lee, E., F. Y. Yen and J. P. Hsu, "Dynamic Electrophoretic Mobility of Concentrated Spherical Dispersions," J. Phys. Chem. B 105, 7239-7245 (2001).

Lee, E., C. H. Fu and J. P. Hsu, "Dynamic Electrophoretic Mobility of a Concentrated Dispersion of Particles with a Charge-Regulated Surface at Arbitrary Potential," J. Colloid Interface Sci. 250, 327-336 (2002).

Levine, S. and G. H. Neale, "The Prediction of Electrokinetic Phenomena within Multiparticle Systems. I. Electrophoresis and Electroosmosis," J. Colloid Interface Sci. 47, 520-529 (1974).

Levine, S. and G. H. Neale and N. J. Epstein, "The Prediction of Electrokinetic Phenomena within Multiparticle Systems: II. Sedimentation Potential," J. Colloid Interface Sci. 57, 424-437 (1976).

Lynch, E. D. and E. Herbolzheimer, "Formation of Microscale Structure in Sedimenting Suspensions," Bull. Am. Phys. Soc. 28, 1365 (1983).

Masliyah, J. H., "Viscous Flow across Banks of Circular and Elliptic Cylinders: Momentum and Heat Transfer," Can. J. Chem. Eng. 51, 550-555 (1973).

Masliyah, J. H. and T. G. M. van de Ven, "A Two Dimensional Model of the Flow of Ordered Suspensions of Rods," Int. J. Multiphase Flow 12, 791-806 (1986).

Masliyah, J. H. and S. Bhattacharjee, "Electrokinetic and Colloid Transport Phenomena," John Wiley and Sons, Inc., Hoboken, NJ (2006) p. 707.

Maxey, M. R. and B. K. Patel, "Localized Force Representations for Particles Sedimenting in Stokes Flow," Int. J. Multiphase Flow 27, 1603-1626 (2001).

Maxwell J. C., "A Treatise of Electricity and Magnetism," 1, Clarendon Press, Oxford, 2nd ed., 435 (1881).

Mehta, G. D. and T. F. Morse, "Flow through Charged Membranes," J. Chem. Phys. 63, 1878-1889 (1975).

Neale, G. H. and W. K. Nader, "Prediction of Transport Processes within Porous Media: Creeping Flow Relative to a Fixed Swam of Spherical Particles," AICJE J. 20, 530-538 (1974).

Neale, G. H. and J. H. Masliyah, "Flow Perpendicular to Mats of Randomly Arranged Cylindrical Fibers (Importance of Cell Models)," AICHE J. 21, 805-807 (1975).

Ohshima, H., "Electrophoretic Mobility of Spherical Colloidal Particles in Concentrated Suspensions," J. Colloid Interface Sci. 188, 481-485 (1997a).

Ohshima, H., "Dynamic Electrophoretic Mobility of Spherical Colloidal Particles in Concentrated Suspensions," J. Colloid Interface Sci. 195, 137-148 (1997b).

Ohshima, H., "Sedimentation Potential of a Concentrated Suspension of Spherical Colloidal Particles," J. Colloid Interface Sci. 208, 295-301 (1998).

Ohshima, H., "Sedimentation Potential and Velocity in a Concentrated Suspension of Soft Colloidal Particles," J. Colloid Interface Sci. 229, 140-147 (2000a).

Ohshima, H., "Cell Model Calculation for Electrokinetic Phenomena in Concentrated Suspensions: An Onsager Relation between Sedimentation Potential and Electrophoretic Mobility," Adv. Colloid Interface Sci. 88, 1-18 (2000b).

Richardson, J. F. and W. N. Zaki, "Sedimentation and Fluidization: Part 1," Trans. Inst. Chem. Eng. 32, 35-53 (1954).

Ruiz-Reina, E., F. Carrique, F. J. Rubio-Hernandez, A. I. Gomez-Merino and P. Garcia-Sanchez, "Electroviscous Effect of Moderately Concentrated Colloidal Suspensions," J. Phys. Chem. B 107, 9528-9534 (2003).

Rubio-Hernandez, F. J., F. Carrique and E. Ruiz-Reina, "The Primary Electroviscous Effect in Colloidal Suspensions," Adv. in Colloid and Interface Sci. 107, 51-60 (2004).

Sangani, A. S. and A. Acrivos, "Slow Flow through a Periodic Array of Spheres," Int. J. Multiphase Flow 8, 343-360, (1982).

Shilov, V. N., N. I. Zharkikh and Y. B. Borkovskaya, "Theory of Non-Equilibrium Electro-Surface Phenomena in Concentrated Disperse Systems: 1. Application of Non-Equilibrium Thermodynamics to Cell Model of Concentrated Dispersions," Colloid J. U.S.S.R. (English translation) 43, 434-438 (1981).

Shilov, V. N., Y. B. Borkovskaya and A. S. Dukhin, "Electroacoustic Theory for Concentrated Colloids with Overlapped DLs at Arbitrary [kappa]a: I. Application to Nanocolloids and Nonaqueous Colloids," J .Colloid Interface Sci. 277, 347-358 (2004).

Simha, R., "A Treatment of Viscosity in Concentrated Suspensions," J. Appl. Phys 23, 1020-1024 (1952).

Smith, T. N., "The Differential Sedimentation of Particles of Various Species," Trans. Inst. Chem. Eng. 44, T153-T177 (1966).

Smith, T. N., "The Sedimentation of Particles Having Dispersion of Sizes," Trans. Inst. Chem. Eng. 45, T311-313 (1967).

Smoluchowski, M., "Contribution a la theorie de l'endosmose electroque et de quelques phenomenes correlatifs," Bull. Int. Academie des Sciences de Cracovie 8, 183-200 (1903).

Sparrow, E. M. and A. L. Loeffler, "Longitudinal Laminar Flow between Cylinders Arranged in Regular Array," AICHE J. 5, 325- 330 (1959).

Uchida, S., "Slow Viscous Flow through a Mass of Particles," Rep. Inst. Sci. Tech. Univ. Tokyo 3, 97-100 (1949).

Yaron, J. and B. Gal-Or, "Relative Velocities and Pressure Drop in Clouds of Drops, Bubbles or Solid Particles," AICHE J. 17, 1064-1074 (1971).

Yu, C. P. and T. T. Soong, "A Random Cell Model for Pressure Drop in Fibrous Filters," J. Appl. Mech.-Trans. ASME 42, 301-304 (1975).

Zharkikh, N. I. and V. N. Shilov, "Theory of Collective Electrophoresis of Spherical Particles in the Henry Approximation," Colloid J. U.S.S.R. (English translation) 43, 865-870 (1981).

Zholkovskiy, E. K., O. B. Adeyinka and J. H. Masliyah, "Spherical Cell Approach for the Effective Viscosity of Suspensions," J. Phys. Chem. B. 110, 19726-19734 (2006).

Ziemer, W. and B. Foerester, "A New Cell Model Filtration Theory," J. Aerosol Sci. 27, Supplement 1, S605-S606 (1996).

Zick, A. A. and G. M. Homsy, "Stokes Flow through Periodic Arrays of Spheres," J. Fluid Mech. 115, 13-26 (1982).

Emiliy K. Zholkovskiy (1) *, Vladimir N. Shilov (1), Jacob H. Masliyah (2) and Mykola P. Bondarenko (1)

(1.) Institute of Bio-Colloid Chemistry of Ukrainian Academy of Sciences, Vernadskogo, 42, 03142, Kiev, Ukraine

(2.) Department of Chemical and Materials Engineering, University of Alberta, 536 Chemical and Materials Engineering Building, Edmonton, AB, Canada T6G 2G6