Numerical Study on the Charge Transport in a Space between Concentric Circular Cylinders

Abstract

Electrification is one of the key factors to be considered in the design of power transformers utilizing dielectric liquid as a coolant. Compared with enormous quantity of experimental and analytical studies on electrification, numerical simulations are very few. This paper describes essential elements of numerical solution methods for the charge transport equations in a space between concentric cylinders. It is found that maintaining the conservation property of the convective terms in the governing equations is of the uttermost importance for numerical accuracy, in particular at low reaction rates. Parametric study on the charge transport on the axial plane of the annular space with a predetermined velocity shows that when the convection effect is weak the solutions tend to a one-dimensional nature, where diffusion is simply balanced by conduction. As the convection effect is increased the contours of charge distribution approach the fluid streamlines. Thus, when the conduction effect is weak, charge distribution tends to be uniform and the role of the convection effect becomes insignificant. At an increased conduction effect, on the other hand, the fluid motion transports the charge within the electric double layers toward the top and bottom boundaries leading to an increased amount of total charge in the domain.

1. Introduction

When a dielectric liquid containing impurities is in contact with a solid surface, a certain physicochemical process occurs at the interface yielding free ions near the surface of the liquid. Usually negative ions are adsorbed to the solid surface and the positive ions are diffused away forming the electric double layer (to be referred to as EDL). Since the positive ions are mobile, they are convected by the fluid flow giving rise to the streaming current, which is called flow electrification. Problem occurs when they are accumulated in a certain location downstream resulting in locally high electric-field intensity which can cause electrical discharge, breakdown, and local failure of the device employing the liquid transport.

Electrification becomes one of the key factors to be considered in the design of electrical devices utilizing dielectric liquid (mineral or ester oil) as a coolant, such as power transformers. Demand for increased capacity from the users of power transformers tempts designers to increase the oil flow rate for increased cooling capacity, which, however, brings increased electrification and makes the device more susceptible to the electrical failure.

Studies on electrification and discharge with full-scale transformers were carried out by Higaki et al. [1, 2] and Tamura et al. [3]. They inserted numerous sensors to measure the charge distribution within transformers and measured the local leakage current through the solid surfaces. Higaki et al. [1, 2] demonstrated that the point of maximum electric field on the solid surface, calculated by solving the potential equation with the leakage currents being used as the boundary conditions, was consistent with the point of discharge actually observed from the experiment. Tamura et al. [3] presented a diagram in the parameter space where high flow rate was shown to lead to electrical discharge.

In order to perform more fundamental studies, researchers have considered simple experimental apparatus other than actual transformers, which is easier to build and easier to measure data with, such as electrical charge tendency (ECT). In addition, how to interpret the measured data in relation to the actual transformers is also an important issue in selecting suitable geometries for study. Most of the initial studies focused on the flow between parallel plates and circular pipes [4–7]. A spinning disk system was used by Kedzia and Willner [8] and Gibbings [9] in their experimental and theoretical studies on electrification. In order to attain a fully developed flow in a compact space Washabaugh [10] and Washabaugh and Zahn [11, 12] used a circular Couette system. Moreau and Touchard [13] conducted an experimental study to show that an impinging jet can result in a surface current with an order of magnitude larger than that with the parallel flow. The swing cylinder system, where the inner cylinder shows back-and-forth rotation, was also used in [14] to study ECT.

With the simple flow apparatus in hand, researchers can perform investigations on the effect of various factors on the electrification independently. There are many factors influencing electrification or ECT. They can be categorized into two kinds, fluid/flow properties and electrical properties. Included in the former are flow rate, geometrical features determined by the fluid path, and fluid viscosity, while in the latter electrical conductivity and permittivity are the key elements; the operating temperature and material degradation may influence many of these properties, such as viscosity, conductivity and permittivity. Properties of the pressboard such as chemical composition of the material and surface roughness may also influence the ECT. General understanding of the effect of various parameters on ECT and design aspect for avoiding discharge in transformers was given in [15, 16].

Touchard [17] and Touchard [18] included detailed kinetics of wall surface reaction in the formulation of the charge flux through planar and circular duct to attain a reasonable matching with experimental data. Moreau et al. [19] measured ECT for a flow passing through a filter made of degraded pressboard to show that a degraded surface enhances the charge accumulation and that small amount of additive (i.e., BTA) can reduce ECT. The effect of additive was further studied in [20, 21], and by using a flow-loop apparatus, Bourgeois et al. [21] addressed the mechanism of the enhancement of ECT in terms of the carboxyl group. Aksamit and Zmarzly [22] also studied the inhibition of flow electrification with the additive 𝐶₆₀. Developing advanced models for the wall reaction that can well fit the experimental data is one of the most important issues in the study of electrification. Cabaleiro et al. [23] and El-Adawy et al. [24] for instance demonstrated that the wall-reaction constant must be varied depending on the other parameters such as flow rate in order to fit the experimental data. Cabaleiro et al. [25] performed analysis on ECT within a shallow rectangular duct with a more complex model for the wall reaction. Paillat et al. [26] showed that inclusion of the effect of the fluid shear stress in the physicochemical process at the interface provides a much better agreement with experimental data at high laminar Reynolds numbers. El-Adawy et al. [27] conducted numerical simulation for the foundation of EDL without the convection effect and calculated ECT in the presence of fluid flow with source terms representing the ion dissociation and recombination in the bulk. Kobayashi et al. [28] paid attention to the competitive role of oil and pressboard in the electrification process. Okabe et al. [29] investigated the effect of the compounds in oil on ECT and showed that increase of ECT was mostly caused by the oxidation of sulfides. In [30], both the transient and steady-state data of the electrification experiment could be matched with the one-dimensional model for ion transport by using a hybrid boundary condition, where a constant flux as well as the flux proportional to the local ionic concentration was employed. Due to environmental problems, the electrical power industry considers using substitutive oils in transformers. Paillat et al. [31] investigated experimentally the electrification property of ester oil compared with conventional mineral oil, in particular in terms of ECT. They confirmed that charge accumulation with ester oil is one or two orders of magnitude larger than that with mineral oil.

Most studies on electrification even with simple geometry have been performed experimentally and/or analytically, and numerical studies, at least in a two-dimensional space, are very few. In [13], numerical simulation was performed for the impinging jet configuration by using the authors’ in-house code. In [27], two-dimensional numerical simulation was performed to investigate the transient development of EDL and surface current by using comprehensive models incorporating ion dissociation and recombination of molecules not only in liquid but also in solid. The authors also included in their simulation the wall reaction model representing the combination of cation in the solid and anion in the fluid. In the second stage in their simulation they considered a fully developed parabolic velocity profile to investigate the time-dependent flow electrification. Their numerical results are qualitatively in line with experimental results and knowledge.

The main purpose of the present study is to develop a two-dimensional numerical code and perform simulations for charge transport in a confined space under a various range of parameters. In particular, we select as the computational domain the annulus between concentric cylinders, following [10–12]. We are concerned with the axisymmetric secondary flow developed on the axial plane at supercritical Reynolds numbers. One of the most important issues to be addressed in this study is the effect of nondimensional parameters on the numerical solutions as well as their accuracy. Aside from the one- and two-dimensional in-house codes, we also use the commercial software COMSOL and exact and approximate analytical solutions of the one-dimensional problem for verification of the numerical solutions. The present in-depth analysis of the characteristics of the equations governing the charge transport in relation to the numerical solutions may play an important role in the development of more practical simulation codes and in interpretation of the numerical results obtained either by an in-house or a commercial code.

2. Mathematical Formulation

We consider transport of a space charge density distributed in an annulus space between two concentric circular cylinders of radii 𝑅₁ and 𝑅₂, respectively, caused by diffusion, convection, and electrical conduction. The governing equations for the problem can be written as follows:𝜕𝑞^∗𝜕𝑡^∗+∇_∗⋅𝐉^∗∇=0,_∗⋅𝜀𝐄^∗=𝑞^∗,𝜕𝐮^∗𝜕𝑡^∗+𝐮^∗⋅∇_∗𝐮^∗1=−𝜌∇_∗𝑝^∗+𝜈∇²_∗𝐮^∗+𝑞^∗𝐄^∗,∇_∗⋅𝐮^∗=0, (1a)(1b)(1c)(1d)

where 𝑡^∗ is the time, 𝑞^∗ the charge density, 𝐮^∗ the fluid velocity, 𝑝^∗ the pressure, 𝐄^∗ the electric field related to the electric potential 𝜙^∗ as 𝐄^∗=−∇_∗𝜙^∗, ∇_∗ the gradient operator, and 𝐉^∗=−𝐷∇_∗𝑞^∗+𝑞^∗𝐮^∗+𝜎𝐄^∗ (2) the current density (or charge flux). Further, 𝜌 is the fluid density, 𝜈 the kinematic viscosity of the fluid, 𝐷 the diffusivity of the species (i.e., charge carriers), 𝜎 the electrical conductivity, and 𝜀 the electrical permittivity, all of which are assumed to be constant in this study.

It is assumed that metals or pressboards in contact with dielectric liquid create charges (or they may be adsorbed) by certain chemical reaction, and we simply employ the model used in [10–12] dictating that the charge flux from the electrode surface to the liquid is proportional to the local space charge density, reading𝐽^∗_𝑟=−𝑘^∗₁𝑞^∗−𝑞^∗_1𝑤at𝑟^∗=𝑅₁,𝐽^∗_𝑟=𝑘^∗₂𝑞^∗−𝑞^∗_2𝑤at𝑟^∗=𝑅₂, (3a)(3b)

where 𝐽^∗_𝑟 denotes the radial component of 𝐉^∗, 𝑞^∗_𝑖𝑤 the wall charge density at 𝑟^∗=𝑅_𝑖 (𝑖=1 for the inner and 𝑖=2 for the outer cylinder surface, resp.), and 𝑘^∗_𝑖 the reaction rate at the surface 𝑖. Between the two cylinder surfaces, we could connect a resistor and/or apply a potential difference, but in this study we only consider a short circuit. So, the boundary conditions for the potential read 𝜙^∗=0at𝑟^∗=𝑅₁,𝑟^∗=𝑅₂. (4) We may apply no-slip and impermeable conditions on the solid walls surrounding the fluid. Boundary conditions on the upper and lower ends of the domain will be addressed after dimensionless equations are presented.

The fluid flow within the annulus can be assumed to be created by two kinds of forcing; one is by the rotation of the inner cylinder and the other by the so-called induced charge electroosmotic effect. While the former is driven by the boundary condition for the Navier-Stokes equation (1c), the latter comes from the last term in (1c), that is, the Coulomb-force term, which couples the fluid flow and the charge transport problem. Apparently, the fluid velocity created by the two effects contributes to the convection of the charge, that is, the second term in (2).

As a first step in our series of studies on the charge transport within an annulus, we in this paper focus on steady and axisymmetric solutions. Then the azimuthal component of the fluid velocity does not contribute to the charge transport. Thus, even the circular Couette flow driven by the inner cylinder’s rotation has no effect on the charge transport when it is stable, which is relevant at low Reynolds numbers exhibiting only the primary azimuthal flow (referred to as steady circular Couette flow; see, e.g., Liao et al. [32]); here, the Reynolds number may be based on the tangential velocity of the inner cylinder and the gap between the two cylinders, 𝑑≡𝑅₂−𝑅₁. As the Reynolds number is increased, the primary flow becomes unstable and creates a secondary flow in the axial plane (referred to as “steady axisymmetric Taylor vortex flow”; Liao et al. [32]), which now plays an important role in the convective transport of charge. It is this type of flow that is used in the calculation of the convective terms in the charge transport equation (1a); for more complex supercritical flow regimes in the circular Couette flow system, such as nonaxisymmetric flow and travelling waves, one can refer to Koschmieder [33].

In the present study, convection due to the induced charge electroosmotic flow effect is assumed to be negligible compared with the effect of the secondary flow caused by the flow instability mentioned above following El-Adawy et al. [27]. This is valid in particular at low diffusivity of the charge carriers, 𝐷. Thus, we can decouple the fluid flow from the charge transport problem and are allowed to impose an arbitrary velocity field (but without losing the physical relevance, of course).

Based on the above reasoning, we can write the governing equations for the dimensionless charge density 𝑞 and the dimensionless potential 𝜙 in dimensionless form as follows:−1𝑟𝜕𝑟𝜕𝑟𝜕𝑞+𝜕𝜕𝑟²𝑞𝜕𝑧²1+𝑃𝑟𝜕+𝜕𝜕𝑟𝑟𝑢𝑞1𝜕𝑧𝑤𝑞+𝑄𝑞=0,𝑟𝜕𝑟𝜕𝑟𝜕𝜙+𝜕𝜕𝑟²𝜙𝜕𝑧²=−𝑞. (5a)(5b)

Here, (𝑟,𝑧)=(𝑟^∗,𝑧^∗)/𝑑, (𝑢,𝑤)=(𝑢^∗,𝑤^∗)/𝑈, 𝑞=𝑞^∗/𝑞^∗_1𝑤, 𝜙=𝜙^∗/[𝑑²𝑞^∗_1𝑤/𝜀], and 𝐄=−∇𝜙=𝐄^∗/[𝑑𝑞^∗_1𝑤/𝜀] are dimensionless variables. Further, we define 𝐉=𝐉^∗/(𝐷𝑞^∗_1𝑤/𝑑). The dimensionless parameters 𝑃=𝑈𝑑/𝐷 and 𝑄=𝜎𝑑²/(𝜀𝐷) may be considered as the ratio of time scales; 𝑃=𝜏_dif/𝜏_conv and 𝑄=𝜏_dif/𝜏_cond, where 𝜏_dif≡𝑑²/𝐷 is the diffusion time scale, 𝜏_conv≡𝑑/𝑈 the convection time scale, and 𝜏_conv≡𝜀/𝜎 the conduction time scale.

Boundary conditions for 𝑞 can be written in terms of the dimensionless radial component of the charge flux, 𝐽_𝑟≡−𝜕𝑞/𝜕𝑟−𝑄𝜕𝜙/𝜕𝑟, as follows:𝐽_𝑟=−𝑘₁𝑞−1at𝑟=𝑟₁,𝐽_𝑟=𝑘₂𝑞−𝑞_2𝑤at𝑟=𝑟₂, (6a)(6b)

where 𝑘_𝑖=𝑘^∗_𝑖/(𝐷/𝑑), 𝑟_𝑖=𝑅_𝑖/𝑑, and 𝑞_2𝑤=𝑞^∗_2𝑤/𝑞^∗_1𝑤; note that 𝑞_1𝑤=1. Boundary conditions for 𝜙 are simply 𝜙=0at𝑟=𝑟₁,𝑟=𝑟₂. (7)

As for the conditions on the upper and lower boundaries we apply zero gradient for 𝑞 and 𝜙:𝜕𝑞=𝜕𝑧𝜕𝜙𝜕𝑧=0at𝑧=0,𝑧=ℎ, (8) where ℎ corresponds to the dimensionless height of the domain, that is, the recirculating flow cell, whose dimensional quantity is defined as 𝐻; that is, ℎ=𝐻/𝑑. In this study ℎ is set as 𝐻=𝑑 or ℎ=1.

The fact that 𝜙 does not explicitly appear in (5a) may mislead to the conclusion that (5a) is decoupled from (5b), but it is not the case because it appears in boundary conditions (6a) and (6b).

3. Analytic Solutions of 1D Transport Equation

When the convection effect is neglected, we may well assume, in view of the boundary conditions, that solutions are independent of 𝑧. Then the system of equations reduces to a one-dimensional transport problem, where only the diffusion and conduction effects are present:𝑑²𝑞𝑑𝑟²+1𝑟𝑑𝑞𝑑𝑑𝑟−𝑄𝑞=0,₂𝜙𝑑𝑟²+1𝑟𝑑𝜙𝑑𝑟+𝑞=0. (9a)(9b)

The general solution to (9a) takes the following form: 𝑞=𝐶₁𝐼₀𝑟𝜆+𝐶₂𝐾₀𝑟𝜆, (10) where 𝐼₀ and 𝐾₀ denote the modified Bessel function of the first and second kind, respectively, of order zero, and 𝜆=1/𝑄 corresponds to the dimensionless thickness of EDL adjacent to the cylinder surfaces where charge is mainly distributed. Eliminating the third terms on the left-hand side of (9a) and (9b) and integrating the result twice with respect to 𝑟 yield 𝜙=−𝜆²𝑞−𝐶₃ln𝑟+𝐶₄. (11) Four constants, 𝐶₁ through 𝐶₄, in (10) and (11) can be obtained from the boundary conditions (6a), (6b), and (7) as follows: 𝐶₁=−𝐾₀₂+𝐾₀₁𝑞_2𝑤+𝐾₀₂/𝑠₁+𝐾₀₁/𝑠₂/𝜆²𝐶₃𝐼₀₂𝐾₀₁−𝐼₀₁𝐾₀₂,𝐶₂=𝐼₀₂−𝐼₀₁𝑞_2𝑤−𝐼₀₂/𝑠₁+𝐼₀₁/𝑠₂/𝜆²𝐶₃𝐼₀₂𝐾₀₁−𝐼₀₁𝐾₀₂,𝐶₃=𝑠₁𝑠₂𝜆²1−𝑞_2𝑤𝑠₁𝑠₂ln𝑟₂/𝑟₁+𝑠₁+𝑠₂,𝐶₄=𝜆²𝑠₁𝑠₂ln𝑟₂/𝑟₁+𝑠₁+𝑠₂𝑞_2𝑤𝑠₁𝑠₂ln𝑟₂/𝑟₁+𝑠₁+𝑠₂, (12) where 𝑠_𝑖=𝑘_𝑖𝑟_𝑖, 𝐼_0𝑖=𝐼₀(𝑟_𝑖/𝜆), and 𝐾_0𝑖=𝐾₀(𝑟_𝑖/𝜆) are constants.

For the case with 𝜆≪1, we are allowed to approximate the governing equations and obtain the solutions in terms of more familiar functions. We let 𝑞=𝑓(𝑟)/√𝑟 and then (9a) becomes 𝑑²𝑓𝑑𝑟²−1𝜆²−14𝑟²𝑓=0. (13) Under the assumption of 𝜆≪𝑟 (thin-layer approximation), we can ignore the second term in the bracket. Solving the resultant equation is straight forward, and we arrive at 1𝑞=√𝑟𝑞₁√𝑟₁sinh𝑟₂−𝑟/𝜆()+𝑞sinh1/𝜆₂√𝑟₂sinh𝑟−𝑟₁/𝜆(),sinh1/𝜆 (14) where the unknown constant 𝑞_𝑖 stands for the value of 𝑞 at 𝑟=𝑟_𝑖. Solution (14) is shown to be the same as that provided in [10]. In fact, the solution form (14) can also be derived from the leading order terms in the asymptotic expansion of 𝐼₀ and 𝐾₀ in (10) for large argument [34]. Equation (11) is still used for obtaining 𝜙 after 𝑞 is obtained from (14). Applying boundary conditions (6a), (6b), and (7), we get𝑞₁=1−1−𝑞_2𝑤𝑠₁𝐵,𝑞₂=𝑞_2𝑤+1−𝑞_2𝑤𝑠₂𝐵, (15a)(15b)

where 𝐵=(𝑠₁+𝑠₂)/(𝑠₁𝑠₂)+ln(𝑟₂/𝑟₁).

4. Numerical Solutions of Full 2D Transport Equations

When the velocity field is arbitrarily imposed, we must use numerical methods to solve the 2D charge transport problem governed by (5a) and (5b) under the boundary conditions (6a)–(8). We performed two kinds of numerical simulation, one using a self-developed (in-house) code and the other using the commercial code COMSOL.

We briefly address first the numerical method employed in the in-house code. Although the steady-state solutions are our primary concern, we add the transient term to the left-hand side of (5a) in order to facilitate the relaxation method. There are two key factors in developing the numerical schemes which must be borne in mind for the successful run of the simulation, numerical instability and accuracy. When convection is dominant, as is common in engineering applications, the convection terms themselves become the primary cause of the numerical instability when the central difference schemes are used. For this reason, in this study we employed the second-order upwind method to discretize the convection terms. Secondly, numerical accuracy can be maintained by constructing a variable grid system. Another key factor affecting numerical accuracy turns out to be the conservation property of the numerical schemes employed for discretization of the equations. Using the conservative form of the convection terms and employing the finite volume method in discretization turns out to be of the uttermost importance for maintaining the numerical accuracy.

Since 𝜆 usually remains very small, we expect thin layers, that is, EDLs, near the surfaces 𝑟=𝑟₁ and 𝑟=𝑟₂. To resolve such thin layers, we construct fine grids there. Along the radial direction, for instance, we use the algebraic function 𝑟=𝑟₁+3+2𝑎_𝑟𝜉(3−2𝜉)𝜉3+2𝑎_𝑟 (16) for the range 0≤𝜉≤1, where 𝑎_𝑟 is a control parameter for the variable grid along the radial direction; the case with 𝑎_𝑟=0 corresponds to no grid refinement and a larger value of 𝑎_𝑟 means finer grids near the surfaces. Usually we take 𝑎_𝑟=3 or larger. We use a similar function for the variable grids along the 𝑧 direction, where the refinement of the grids is controlled by 𝑎_𝑧, which is usually taken as 3.

Both 𝑞 and 𝜙 are defined at the point “0” (to be referred to as 𝑞_𝑖,𝑗 and 𝜙_𝑖,𝑗), the center of the grid cell being of the size Δ𝑟⁰_𝑖×Δ𝑧⁰_𝑗 as shown in Figure 1. Before discretizing the governing equation (5a), we first integrate it over the grid cell. Then we use the central difference scheme in discretization of the first-order derivatives for the diffusion terms at the surrounding four points denoted as “𝑒,” “𝑛,” “𝑤,” and “𝑠” in Figure 1. We then need to evaluate the values of 𝑞 at those four points arising from the convection terms. As mentioned before, choosing upwind methods in that evaluation is very important in establishing numerical stability. In this work, we employ the second-order upwind algorithm. For instance, when the radial component of the velocity at the point “𝑒,” 𝑢_𝑒, is positive, we construct the second-order polynomial in 𝑟 with 𝑞_{𝑖−1,𝑗}, 𝑞_𝑖,𝑗, and 𝑞_𝑖+1,𝑗 and evaluate the result at the point “𝑒” to get 𝑞_𝑒=−Δ𝑟^𝑏_𝑖²𝑞_{𝑖−1,𝑗}4Δ𝑟^𝑏_𝑖−1Δ𝑟^𝑏_𝑖−1+Δ𝑟^𝑏_𝑖+2Δ𝑟^𝑏_𝑖−1+Δ𝑟^𝑏_𝑖𝑞_𝑖,𝑗4Δ𝑟^𝑏_𝑖−1+2Δ𝑟^𝑏_𝑖−1+Δ𝑟^𝑏_𝑖𝑞_𝑖+1,𝑗4Δ𝑟^𝑏_𝑖−1+Δ𝑟^𝑏_𝑖. (17) Similarly, when 𝑢_𝑒 is negative, the polynomial is constructed with 𝑞_𝑖,𝑗, 𝑞_𝑖+1,𝑗, 𝑞_𝑖+2,𝑗, and so on. The Poisson equation (5b) for 𝜙 is treated in the same way as for the diffusion terms in (5a). We use the backward Euler method to discretize the transient term in (5a), which corresponds to the simplest stable algorithm. Since the accuracy of the solutions is independent of the time step Δ𝑡 chosen so as to make the solutions converged, we take the larger time step if possible to speed up the calculation.

Figure 1 

Notation for grid sizes, coordinates, and points in the variable grid system.

The two algebraic systems of equations constructed in this way are solved by using the SOR (successive-over relaxation) method in a coupled manner. Relaxation parameter for (5a) is usually taken as smaller than that for (5b) due to the convection terms.

In the use of the commercial software COMSOL, we employ two models, “transport of a diluted species” and “electrostatics.” The original form of the model however leads to numerical instability due to the fact that the conductivity is set to be proportional to the charge density in the original model, whereas in this study the conductivity is set to be constant. So, we modified the model in such a way that the conduction term is excluded from the charge flux 𝐉 but treated as a source. Boundary conditions (3a) and (3b), written in terms of the flux, must also be modified for the same reason. On the other hand, COMSOL allows us to use not only the conservative but also the nonconservative form for the convection terms. We will see that the nonconservative form leads to significant errors compared with the conservative form. The grid system is constructed on the same principle as applied in developing the in-house code; that is, fine grids are built near the cylinder surfaces to resolve thin EDLs.

The velocity field we are interested in is the secondary Taylor-vortex flow observed in the axial plane caused by hydrodynamic instability. Instead of using the exact solution of the secondary flow given from the numerical simulation of the Navier-Stokes equations or the experimental measurement, we set the flow in an arbitrary manner but with physical relevance if possible. For this, we assume that the axial plane between the coaxial cylinders is occupied by the series of spatially periodic flow cells. The velocity components 𝑢 and 𝑤 then can be written in a separation-of-variable form like (𝑟)𝑢=𝜋𝜓ℎ𝑟cos𝜋𝑧ℎ1,𝑤=−𝑟𝑑𝜓(𝑟)𝑑𝑟sin𝜋𝑧ℎ, (18) where 𝜓 is a kind of stream function for the axisymmetric velocity field. To determine the functional form of 𝜓(𝑟), we assume that 𝜓(𝑟) is a quadratic function (or 𝑟𝑤 is a linear function) of 𝑟 in the bulk, while 𝑟𝑤 exponentially approaches zero as 𝑟→𝑟₁ or 𝑟→𝑟₂. Then we let 𝜓(𝑟)=𝑐₀+𝑐₁𝑟+𝑐₂𝑟²+𝑐₃exp−𝛽₁𝑟−𝑟₁+𝑐₄exp−𝛽₂𝑟₂−𝑟, (19) where five unknown constants 𝑐₀--𝑐₄ are determined from the four restrictions 𝜓𝑟₁=𝜓𝑟₂=0,𝑑𝜓𝑟𝑑𝑟¹=𝑑𝜓𝑟𝑑𝑟²=0, (20) and the normalization condition for 𝑤 on 𝑧=ℎ/2, 𝜓𝑟_𝑚𝑟=−^2𝑚−𝑟²₁2. (21) Here 𝑟_𝑚 denotes the radial coordinate of the point in the bulk on 𝑧=ℎ/2 where 𝑤 vanishes; that is, 𝑑𝜓𝑟𝑑𝑟^𝑚=0. (22) Note that the spatially averaged vertical velocity component at 𝑧=ℎ/2 (averaged over 𝑟₁≤𝑟≤𝑟_𝑚) is now 1. The two parameters 𝛽₁ and 𝛽₂ control the boundary-layer thickness near the cylinder walls, where the steep distribution of 𝑤 along the normal to the wall is expected; larger 𝛽_𝑖 means a thinner layer. Figure 2 illustrates typical profiles of the velocity component 𝑤(𝑟,ℎ/2) for different sets of 𝛽₁ and 𝛽₂. It clearly shows that larger 𝛽_𝑖 yields a thinner layer and each wall-layer thickness can be controlled separately.

Figure 2 

Typical profiles of the velocity component 𝑤 on 𝑧=ℎ/2 given with three different sets of 𝛽1 and 𝛽2.

We have also prepared 1D code applicable to the case where no fluid motion exists so that the convection terms vanish. The numerical schemes are identical to those employed in the 2D code except that the variables’ dependence on 𝑧 has been removed in the 1D code.

5. Results and Discussion

The standard parameter set is given as follows [10]: 𝑅₁=76.2[mm],𝑅₁=101.6[mm],𝜀=2×8.85×10⁻¹²[],𝑘s⋅S/m_∗1≡𝑘^∗₁₀=2.84×10⁻⁶[],𝑘m/s_∗2≡𝑘^∗₂₀=3.12×10⁻⁶[m/s],𝜎=𝜎₀≡20×10⁻¹²[],𝑞S/m_∗1𝑤=7.71×10⁻³C/m³,𝑞^∗_2𝑤=1.84×10⁻³C/m³. (23) The diffusivity is set at 𝐷=8.86×10⁻¹¹ [m²/s] in Washabaugh [10]. However, in laminar flow with such a very low diffusivity, numerical simulation in general requires a very long time and the EDLs near the cylinder surfaces are too thin to be effectively resolved by a reasonably fine grid system. It also turned out that the conservation property cannot be established with such low diffusivity. Since we are concerned with laminar flow in this study, we assume much higher diffusivity than the original value, usually in the range 𝑂(10⁻⁵)∼𝑂(10⁻⁶) [m²/s]. So, in this study we set 𝐷₀≡10⁻⁶ [m²/s] as the standard diffusivity. Increasing the diffusivity is equivalent to decreasing the geometric scale as can be seen from the definition of the two main dimensionless parameters, 𝑃 and 𝑄 appearing in (5a). Other parameters to be varied here are 𝑘^∗_𝑖 and 𝜎. So, we set the following: 𝐷=𝛼_𝐷𝐷₀, 𝑘^∗_𝑖=𝛼_𝑘𝑖𝑘^∗_𝑖0, and 𝜎=𝛼_𝜎𝜎₀, where 𝛼_𝐷, 𝛼_𝑘𝑖, and 𝛼_𝜎 are multiplying factors. The reference velocity will vary in the range 𝑈=10⁻⁶∼1 [m/s], and the control parameters for the velocity profiles will be set at 𝛽₁=𝛽₂=20, except where otherwise mentioned.

Figure 3 shows the sensitive dependence of the numerical solutions given by the COMSOL simulations on the conservation property of the convection terms and the form of the source term. The diffusivity factor is set at 𝛼_𝐷=88.6 in the simulations; then, we get 𝑃=287 and 𝑄=8.23. Here CC means that the convection terms are treated with the conservative form, that is, with the form shown in (5a) (second term). CN means that the convection terms are treated with the nonconservative form, that is, with 𝑃[𝑢𝜕𝑞/𝜕𝑟+𝑤𝜕𝑞/𝜕𝑧]. Further, SC means that the source term is written as 𝑄𝑞 (like the third term in (5a)), whereas SE means that it is written as 𝑄∇⋅𝐄. When the nonconservative form for the convection terms (CN) is used, the results vary enormously with 𝛼_𝑘 (≡𝛼_𝑘1=𝛼_𝑘2), and at low values of 𝛼_𝑘 it can even be negative, which is physically irrelevant. We also conducted the grid-dependence test with different grids. At 𝛼_𝑘=0.1, the in-house code gives 𝑞_1𝑐=0.534 with grids 51×51, whereas it gives 𝑞_1𝑐=0.536 with 101×101, showing a very small change. The scheme CC-SE in COMSOL on the other hand gives 𝑞_1𝑐=0.399 with the number of grid elements 12,000, but it gives 𝑞_1𝑐=0.502 with 33,000 elements and 𝑞_1𝑐=0.524 with 61,000 elements, indicating that as the grids are refined the data approaches the value given by the in-house code. We also confirm that the in-house code yields the same result regardless of 𝛼_𝑘. We also developed an in-house code which uses the nonconservative convection terms. For the same parameter set given above and for 𝛼_𝑘=1, we get 𝑞_1𝑐=0.025 with the grids 51×51 and 𝑞_1𝑐=0.113 with the grids 101×101, both being much smaller than the correct value 𝑞_1𝑐=0.534. This implies that the conservative property of the convection terms is one of the most important factors regarding numerical accuracy.

Figure 3 

Dependence of 𝑞1𝑐, the value of 𝑞1 at 𝑧=ℎ/2, on the conservation property of the convection terms and the specific form of the source terms in the charge transport equation given from COMSOL simulation with the number of grid elements being 12,000. The diffusivity factor is set at 𝛼𝐷=88.6 and the reference velocity at 𝑈=1 [m/s]. Solid lines with square symbols are given by CN-SC (convection terms in a nonconservative form and the source term written in terms of the charge), solid lines with circles by CN-SE (convection terms in a nonconservative form and the source term written in terms of the electric field), solid line with deltas by CC-SC (convection terms in a conservative form and the source term written in terms of the charge), and solid line with diamond symbols by CC-SE (convection terms in a conservative form and the source term written in terms of the electric field). The solid line without symbols, 𝑞1𝑐=0.534, denotes the numerical solution obtained by the 2D in-house code.

The fact that the numerical solutions are sensitively dependent on the form of the convection or conduction (source) term in the governing equations implies that a small error in the equations can yield a significantly different solution. In order to explore the reason, we perform a simple analysis with the one-dimensional equation without convection effect:𝑞^′′−𝜆²𝜙𝑞=𝛿,_′′+𝑞=0, (24a)(24b)

where the prime denotes differentiation with respect to the new variable 𝑥=𝑟−𝑟₁ and 𝛿 is an arbitrary small error which is supposed to be contained in the charge transport equation due to the use of different forms of each term. In fact, the above equations can be derived from (9a) and (9b) under the limit of very large cylinders, 𝑟→∞. The boundary conditions to be satisfied are𝐽_𝑟≡−𝑞^′+𝜆²𝜙^′=−𝑘₁−𝑞𝑞−1at𝑥=0,_′+𝜆²𝜙^′=𝜅𝑘₁𝑞−𝑞_2𝑤at𝑥=1,𝜙=0at𝑥=0,1. (25a)(25b)(25c)

Here, 𝑘₁ is taken to be very small, and 𝜅≡𝑘₂/𝑘₁=𝑂(1); we calculate 𝑘₁=8.1×10⁻⁴ with the standard parameter set and 𝐷=8.86×10⁻⁵ [m²/s]. Then, after some algebraic work we derive 𝑞₁=𝑞₂=1+𝜅𝑞_2𝑤−𝛿/𝑘₁.1+𝜅 (26) This gives 𝑞₁=𝑞₂=0.601 for 𝛿=0 (no error), which is not so much different from 𝑞_1𝑐=0.534 given from the two-dimensional simulation (see Figure 3). More importantly, we see from (26) that even a small value of 𝛿 can bring a significantly different value of 𝑞₁ as long as 𝑘₁ remains small. It can even produce a negative value of 𝑞₁ when 𝛿>𝑘₁(1+𝜅𝑞_2𝑤).

On the other hand, the numerical data of 𝑞_1𝑐 given from the 2D in-house code with the conservative convection terms are almost invariant of the parameter 𝑘₁ as shown in Figure 3. This can be understood from the expansion of (15a) for small 𝑘₁ reading 𝑞₁=1−1−𝑞_2𝑤𝜅𝑟₂ln𝑟₂/𝑟₁𝑟₁+𝜅𝑟₂+𝑂𝑘₁. (27) Equation (27) clearly indicates that 𝑞₁, which can also be shown to be the same as 𝑞₂ for small 𝑘₁, is independent of 𝑘₁ to the leading order. In passing, evaluation of (27) gives 𝑞₁=0.547, which differs only 2% from the 2D result 𝑞_1𝑐=0.534 (see Figure 3). This agreement is remarkable considering that in 2D simulations 𝑞₁ is also a function of 𝑧.

Figure 4(a) shows the distribution of 𝑞(𝑟) and 𝑄𝜙(𝑟) obtained numerically from the 1D code for the case without convection effect; the graph is almost indistinguishable from the one given by the analytical solutions (11) and (14). We confirm an almost linear relationship between 𝑞 and 𝜙 as described in (11); 𝐶₃ is small because 𝑘₁ is small. Since the charge is positive everywhere, the second derivative of the potential is negative, as can be seen from (24b), being consistent with Figure 4(a). Figure 4(b) shows the distribution of the total amount of the charge flux times the radius, (𝑟𝐽_𝑟)_tot, and the two contributions, (𝑟𝐽_𝑟)_dif and (𝑟𝐽_𝑟)_cond obtained from the 1D code. Since the charge is large near both walls and small in between the two, the diffusion must occur from the walls to the central region so that (𝑟𝐽_𝑟)_dif>0 on the left-hand side and (𝑟𝐽_𝑟)_dif<0 on the right-hand side as shown in Figure 4(b). On the other hand, since the electric field is directed from the central region toward the walls, the charge receives Coulomb force in the same direction as the field vector, and so the sign of (𝑟𝐽_𝑟)_cond is the reverse of that of (𝑟𝐽_𝑟)_dif as shown in Figure 4(b). We will see below that those two can make balance with each other independently of the charge input or output through the walls. First, we note from Figure 4(b) that (𝑟𝐽_𝑟)_tot is much smaller than the other two contributions. As addressed before, this is caused by the smallness of 𝑘_𝑖. Since (𝑟𝐽_𝑟)_tot is positive, the charge is transmitted from the inner to the outer cylinder side, which is physically correct because 𝑞_1𝑤=1 is higher than 𝑞_2𝑤. This implies that for a small value of 𝑘_𝑖 the diffusive charge flux is balanced by the conductive flux. We confirm from the 1D simulation that setting 𝑞_2𝑤=𝑞_1𝑤=1 produces (𝑟𝐽_𝑟)_tot=0 although the distributions of (𝑟𝐽_𝑟)_dif and (𝑟𝐽_𝑟)_cond are similar to those of Figure 4(b). This means that the nonzero charge distribution can be established as long as the wall charge is assigned with a nonzero value.

(a)

(b)

Figure 4 

Numerical results of 1D simulation with 𝛼𝐷=88.6; (a) distribution of 𝑞(𝑟) (solid line) and 𝑄𝜙(𝑟) (dashed line) and (b) distribution of (𝑟𝐽𝑟)dif (dashed line), (𝑟𝐽𝑟)cond (dash-dot line), and 1000(𝑟𝐽𝑟)tot (solid line).

Two important dimensionless parameters explicitly appearing in the governing equations (5a) and (5b) are 𝑃 and 𝑄, and exploring the influence of these parameters on the solutions’ behavior is the main purpose of this study. The former represents the importance of the convection terms compared with the diffusion terms, whereas the latter is related to the thickness of EDL as 𝜆=1/𝑄. Under the standard parameter set and 𝑈=1 [m/s], we calculate 𝑃=2.54×10⁴ and 𝑄=729. In the following parametric study, we will start with small values of 𝑃 and 𝑄, and see how the solution structures vary with increase of these parameters.

Figure 5 shows the numerical results of the 2D code given at 𝛼_𝜎=0.001 and 𝑈=10⁻⁵ [m/s] at which 𝑃=0.254 and 𝑄=0.729. 𝑃 remains small enough that the distributions of 𝑞 (Figure 5(a)) and 𝜙 (Figure 5(b)) are almost one-dimensional being invariant of 𝑧. The diffusive charge-flux vector, 𝐉_dif≡−∇𝑞 (Figure 5(c)), and the conductive charge-flux vector, 𝐉_cond≡−𝑄∇𝜙 (Figure 5(d)), also reveal a one-dimensional nature; they are almost heading for the radial direction and are balanced by each other. On the other hand, the convective charge-flux vector, 𝐉_conv≡𝑃𝐮𝑞 (Figure 5(e)), is almost in the same pattern as the velocity vector, 𝐮, because the charge distribution is nearly uniform as shown in Figure 5(a). The total charge-flux vector, 𝐉_tot≡𝐉_dif+𝐉_conv+𝐉_cond (Figure 5(f)), then looks not so much different from 𝐉_conv, because although 𝐉_dif and 𝐉_cond are in the same order of magnitude as 𝐉_conv, they combine to become much smaller than each one, as discussed above. The stream traces of the total charge flux (Figure 5(f)) show that the wall charge at 𝑟=𝑟₁ is diffused into the domain, convected along the streamlines, and finally diffused into the wall at 𝑟=𝑟₂. In the bulk, the total charge-flux vector plot shows a recirculating pattern like the fluid velocity vector. Figure 6 shows a sketch of the passage through which the charge is transmitted from the inner wall at 𝑟=𝑟₁ to the outer wall at 𝑟=𝑟₂ together with the definition of the thicknesses of the passage, 𝛿_𝐿𝑐 and 𝛿_𝑇𝑐 at particular places.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5 

Numerical results of 2D simulation at 𝛼𝜎=10−3 and 𝑈=10−5 [m/s]; (a) contours of 𝑞, (b) contours of 𝜙, (c) 𝐉dif, (d) 𝐉cond, (e) 𝐉conv, and (f) 𝐉tot.

Figure 6 

Passage of the charge transmission from the inner cylindrical wall at 𝑟=𝑟1 to the outer cylindrical wall at 𝑟=𝑟2. Also shown are the thickness of the passage at 𝑧=ℎ/2 on the inner wall, 𝛿𝐿𝑐, and that at 𝑟=𝑟𝑐 on the top wall, 𝛿𝑇𝑐.

At a reference velocity 10 times higher, 𝑈=10⁻⁴ [m/s], we get 𝑃=2.54 and 𝑄=0.729. Distribution of 𝑞 (Figure 7(a)) shows deviation from the one-dimensional structure, while 𝜙 (Figure 7(b)) still keeps its one-dimensional nature. Accordingly, 𝐉_dif is somewhat deteriorated (Figure 7(c)), but 𝐉_cond is still heading for the radial direction (Figure 7(d)). The sum of the two flux vectors, 𝐉_dif+𝐉_cond (Figure 7(e)), now shows much decreased level caused by their mutual balance but with complex stream-trace structure. The total flux vector, 𝐉_tot (Figure 7(f)), then shows that most of the bulk region is composed of the recirculating pattern and simultaneously the passage near the walls responsible for the wall-to-wall charge transmission is more narrowed than in Figure 5(f) with smaller 𝑈.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7 

Numerical results of 2D simulation at 𝛼𝜎=10−3 and 𝑈=10−4 [m/s]; (a) contours of 𝑞, (b) contours of 𝜙, (c) 𝐉dif, (d) 𝐉cond, (e) 𝐉dif+𝐉cond, and (f) 𝐉tot.

Further increase of 𝑈 (at 𝑈=10⁻³ [m/s]; 𝑃=25.4, 𝑄=0.729) tends to make the contours of 𝑞 the closed-curve style in particular in the bulk (Figure 8(a)). Thus the diffusive flux vectors (Figure 8(c)) focus on the center of the contours of 𝑞, where 𝑞 is minimized. However variation of 𝑞 is not so significant that the contours of 𝜙 still show the one-dimensional nature (Figure 8(b)). Thus, the conductive flux vectors (Figure 8(d)) are heading for the radial direction. The pattern of 𝐉_dif+𝐉_cond (Figure 8(e)) is now more complex than Figure 7(e) but its magnitude still remains at low value. The charge transmission passages near the walls at 𝑟=𝑟₁ and 𝑟=𝑟₂ become thinner than Figure 7(f) because the convective effect is more pronounced.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8 

Numerical results of 2D simulation at 𝛼𝜎=10−3 and 𝑈=10−3 [m/s]; (a) contours of 𝑞, (b) contours of 𝜙, (c) 𝐉dif, (d) 𝐉cond, (e) 𝐉dif+𝐉cond, and (f) 𝐉tot.

When 𝑈 is further increased, the contour plot of 𝑞 tends to the streamline pattern of the fluid flow and the diffusive flux vector plot shows the radial inward pattern more clearly than Figure 8(c), but the distribution of 𝜙 and the conductive flux vector plot are qualitatively the same as Figures 8(b) and 8(d). The total charge-flux vector plot also tends more closely to the convective flux vector plot, because the convective effect is more dominant at higher values of 𝑈 (or equivalently at larger values of 𝑃). The thickness of the transmission passage becomes smaller at larger 𝑈. We measure the dimensionless thickness 𝛿_𝐿𝑐 as 0.094 at 𝑈=10⁻⁵ [m/s] (from Figure 5(f)), 0.024 at 𝑈=10⁻⁴ [m/s] (from Figure 7(f)), 0.0070 at 𝑈=10⁻³ [m/s] (from Figure 8(f)), 0.0016 at 𝑈=10⁻² [m/s], and so forth. We also measure 𝛿_𝑇𝑐 as 0.127 at 𝑈=10⁻⁵ [m/s], 0.014 at 𝑈=10⁻⁴ [m/s], 0.0014 at 𝑈=10⁻³ [m/s], 0.00014 at 𝑈=10⁻² [m/s], and so forth. This indicates that the amount of charge transmitted through the walls per unit area and unit time is almost independent of 𝑈 at high values of 𝑈.

As 𝑈 is increased at a fixed 𝜎, the contour of 𝑞 tends to follow the fluid streamline pattern because, when the fluid circulates along its closed path with high velocity, the fluid circulation time, that is, 𝜏_conv, becomes so short compared with 𝜏_dif and 𝜏_cond that the diffusive or conductive action does not have enough time to modify the value of 𝑞 but must leave it to remain at a constant value specific to the circulation path. That is, we can estimate for large value of 𝑈 the relation 𝑞=𝑞(Ψ), where Ψ is the stream function reading Ψ=𝜓(𝑟)sin(𝜋𝑧/ℎ) in this study. This can be also proved from (5a). At the limit 𝑃→∞, (5a) becomes ∇⋅(𝐮𝑞)≈0 or equivalently (𝐮∇)𝑞≈0 (from the continuity equation). In terms of the material derivative, this can be written as 𝐷𝑞/𝐷𝑡=0, indicating that 𝑞=constant along the streamline, on which Ψ is set as constant.

In order to confirm the above reasoning, we calculate the Lagrangian variation of 𝑞 of a fluid particle while it flows along the given streamline. Figure 9 shows the numerical results obtained at various reference velocities and at fixed 𝛼_𝜎=10⁻³. The abscissa 𝑠 indicates the distance travelled by the fluid particle along the closed streamline specified by the initial point, normalized by the total length of the streamline. So, 𝑠=0 means the starting point and 𝑠=1 the final point, which is set to be the same as the initial point. The figure clearly demonstrates that as the reference velocity increases, the charge density variation is more uniform because the time taken for one complete circulation is decreased accordingly. It also reveals that the inner streamline starting at (𝑟,𝑧)=(𝑟_𝑐,0.2) (Figure 9(b)) shows smaller variation of 𝑞 than the outer one starting at (𝑟,𝑧)=(𝑟_𝑐,0.02) (Figure 9(a)) for the same parameter set; this is also consistent with the above reasoning because the former requires a smaller circulation time than the latter. The overall level of 𝑞 is observed to increase as 𝑈 is increased in particular near 𝑠=0.5, because the contour of 𝑞 tends to the plume structure when 𝑈 is increased (partly shown in Figure 8(a) near 𝑧=0 and 1).

(a)

(b)

Figure 9 

Lagrangian variation of the charge density versus the normalized distance travelling from the initial point of a fluid particle while flowing along the streamline, numerically given at the four reference velocities indicated and at 𝛼𝜎=10−3. Initial points are (a) (𝑟,𝑧)=(𝑟𝑐,0.02) and (b) (𝑟,𝑧)=(𝑟𝑐,0.2).

Now we investigate the effect of 𝜎 on the charge-transport behavior at fixed 𝑈=10⁻⁴ [m/s]. Figure 10 shows the numerical results given at 𝛼_𝜎=0.1 (𝑃=2.54, 𝑄=72.9). Compared with Figure 7(a) for 𝜎 being 100 times smaller, the level of 𝑞 is significantly decreased in the bulk, because 𝜆 is decreased by the factor 0.1. As a result, the level of 𝜙 is also decreased (Figure 10(b)). The vector plots of 𝐉_dif (Figure 10(c)) and 𝐉_cond (Figure 10(d)) show the one-dimensional nature and they are not qualitatively so much different from Figures 7(c) and 7(d), respectively. However, their magnitude is increased significantly, in particular near the walls, of course. The level of their sum, 𝐉_dif+𝐉_cond (Figure 10(e)), is accordingly increased but not so much as 𝐉_dif and 𝐉_cond. On the other hand, 𝐉_tot (Figure 10(f)) differs from Figure 7(f) significantly. Its magnitude is about 3 times smaller than Figure 7(f), because at 𝛼_𝜎=0.1 the charge density 𝑞 decreases sharply with the normal distance from the walls, which causes 𝐉_conv to be reduced; note that in this parameter set, too, the contribution of 𝐉_conv to 𝐉_tot is dominant over the other fluxes. On the other hand, in the central region near 𝑟=𝑟_𝑐, the level of 𝑞 is much lowered so that the levels of 𝐉_conv and 𝐉_tot must also be low there. Thus, in order to satisfy the charge conservation, the passage of the charge transmission near the top wall must be wider than that near the side walls, 𝛿_𝑇𝑐>𝛿_𝐿𝑐. In other words, at a high enough value of 𝜎, the central region is dominated by the diffusive and conductive flux but not by the convective flux, whereas the EDLs are dominated by the convective flux, except for the very thin EDLs closer to the walls where the diffusion terms should be more effective because the fluid velocity remains very small.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 10 

Numerical results of 2D simulation at 𝛼𝜎=0.1 and 𝑈=10−4 [m/s]; (a) contours of 𝑞, (b) contours of 𝜙, (c) 𝐉dif, (d) 𝐉cond, (e) 𝐉dif+𝐉cond, and (f) 𝐉tot.

As 𝛼_𝜎 is further increased, the EDLs are clearly distinguished from the bulk as shown in Figure 11(a) for the radial distribution of 𝑞 and 𝜙 at 𝑧=ℎ/2. Since 𝑞≈0 in the bulk (to be referred to as “charge depletion zone”), 𝜙 must be a logarithmic function of 𝑟 as indicated in (11) and shown in Figure 11(a). Thus the electric field is nonzero in the bulk causing 𝐉_cond≡𝑄𝜎𝐄 to be finite there. The diffusive flux and convective flux must vanish in the bulk, because 𝑞≈0 there. Then, the total flux must be dominated by the conductive flux, which can be confirmed from Figures 11(b), 11(c), and 11(f), where the patterns of stream trace of the fluxes look similar to each other near the central region, 𝑟=𝑟_𝑐. On the other hand, the EDLs are dominated by the convective flux. At this parameter set, the level of 𝐉_tot shown in Figures 11(e) and 11(f) is further decreased from Figure 10(f) because of the decreased EDL thickness. We also observe that 𝛿_𝑇𝑐 is increased significantly with increase of 𝜎. Measurement of 𝛿_𝑇𝑐 from the simulation results gives 𝛿_𝑇𝑐=0.021 at 𝛼_𝜎=0.01, 0.079 at 𝛼_𝜎=0.1, and 0.386 at 𝛼_𝜎=1, while 𝛿_𝐿𝑐 is not so much changed upon 𝜎; we measure 𝛿_𝐿𝑐=0.028 at 𝛼_𝜎=0.01, 0.034 at 𝛼_𝜎=0.1, and 0.038 at 𝛼_𝜎=1.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 11 

Numerical results of 2D simulation at 𝛼𝜎=1 and 𝑈=10−4 [m/s]; (a) radial distributions of 𝑞 and 𝜙 on 𝑧=ℎ/2, (b) 𝐉cond, (c) 𝐉dif+𝐉cond, (d) 𝐉conv, (e) 𝐉tot showing its 𝑟-profile, and (f) 𝐉tot showing its 𝑧-profile.

We can estimate the dependence of the charge density distribution on 𝜎 from the Lagrangian variation of 𝑞 for various values of 𝛼_𝜎 as shown in Figure 12. On the streamline having the initial point at (0.5,0.02), the overall level of 𝑞 decreases with 𝛼_𝜎, but the range of variation is rather increased up to 𝛼_𝜎=0.1 (Figure 12(a)). The latter result can be understood from the fact that as 𝛼_𝜎 is increased the wall charge 𝑞_𝑖 remains almost constant (we measure only 4% change for the range 10⁻³≤𝛼_𝜎≤1) but the charge depletion zone begins to appear near 𝛼_𝜎=0.1. As 𝛼_𝜎 is further increased, EDLs become thinner while the charge depletion zone becomes wider, which causes the variation of 𝑞 to be smaller as shown in Figure 12(a). We can also give a qualitatively similar description as to the Lagrangian variation of 𝑞 of the fluid particle on the streamline with the initial point at (𝑟_𝑐,0.2) (Figure 12(b)). In particular, the level of 𝑞 vanishes all the way through the streamline at 𝛼_𝜎=1 implying that charge depletion zone should encompass the region surrounded by the streamline.

(a)

(b)

Figure 12 

Lagrangian variation of the charge density versus the normalized distance travelling from the initial point of a fluid particle while flowing along the streamline, numerically given at the three indicated values of 𝛼𝜎 and at 𝑈=10−4 [m/s]. Initial points are (a) (𝑟,𝑧)=()=(𝑟𝑐,0.02) and (b) (𝑟,𝑧)=()=(𝑟𝑐,0.2).

A typical solution structure of the charge transport equations with 𝛼_𝜎 and 𝑈 sufficiently high is shown in Figure 13. We can say that the given value of 𝛼_𝜎 is high enough in that the EDLs are thin and the charge depletion zone is wide. We can also say that the given value of 𝑈 is high enough in that contours of 𝑞 in the regions near 𝑧=0 and 𝑧=1 show the plume structure, which then causes the charge to accumulate there and makes the contour of 𝜙 to be of a saddle type near the center point at (𝑟_𝑐,ℎ/2).

(a)

(b)

Figure 13 

Typical solution structure of the system of charge transport equations given numerically at 𝛼𝜎=0.2 and 𝑈=0.01 [m/s]; contours of (a) 𝑞 and (b) 𝜙.

Electrification is known to be directly related to the amount of charge accumulated in the bulk, which in this study is quantified by the averaged charge density 𝑞_av=2ℎ𝑟²₂−𝑟²_1𝐴𝑟𝑞𝑑𝑟𝑑𝑧. (28) Figure 14 shows variation of 𝑞_av obtained from the 2D code for various 𝑈 values. In the limit 𝑈→0, 𝑞_av asymptotes to a constant value, because the charge distribution tends to the 1D structure (see, e.g., Figure 5(a)). As 𝑈 increases, 𝑞_av also increases, because the fluid conveys the charge within the EDLs to the top (𝑧=1) and bottom (𝑧=0) boundaries of the domain, leading to the plume structure there (see, e.g., Figures 8(a) and 13(a)). The amount of increase however depends on the value of 𝜎. At low values of 𝛼_𝜎, the level of 𝑞 is high but its spatial variation is small (Figures 7(a) and 8(a)) so that even the plume structure brings a slight increase of 𝑞_av upon 𝑈 as shown in Figure 14(a). On the other hand, at high values of 𝛼_𝜎, thin EDLs appear distinctively and a nonzero value of 𝑞 is detected only within the EDLs. In this case, existence of the convection effect would sweep the charge carriers within EDLs toward the top and bottom boundaries, giving rise to additional thin layers of nonzero charge (i.e., plume structure) there. Since increase of 𝑈 tends to make the overall level of 𝑞 higher in those additional layers, we expect more increase of 𝑞_av upon 𝑈, as shown in Figure 14(a), than is the case with lower 𝛼_𝜎 (Figure 14(b)).

(a)

(b)

Figure 14 

Volume-averaged charge density within the annulus given numerically at (a) 𝛼𝜎=0.2 and (b) 𝛼𝜎=10−3.

6. Conclusions

We studied the physics of charge transport in an annulus between concentric circular cylinders from theoretical and numerical analysis by using a commercial software COMSOL and 2D in-house code.

We have found that the conservation property of the convective terms in the charge transport equation affects numerical accuracy significantly. In both the COMSOL and in-house code simulations, keeping the convective terms in conservative form is essential in maintaining the numerical accuracy. In COMSOL, the conductive terms being treated as sources must also be written in the gradient-of-field form, not in the form of charge so as not to deteriorate the numerical accuracy. Such sensitive dependence of the numerical solutions’ accuracy on a small error in the governing equations can be explained in terms of 1D simplified equations.

In the absence of the convection effect, the analytical and numerical solutions of the 1D equations show that the diffusive charge flux is balanced by the conductive flux and the sum of the two fluxes yields the total flux which remains much smaller than the two fluxes for small values of 𝑘_𝑖.

The effect of two dimensionless parameters, 𝑃 and 𝑄, on the solution structure is then studied from the simulation with 2D in-house code. At low values of 𝑃, the solutions tend to a one-dimensional nature being independent of 𝑧. Increase of 𝑃 caused by increase of 𝑈 tends to make the convective effect dominate over the diffusive or conductive effect and the contour of 𝑞 tends to follow the fluid’s streamline pattern, because the time taken by the fluid particle for one complete circulation along the streamline becomes so short that the diffusive or conductive action is ineffective. At low values of 𝑄, on the other hand, 𝑞 tends to be uniformly distributed, because the diffusion effect is strong compared with conduction. Accordingly, 𝑃 is less effective at low values of 𝑄; spatially averaged value of 𝑞 increases only 2.2% when 𝑈 increases from 10⁻⁵ [m/s] to 0.1 [m/s] at 𝛼_𝜎=10⁻³ fixed, whereas it increases 42% at 𝛼_𝜎=0.2 for the same increase of 𝑈.

Increase of 𝑄 caused by increase of 𝜎 makes the EDLs thinner and the charge depletion zone wider. Thickness of the passage near the inner and outer walls through which the charge is transmitted from one wall to the other also decreases as 𝑄 is increased.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by the National Research Foundation of Korea (NRF) Grant funded by the Korean government (MSIP) (no. 2009-0083510). This work was also supported by the Human Resources Development Program (no. 20114030200030) of the Korean Institute of Energy Technology Evaluation and Planning (KETEP) grant funded by the Korean government Ministry of Trade, Industry and Energy. This paper has been read by Professor M. Duffy.