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

Authors: Y. K. Suh; K. H. Baek
Journal: Journal of Applied Mathematics
 (2014)

---

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


---

## Body

## 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 C60. 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 radiiR1 and R2, respectively, caused by diffusion, convection, and electrical conduction. The governing equations for the problem can be written as follows:(1a)∂q*∂t*+∇*·J*=0,(1b)∇*·εE*=q*,(1c)∂u*∂t*+(u*·∇*)u*=-1ρ∇*p*+ν∇*2u*+q*E*,(1d)∇*·u*=0,wheret* is the time, q* the charge density, u* the fluid velocity, p* the pressure, E* the electric field related to the electric potential ϕ* as E*=-∇*ϕ*, ∇* the gradient operator, and
(2)J*=-D∇*q*+q*u*+σE*
the current density (or charge flux). Further, ρ is the fluid density, ν the kinematic viscosity of the fluid, D 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(3a)Jr*=-k1*(q*-q1w*)atr*=R1,(3b)Jr*=k2*(q*-q2w*)atr*=R2,whereJr* denotes the radial component of J*, qiw* the wall charge density at r*=Ri (i=1 for the inner and i=2 for the outer cylinder surface, resp.), and ki* the reaction rate at the surface i. 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
(4)ϕ*=0atr*=R1,r*=R2.
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, d≡R2-R1. 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, D. 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 densityq and the dimensionless potential ϕ in dimensionless form as follows:(5a)-[1r∂∂r(r∂q∂r)+∂2q∂z2]+P[1r∂∂r(ruq)+∂∂z(wq)]+Qq=0,(5b)1r∂∂r(r∂ϕ∂r)+∂2ϕ∂z2=-q.Here,(r,z)=(r*,z*)/d, (u,w)=(u*,w*)/U, q=q*/q1w*, ϕ=ϕ*/[d2q1w*/ε], and E=-∇ϕ=E*/[dq1w*/ε] are dimensionless variables. Further, we define J=J*/(Dq1w*/d). The dimensionless parameters P=Ud/D and Q=σd2/(εD) may be considered as the ratio of time scales; P=τdif/τconv and Q=τdif/τcond, where τdif≡d2/D is the diffusion time scale, τconv≡d/U the convection time scale, and τconv≡ε/σ the conduction time scale.Boundary conditions forq can be written in terms of the dimensionless radial component of the charge flux, Jr≡-∂q/∂r-Q∂ϕ/∂r, as follows:(6a)Jr=-k1(q-1)atr=r1,(6b)Jr=k2(q-q2w)atr=r2,whereki=ki*/(D/d), ri=Ri/d, and q2w=q2w*/q1w*; note that q1w=1. Boundary conditions for ϕ are simply
(7)ϕ=0atr=r1,r=r2.As for the conditions on the upper and lower boundaries we apply zero gradient forq and ϕ:(8)∂q∂z=∂ϕ∂z=0atz=0,z=h,
where h corresponds to the dimensionless height of the domain, that is, the recirculating flow cell, whose dimensional quantity is defined as H; that is, h=H/d. In this study h is set as H=d or h=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 ofz. Then the system of equations reduces to a one-dimensional transport problem, where only the diffusion and conduction effects are present:(9a)d2qdr2+1rdqdr-Qq=0,(9b)d2ϕdr2+1rdϕdr+q=0.The general solution to (9a) takes the following form:
(10)q=C1I0(rλ)+C2K0(rλ),
where I0 and K0 denote the modified Bessel function of the first and second kind, respectively, of order zero, and λ=1/Q 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 r yield
(11)ϕ=-λ2q-C3ln⁡r+C4.
Four constants, C1 through C4, in (10) and (11) can be obtained from the boundary conditions (6a), (6b), and (7) as follows:
(12)C1=-K02+K01q2w+[(K02/s1+K01/s2)/λ2]C3I02K01-I01K02,C2=I02-I01q2w-[(I02/s1+I01/s2)/λ2]C3I02K01-I01K02,C3=s1s2λ2(1-q2w)s1s2ln⁡(r2/r1)+s1+s2,C4=λ2[s1s2ln⁡(r2/r1)+s1+s2q2w]s1s2ln⁡(r2/r1)+s1+s2,
where si=kiri, I0i=I0(ri/λ), and K0i=K0(ri/λ) 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 q=f(r)/r and then (9a) becomes
(13)d2fdr2-(1λ2-14r2)f=0.
Under the assumption of λ≪r (thin-layer approximation), we can ignore the second term in the bracket. Solving the resultant equation is straight forward, and we arrive at
(14)q=1r{q1r1sinh⁡[(r2-r)/λ]sinh⁡(1/λ)+q2r2sinh⁡[(r-r1)/λ]sinh⁡(1/λ)},
where the unknown constant qi stands for the value of q at r=ri. 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 I0 and K0 in (10) for large argument [34]. Equation (11) is still used for obtaining ϕ after q is obtained from (14). Applying boundary conditions (6a), (6b), and (7), we get(15a)q1=1-(1-q2w)s1B,(15b)q2=q2w+(1-q2w)s2B,whereB=(s1+s2)/(s1s2)+ln⁡(r2/r1).


## 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 r=r1 and r=r2. To resolve such thin layers, we construct fine grids there. Along the radial direction, for instance, we use the algebraic function
(16)r=r1+[3+2arξ(3-2ξ)]ξ3+2ar
for the range 0≤ξ≤1, where ar is a control parameter for the variable grid along the radial direction; the case with ar=0 corresponds to no grid refinement and a larger value of ar means finer grids near the surfaces. Usually we take ar=3 or larger. We use a similar function for the variable grids along the z direction, where the refinement of the grids is controlled by az, which is usually taken as 3.Bothq and ϕ are defined at the point “0” (to be referred to as qi,j and ϕi,j), the center of the grid cell being of the size Δri0×Δzj0 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 “e,” “n,” “w,” and “s” in Figure 1. We then need to evaluate the values of q 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 “e,” ue, is positive, we construct the second-order polynomial in r with qi-1,j, qi,j, and qi+1,j and evaluate the result at the point “e” to get
(17)qe=-(Δrib)2qi-1,j4Δri-1b(Δri-1b+Δrib)+(2Δri-1b+Δrib)qi,j4Δri-1b+(2Δri-1b+Δrib)qi+1,j4(Δri-1b+Δrib).
Similarly, when ue is negative, the polynomial is constructed with qi,j, qi+1,j, qi+2,j, 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 Δt 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 fluxJ 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 componentsu and w then can be written in a separation-of-variable form like
(18)u=πψ(r)hrcos⁡(πzh),w=-1rdψ(r)drsin(πzh),
where ψ is a kind of stream function for the axisymmetric velocity field. To determine the functional form of ψ(r), we assume that ψ(r) is a quadratic function (or rw is a linear function) of r in the bulk, while rw exponentially approaches zero as r→r1 or r→r2. Then we let
(19)ψ(r)=c0+c1r+c2r2+c3exp⁡[-β1(r-r1)]+c4exp⁡[-β2(r2-r)],
where five unknown constants c0–c4 are determined from the four restrictions
(20)ψ(r1)=ψ(r2)=0,dψdr(r1)=dψdr(r2)=0,
and the normalization condition for w on z=h/2,
(21)ψ(rm)=-rm2-r122.
Here rm denotes the radial coordinate of the point in the bulk on z=h/2 where w vanishes; that is,
(22)dψdr(rm)=0.
Note that the spatially averaged vertical velocity component at z=h/2 (averaged over r1≤r≤rm) is now 1. The two parameters β1 and β2 control the boundary-layer thickness near the cylinder walls, where the steep distribution of w along the normal to the wall is expected; larger βi means a thinner layer. Figure 2 illustrates typical profiles of the velocity component w(r,h/2) for different sets of β1 and β2. It clearly shows that larger βi yields a thinner layer and each wall-layer thickness can be controlled separately.Figure 2
Typical profiles of the velocity componentw on z=h/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 onz has been removed in the 1D code.


## 5. Results and Discussion

The standard parameter set is given as follows [10]:
(23)R1=76.2[mm],R1=101.6[mm],ε=2×8.85×10-12[s·S/m],k1*≡k10*=2.84×10-6[m/s],k2*≡k20*=3.12×10-6[m/s],σ=σ0≡20×10-12[S/m],q1w*=7.71×10-3[C/m3],q2w*=1.84×10-3[C/m3].
The diffusivity is set at D=8.86×10-11 [m2/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 O(10-5)~O(10-6) [m2/s]. So, in this study we set D0≡10-6 [m2/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, P and Q appearing in (5a). Other parameters to be varied here are ki* and σ. So, we set the following: D=αDD0, ki*=αkiki0*, and σ=ασσ0, where αD, αki, and ασ are multiplying factors. The reference velocity will vary in the range U=10-6~1 [m/s], and the control parameters for the velocity profiles will be set at β1=β2=20, except where otherwise mentioned.Figure3 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 αD=88.6 in the simulations; then, we get P=287 and Q=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 P[u∂q/∂r+w∂q/∂z]. Further, SC means that the source term is written as Qq (like the third term in (5a)), whereas SE means that it is written as Q∇·E. When the nonconservative form for the convection terms (CN) is used, the results vary enormously with αk (≡αk1=αk2), and at low values of αk it can even be negative, which is physically irrelevant. We also conducted the grid-dependence test with different grids. At αk=0.1, the in-house code gives q1c=0.534 with grids 51×51, whereas it gives q1c=0.536 with 101×101, showing a very small change. The scheme CC-SE in COMSOL on the other hand gives q1c=0.399 with the number of grid elements 12,000, but it gives q1c=0.502 with 33,000 elements and q1c=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 αk. We also developed an in-house code which uses the nonconservative convection terms. For the same parameter set given above and for αk=1, we get q1c=0.025 with the grids 51×51 and q1c=0.113 with the grids 101×101, both being much smaller than the correct value q1c=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 ofq1c, the value of q1 at z=h/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 αD=88.6 and the reference velocity at U=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, q1c=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:(24a)q′′-λ2q=δ,(24b)ϕ′′+q=0,where the prime denotes differentiation with respect to the new variablex=r-r1 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, r→∞. The boundary conditions to be satisfied are(25a)Jr≡-(q′+λ2ϕ′)=-k1(q-1)atx=0,(25b)-(q′+λ2ϕ′)=κk1(q-q2w)atx=1,(25c)ϕ=0atx=0,1.Here,k1 is taken to be very small, and κ≡k2/k1=O(1); we calculate k1=8.1×10-4 with the standard parameter set and D=8.86×10-5 [m2/s]. Then, after some algebraic work we derive
(26)q1=q2=1+κq2w-δ/k11+κ.
This gives q1=q2=0.601 for δ=0 (no error), which is not so much different from q1c=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 q1 as long as k1 remains small. It can even produce a negative value of q1 when δ>k1(1+κq2w).On the other hand, the numerical data ofq1c given from the 2D in-house code with the conservative convection terms are almost invariant of the parameter k1 as shown in Figure 3. This can be understood from the expansion of (15a) for small k1 reading
(27)q1=1-(1-q2w)κr2ln⁡(r2/r1)r1+κr2+O(k1).
Equation (27) clearly indicates that q1, which can also be shown to be the same as q2 for small k1, is independent of k1 to the leading order. In passing, evaluation of (27) gives q1=0.547, which differs only 2% from the 2D result q1c=0.534 (see Figure 3). This agreement is remarkable considering that in 2D simulations q1 is also a function of z.Figure4(a) shows the distribution of q(r) and Qϕ(r) 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 q and ϕ as described in (11); C3 is small because k1 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, (rJr)tot, and the two contributions, (rJr)dif and (rJr)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 (rJr)dif>0 on the left-hand side and (rJr)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 (rJr)cond is the reverse of that of (rJr)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 (rJr)tot is much smaller than the other two contributions. As addressed before, this is caused by the smallness of ki. Since (rJr)tot is positive, the charge is transmitted from the inner to the outer cylinder side, which is physically correct because q1w=1 is higher than q2w. This implies that for a small value of ki the diffusive charge flux is balanced by the conductive flux. We confirm from the 1D simulation that setting q2w=q1w=1 produces (rJr)tot=0 although the distributions of (rJr)dif and (rJr)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.Numerical results of 1D simulation withαD=88.6; (a) distribution of q(r) (solid line) and Qϕ(r) (dashed line) and (b) distribution of (rJr)dif (dashed line), (rJr)cond (dash-dot line), and 1000(rJr)tot (solid line).
(a)
(b)Two important dimensionless parameters explicitly appearing in the governing equations (5a) and (5b) are P and Q, 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/Q. Under the standard parameter set and U=1 [m/s], we calculate P=2.54×104 and Q=729. In the following parametric study, we will start with small values of P and Q, and see how the solution structures vary with increase of these parameters.Figure5 shows the numerical results of the 2D code given at ασ=0.001 and U=10-5 [m/s] at which P=0.254 and Q=0.729. P remains small enough that the distributions of q (Figure 5(a)) and ϕ (Figure 5(b)) are almost one-dimensional being invariant of z. The diffusive charge-flux vector, Jdif≡-∇q (Figure 5(c)), and the conductive charge-flux vector, Jcond≡-Q∇ϕ (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, Jconv≡Puq (Figure 5(e)), is almost in the same pattern as the velocity vector, u, because the charge distribution is nearly uniform as shown in Figure 5(a). The total charge-flux vector, Jtot≡Jdif+Jconv+Jcond (Figure 5(f)), then looks not so much different from Jconv, because although Jdif and Jcond are in the same order of magnitude as Jconv, 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 r=r1 is diffused into the domain, convected along the streamlines, and finally diffused into the wall at r=r2. 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 r=r1 to the outer wall at r=r2 together with the definition of the thicknesses of the passage, δLc and δTc at particular places.Numerical results of 2D simulation atασ=10-3 and U=10-5 [m/s]; (a) contours of q, (b) contours of ϕ, (c) Jdif, (d) Jcond, (e) Jconv, and (f) Jtot.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
Passage of the charge transmission from the inner cylindrical wall atr=r1 to the outer cylindrical wall at r=r2. Also shown are the thickness of the passage at z=h/2 on the inner wall, δLc, and that at r=rc on the top wall, δTc.At a reference velocity 10 times higher,U=10-4 [m/s], we get P=2.54 and Q=0.729. Distribution of q (Figure 7(a)) shows deviation from the one-dimensional structure, while ϕ (Figure 7(b)) still keeps its one-dimensional nature. Accordingly, Jdif is somewhat deteriorated (Figure 7(c)), but Jcond is still heading for the radial direction (Figure 7(d)). The sum of the two flux vectors, Jdif+Jcond (Figure 7(e)), now shows much decreased level caused by their mutual balance but with complex stream-trace structure. The total flux vector, Jtot (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 U.Numerical results of 2D simulation atασ=10-3 and U=10-4 [m/s]; (a) contours of q, (b) contours of ϕ, (c) Jdif, (d) Jcond, (e) Jdif+Jcond, and (f) Jtot.
(a)
(b)
(c)
(d)
(e)
(f)Further increase ofU (at U=10-3 [m/s]; P=25.4, Q=0.729) tends to make the contours of q 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 q, where q is minimized. However variation of q 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 Jdif+Jcond (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 r=r1 and r=r2 become thinner than Figure 7(f) because the convective effect is more pronounced.Numerical results of 2D simulation atασ=10-3 and U=10-3 [m/s]; (a) contours of q, (b) contours of ϕ, (c) Jdif, (d) Jcond, (e) Jdif+Jcond, and (f) Jtot.
(a)
(b)
(c)
(d)
(e)
(f)WhenU is further increased, the contour plot of q 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 U (or equivalently at larger values of P). The thickness of the transmission passage becomes smaller at larger U. We measure the dimensionless thickness δLc as 0.094 at U=10-5 [m/s] (from Figure 5(f)), 0.024 at U=10-4 [m/s] (from Figure 7(f)), 0.0070 at U=10-3 [m/s] (from Figure 8(f)), 0.0016 at U=10-2 [m/s], and so forth. We also measure δTc as 0.127 at U=10-5 [m/s], 0.014 at U=10-4 [m/s], 0.0014 at U=10-3 [m/s], 0.00014 at U=10-2 [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 U at high values of U.AsU is increased at a fixed σ, the contour of q 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 q but must leave it to remain at a constant value specific to the circulation path. That is, we can estimate for large value of U the relation q=q(Ψ), where Ψ is the stream function reading Ψ=ψ(r)sin(πz/h) in this study. This can be also proved from (5a). At the limit P→∞, (5a) becomes ∇·(uq)≈0 or equivalently (u∇)q≈0 (from the continuity equation). In terms of the material derivative, this can be written as Dq/Dt=0, indicating that q=constant along the streamline, on which Ψ is set as constant.In order to confirm the above reasoning, we calculate the Lagrangian variation ofq 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-3. The abscissa s 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, s=0 means the starting point and s=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 (r,z)=(rc,0.2) (Figure 9(b)) shows smaller variation of q than the outer one starting at (r,z)=(rc,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 q is observed to increase as U is increased in particular near s=0.5, because the contour of q tends to the plume structure when U is increased (partly shown in Figure 8(a) near z=0 and 1).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) (r,z)=(rc,0.02) and (b) (r,z)=(rc,0.2).
(a)
(b)Now we investigate the effect ofσ on the charge-transport behavior at fixed U=10-4 [m/s]. Figure 10 shows the numerical results given at ασ=0.1 (P=2.54, Q=72.9). Compared with Figure 7(a) for σ being 100 times smaller, the level of q 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 Jdif (Figure 10(c)) and Jcond (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, Jdif+Jcond (Figure 10(e)), is accordingly increased but not so much as Jdif and Jcond. On the other hand, Jtot (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 q decreases sharply with the normal distance from the walls, which causes Jconv to be reduced; note that in this parameter set, too, the contribution of Jconv to Jtot is dominant over the other fluxes. On the other hand, in the central region near r=rc, the level of q is much lowered so that the levels of Jconv and Jtot 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, δTc>δLc. 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.Numerical results of 2D simulation atασ=0.1 and U=10-4 [m/s]; (a) contours of q, (b) contours of ϕ, (c) Jdif, (d) Jcond, (e) Jdif+Jcond, and (f) Jtot.
(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 q and ϕ at z=h/2. Since q≈0 in the bulk (to be referred to as “charge depletion zone”), ϕ must be a logarithmic function of r as indicated in (11) and shown in Figure 11(a). Thus the electric field is nonzero in the bulk causing Jcond≡QσE to be finite there. The diffusive flux and convective flux must vanish in the bulk, because q≈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, r=rc. On the other hand, the EDLs are dominated by the convective flux. At this parameter set, the level of Jtot 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 δTc is increased significantly with increase of σ. Measurement of δTc from the simulation results gives δTc=0.021 at ασ=0.01, 0.079 at ασ=0.1, and 0.386 at ασ=1, while δLc is not so much changed upon σ; we measure δLc=0.028 at ασ=0.01, 0.034 at ασ=0.1, and 0.038 at ασ=1.Numerical results of 2D simulation atασ=1 and U=10-4 [m/s]; (a) radial distributions of q and ϕ on z=h/2, (b) Jcond, (c) Jdif+Jcond, (d) Jconv, (e) Jtot showing its r-profile, and (f) Jtot showing its z-profile.
(a)
(b)
(c)
(d)
(e)
(f)We can estimate the dependence of the charge density distribution onσ from the Lagrangian variation of q 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 q 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 qi remains almost constant (we measure only 4% change for the range 10-3≤ασ≤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 q to be smaller as shown in Figure 12(a). We can also give a qualitatively similar description as to the Lagrangian variation of q of the fluid particle on the streamline with the initial point at (rc,0.2) (Figure 12(b)). In particular, the level of q vanishes all the way through the streamline at ασ=1 implying that charge depletion zone should encompass the region surrounded by the streamline.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 U=10-4 [m/s]. Initial points are (a) (r,z)=()=(rc,0.02) and (b) (r,z)=()=(rc,0.2).
(a)
(b)A typical solution structure of the charge transport equations withασ and U 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 U is high enough in that contours of q in the regions near z=0 and z=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 (rc,h/2).Typical solution structure of the system of charge transport equations given numerically atασ=0.2 and U=0.01 [m/s]; contours of (a) q and (b) ϕ.
(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(28)qav=2h(r22-r12)∫Arqdrdz.
Figure 14 shows variation of qav obtained from the 2D code for various U values. In the limit U→0, qav asymptotes to a constant value, because the charge distribution tends to the 1D structure (see, e.g., Figure 5(a)). As U increases, qav also increases, because the fluid conveys the charge within the EDLs to the top (z=1) and bottom (z=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 q 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 qav upon U as shown in Figure 14(a). On the other hand, at high values of ασ, thin EDLs appear distinctively and a nonzero value of q 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 U tends to make the overall level of q higher in those additional layers, we expect more increase of qav upon U, as shown in Figure 14(a), than is the case with lower ασ (Figure 14(b)).Volume-averaged charge density within the annulus given numerically at (a)ασ=0.2 and (b) ασ=10-3.
(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 ofki.The effect of two dimensionless parameters,P and Q, on the solution structure is then studied from the simulation with 2D in-house code. At low values of P, the solutions tend to a one-dimensional nature being independent of z. Increase of P caused by increase of U tends to make the convective effect dominate over the diffusive or conductive effect and the contour of q 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 Q, on the other hand, q tends to be uniformly distributed, because the diffusion effect is strong compared with conduction. Accordingly, P is less effective at low values of Q; spatially averaged value of q increases only 2.2% when U increases from 10-5 [m/s] to 0.1 [m/s] at ασ=10-3 fixed, whereas it increases 42% at ασ=0.2 for the same increase of U.Increase ofQ 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 Q is increased.



---

*Source: 340916-2014-03-27.xml*