#### 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:

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 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

where denotes the radial component of , the wall charge density at ( for the inner and 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 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:

Here, , , , , and are dimensionless variables. Further, we define . The dimensionless parameters and may be considered as the ratio of time scales; and , where is the diffusion time scale, the convection time scale, and the conduction time scale.

Boundary conditions for can be written in terms of the dimensionless radial component of the charge flux, , as follows:

where , , and ; note that . Boundary conditions for are simply

As for the conditions on the upper and lower boundaries we apply zero gradient for and : 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 .

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:

The general solution to (9a) takes the following form: where and denote the modified Bessel function of the first and second kind, respectively, of order zero, and 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 Four constants, through , in (10) and (11) can be obtained from the boundary conditions (6a), (6b), and (7) as follows: where , , and are constants.

For the case with , we are allowed to approximate the governing equations and obtain the solutions in terms of more familiar functions. We let and then (9a) becomes 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 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

where .

#### 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 for the range , where is a control parameter for the variable grid along the radial direction; the case with corresponds to no grid refinement and a larger value of means finer grids near the surfaces. Usually we take 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 , , and and evaluate the result at the point “” to get Similarly, when is negative, the polynomial is constructed with , , , 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.

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 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 where five unknown constants are determined from the four restrictions and the normalization condition for on , Here denotes the radial coordinate of the point in the bulk on where vanishes; that is, Note that the spatially averaged vertical velocity component at (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 for different sets of and . It clearly shows that larger yields a thinner layer and each wall-layer thickness can be controlled separately.

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]:
The diffusivity is set at [m^{2}/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 [m^{2}/s]. So, in this study we set [m^{2}/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: , , and , where , , and are multiplying factors. The reference velocity will vary in the range [m/s], and the control parameters for the velocity profiles will be set at , 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 in the simulations; then, we get and . 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 (), 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 , the in-house code gives with grids , whereas it gives with , showing a very small change. The scheme CC-SE in COMSOL on the other hand gives with the number of grid elements 12,000, but it gives with 33,000 elements and 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 , we get with the grids and with the grids , both being much smaller than the correct value . This implies that the conservative property of the convection terms is one of the most important factors regarding numerical accuracy.

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:

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

Here, is taken to be very small, and ; we calculate with the standard parameter set and [m^{2}/s]. Then, after some algebraic work we derive
This gives for (no error), which is not so much different from 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 .

On the other hand, the numerical data of 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 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 , which differs only 2% from the 2D result (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, , and the two contributions, and 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 on the left-hand side and 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 is the reverse of that of 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 is much smaller than the other two contributions. As addressed before, this is caused by the smallness of . Since is positive, the charge is transmitted from the inner to the outer cylinder side, which is physically correct because is higher than . 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 produces although the distributions of and 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)**

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 . Under the standard parameter set and [m/s], we calculate and . 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 and [m/s] at which and . 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, (Figure 5(c)), and the conductive charge-flux vector, (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, (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, (Figure 5(f)), then looks not so much different from , because although and are in the same order of magnitude as , 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)**

At a reference velocity 10 times higher, [m/s], we get and . Distribution of (Figure 7(a)) shows deviation from the one-dimensional structure, while (Figure 7(b)) still keeps its one-dimensional nature. Accordingly, is somewhat deteriorated (Figure 7(c)), but is still heading for the radial direction (Figure 7(d)). The sum of the two flux vectors, (Figure 7(e)), now shows much decreased level caused by their mutual balance but with complex stream-trace structure. The total flux vector, (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)**

Further increase of (at [m/s]; , ) 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 (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)**

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 at [m/s] (from Figure 5(f)), at [m/s] (from Figure 7(f)), at [m/s] (from Figure 8(f)), at [m/s], and so forth. We also measure as at [m/s], at [m/s], at [m/s], at [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, , becomes so short compared with and 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 in this study. This can be also proved from (5a). At the limit , (5a) becomes or equivalently (from the continuity equation). In terms of the material derivative, this can be written as , indicating that 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 . 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, means the starting point and 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 (Figure 9(b)) shows smaller variation of than the outer one starting at (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 , because the contour of tends to the plume structure when is increased (partly shown in Figure 8(a) near and 1).

**(a)**

**(b)**

Now we investigate the effect of on the charge-transport behavior at fixed [m/s]. Figure 10 shows the numerical results given at (, ). 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 (Figure 10(c)) and (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, (Figure 10(e)), is accordingly increased but not so much as and . On the other hand, (Figure 10(f)) differs from Figure 7(f) significantly. Its magnitude is about 3 times smaller than Figure 7(f), because at the charge density decreases sharply with the normal distance from the walls, which causes to be reduced; note that in this parameter set, too, the contribution of to 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 and 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)**

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 . Since 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 to be finite there. The diffusive flux and convective flux must vanish in the bulk, because 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 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 at , at , and at , while is not so much changed upon ; we measure at , at , and at .

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

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 , the overall level of decreases with , but the range of variation is rather increased up to (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 ) but the charge depletion zone begins to appear near . 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 (Figure 12(b)). In particular, the level of vanishes all the way through the streamline at implying that charge depletion zone should encompass the region surrounded by the streamline.

**(a)**

**(b)**

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 and 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 .

**(a)**

**(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 Figure 14 shows variation of obtained from the 2D code for various values. In the limit , asymptotes to a constant value, because the charge distribution tends to the 1D structure (see, e.g., Figure 5(a)). As increases, also increases, because the fluid conveys the charge within the EDLs to the top () and bottom () 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 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 upon , as shown in Figure 14(a), than is the case with lower (Figure 14(b)).

**(a)**

**(b)**

#### 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 [m/s] to [m/s] at fixed, whereas it increases 42% at 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.