Numerical Investigation of Corona Plasma Region in Negative Wire-to-duct Corona Discharge

Electrostatic precipitators make use of corona discharge phenomena to remove airborne dust particles. Exact assessment of the electric field and charge density distribution is essential to understand the particle behavior and the flow dynamics inside the electrostatic precipitators. The Poisson and charge conservation equations were solved to evaluate the electric field and charge density distributions in the negative wire-to-duct electrostatic precipitator. In this article, a novel computation method calculating the plasma region thickness was presented with the plasma region model. Instead of the conventional Kaptzov’s hypothesis, a boundary condition for the charge density was suggested as a function of applied voltage. When the computation model and the charge boundary condition above were applied to previous experiments, the results showed good agreements with the experimental data. The estimated plasma region thickness was approximately 1.5–2.5 times greater than the wire radius in the wire radius range of 0.15 mm to 1.6 mm.


INTRODUCTION
In industrial electrostatic precipitators, the space charge density (q) and electric field strength (E) distributions in the domain between electrodes are essential for calculating particle charging and trajectory, which can be used to evaluate the precipitation efficiency.In addition, the solution of q and E enables the predictions of the voltage-current (V -I) relations as well as the ionic wind characteristics.The Poisson and charge conservation equations can describe the corona discharge phenomena.As the unipolar space charge is present, the distribution of electric potential (voltage) is governed by the Poisson equation where  is the air permittivity.The charge conservation law becomes where the density of the ionic current j is defined as follows:  Corresponding author.Tel.: +82 2 2123 2821; Fax: +82 2 312 2821 E-mail address: hwangjh@yonsei.ac.kr Where μ is the ion mobility and E is the electric field strength ( ). Eq. ( 3) is valid for a stationary gas at a uniform temperature where diffusion can be neglected.Once j is obtained, I can be calculated by integrating Eq.
A range of numerical methods have been used in the past to calculate the electric field and charge density distribution in a wire-to-duct precipitation system.These include the Finite Difference Method (FDM) (McDonald et al., 1977;Lawless and Sparks, 1980), Finite Elements Method (FEM) (Cristina et al., 1991), Boundary Element Method (BEM) with the Method Of Characteristics (MOC) (Adamiak, 1994), FDM combined with MOC (Lami et al., 1997;Anagnostopoulos and Bergeles, 2002) and Finite Volume Method (FVM) (Neimarlija et al., 2009).Dirichlet conditions were imposed on the electrodes as voltage boundary conditions.In those studies, the ionization region around the wire electrode was not considered.The charge density on the wire surface was determined iteratively to maintain the electric field strength on the surface at the onset electric field strength (E o ).This procedure was based on Kaptzov's hypothesis that if a corona discharge occurs at some point of the wire and charge is injected, the electric field at that point remains at the value it takes at the onset of the corona discharge (Adamiak and Atten, 2004).
Chen and Davidson numerically examined a one dimensional corona plasma region around a wire considering Kaptzov's hypothesis and the electron activities at a given current (Chen and Davidson, 2003).They defined the plasma region as a region where ions are generated due to electron-impact reactions and reported that the ion densities remained relatively constant inside the plasma boundary layer (Chen and Davidson, 2003).
In this study, a uniform negative charge density and a constant thickness around the wire electrode was assumed in the corona plasma region.Eq. ( 1) and Eq. ( 2) were solved over the entire 2-dimensional space of a wire-toduct geometry to examine the relationship between the supplied voltage and charge density in that region, which has not been reported.At the same time, the thickness of the plasma region was calculated and discussed.For this purpose, CFD commercial code, Fluent, which utilizes the FVM, was used as the solver.Two User Defined Scalar (UDS)s were used to solve Eq. ( 1) and Eq. ( 2).

Model of Plasma Region
In a negative corona discharge model, Chen and Davidson (2003) defined the ionization boundary as the radius (r i ) where the reduced electric field (E/N) is 120 Td (1 Td = 10 -21 Vm 2 ).This corresponds to a mean electron kinetic energy of 2.78 eV, where N is the neutral gas density (Chen and Davidson, 2003;Lowke and Morrow, 1994), which can be obtained from the ideal gas law, P = Nk B T (P: pressure, k B : Boltzmann constant, T: Temperature).They also defined the plasma boundary as the radius (r p ), where the reduced electric field (E/N) is 80 Td.In the region between the wire radius (r w ) and r i , the ionization coefficient is superior to the attachment coefficient, which means electrons are generated (Chen and Davidson, 2003).Beyond the ionization boundary, the electrons are energetic enough to attach to the electronegative gas molecules (for instance O 2 ) in the region between r i and r p (Chen and Davidson, 2003).Therefore, the corona plasma region was taken as the region between r w and r p .
In this study, we assumed that the charge density was uniform in the plasma region.The onset voltage V o was evaluated by solving the Laplace equation, which is the special case of q = 0 in Eq. (1), hence the following Peek's equation (Lawless and Sparks, 1980) was satisfied.
where E o is the electric field strength in V/m on the wire surface when the corona discharge is initiated, f is a factor that accounts for wire roughness, and δ is the air density relative to 1 atm and 25°C.We defined V t as the ionization threshold voltage when the value of E/N on the wire surface reached 120 Td in the Laplacian field.

Computation Method
The main computation procedure is outlined as follows: (a) Determine the electrode geometry and the voltages supplied to the wire and plate.The wire and plate voltages are given as 0 and a positive value, respectively.
(b) Assume r p .(c) Input the negative value of the charge density in the plasma region (q i ).(d) Assuming that the charge density beyond the plasma region is an infinitesimal value, calculate the distributions of E, q, I. (e) Find the value of r p for E/N = 80 Td.(f) The calculation is completed if the difference between the assumed r p in (b) and the calculated r p in (e) is < 1 % of the assumed r p .Else, replace the assumed r p with the calculated r p and repeat calculations (d) and (e).The following relation was used iteratively as the core converging algorithm to solve Eq. ( 1) and Eq.(2) in the above procedure (d) (Cristina et al., 1991): where q k is the charge density used in the calculation step number, k, and q k+1 is the charge density used in the next step.E' and E'' are the magnitudes of electric field strength from Eq. ( 1) and Eq. ( 2), respectively.The procedure above was realized using the User Defined Function (UDF) of Fluent.

RESULTS AND DISCUSSION
The proposed computation model was applied to the single wire-to-duct geometry used in the experiments reported by Ohkubo et al. (1986) and our calculation results were compared with their experimental data.The model was then applied to a multiple wire-to-duct geometry for further validation.The ion mobility was assumed to be 1.6 s V m  / 2 (Hinds, 1999).The roughness factor (f) and air density factor (δ) presented in each experiment were used in our calculations.In the case of no experimental information for f and δ, they were assumed to be 1.0.In the following figures, equations and discussions, every charge density was taken as the absolute value of the actual negative value.

Single Wire-to-duct Geometry
The relationship between the supplied voltage and charge density in the corona plasma region was formulated based upon Ohkubo et al.'s experimental work.
Fig. 1 shows one quarter of the total geometry used in the experiments reported by Ohkubo et al. (1986).After a mesh dependence study, the number of rectangular grids on the wire surface was selected as 40.The growth rate of the grid size and maximum grid size were chosen as 1.03 and 2 mm, respectively.The following parameters were used: D = 100 mm, L = 325 mm.Fig. 2 shows that the results for the electric field strength on the plate (for q ≠ 0) agree well with Ohkubo et al.'s experimental data.Fig. 2 also shows the calculation results for the case of q = 0 (Laplacian field).Therefore, the effect of the charge density on the electric field can be seen.The calculations were also performed using different applied  voltages, wire diameters and electrode spacing.Despite the changes in these parameters, the results can be still expressed as in Fig. 2, and coincided well with the experimental data reported by Ohkubo et al. (1986).The calculation results of the normalized current density distribution (figure not shown) on the plate (cd) followed the Warburg distribution, j/jc = cos m θ (j: charge density, j c : charge density on c, tanθ = x/D) when m = 4.2, and showed good agreement with Ohkubo et al. (1986)'s experimental data.
Fig. 3 shows the changes in the reduced electric field E/N with respect to the charge density, q i , in the vicinity of the wire electrode.The radius of the wire was 0.25 mm.Both of E w /N and r p decrease with increasing q i (E w /N = 350 Td, r p -r w = 0.86 mm for q i = 0.0003 C/m 3 ; E w /N = 244 Td, r p -r w = 0.53 mm for q i = 0.0005 C/m 3 ; E w /N = 142 Td, r p -r w = 0.2 mm for q i = 0.001 C/m 3 ).The maximum charge density (q i,max ) was defined as the value when E w /N reaches 120 Td because any further electron generation could not be expected for E w /N < 120 Td.The current corresponding to q i,max was defined as the maximum current (I max ) for a given voltage.The charge density estimated in the plasma region (q i,e ) for a given supplied voltage was defined as the value when the simulated current coincided with the experimental one.The V-I max relations are plotted in Fig. 4 for two wire radii, where the maximum current was compared with the theoretical saturated current, I sat = 1.62μεV 2 /D 2 (per one plate), proposed by Sigmond (1986).The two currents had the same magnitudes.
Fig. 5 shows the V -q i,max , V -q i,e relations for the wire radii of 0.25 mm and 0.545 mm.The ionization threshold voltages (V t ) were 4.6 and 9.0 kV for r w = 0.25 mm and r w = 0.545 mm, respectively.The corona onset voltages (V o ) were 14.0 and 21.5 kV for r w = 0.25 mm and r w = 0.545 mm, respectively.The following equations could be obtained: where k and α are constants depending on the electrode geometry.α was approximately 2.0 and 1.7 for a wire radius 0.25 mm and 0.545 mm, respectively.Both equations can be unified as follows: q i,max can be obtained using the proposed computation model if only the electrode geometry and supplied voltage are given.Subsequently q i,e can be calculated by Eq. ( 8) and a prediction of the actual current is possible.

Multiple Wire-to-duct Geometry
Fig. 7 shows a multiple wire-to-duct geometry for the computation, which was used in the modeling studies by McDonald et al. (1977), Lawless and Sparks (1980), and experiments by Penney and Matick (1960).The same grid 1.7   q i,e (r w =0.545mm) q i,max (r w =0.545 mm) q i,e (r w =0.25mm) q i,max (r w =0.25 mm)

Simulation Simulation Simulation Simulation
Fig. 5. V -q i,max and V -q i,e relations for r w = 0.25 mm and 0.545 mm.Experiment (Ohkubo et al., 1986) Experiment (Ohkubo et al., 1986) Experiment (Ohkubo et al., 1986) Experiment (Ohkubo et al., 1986) Our simulation generation method as the one used in the single wire geometry was adopted.The calculated voltage distributions along lines cb and de were compared with Penney and Matick's (1960) experimental results in Figs. 8 and 9, respectively.Fig. 10 shows the simulated V -I relations, which were compared with the experimental data reported by McDonald's (1977), Lawless andSparks's (1980), andPenney andMatick's (1960).The predictions from the proposed computation model agreed well with the experimental data.Deviations were observed at low  currents.Lawless and Spark (1980) suggested that even small physical variations between the corona wires could have significant effects on the corona starting voltage.On the other hand, the onset voltage of the center and side wires was 41.2 kV and 39.3 kV, respectively, from the simulation results of the Laplacian field considering all the wires in the five wire configuration reported by Lawless and Sparks (1980) (data not shown).We suppose that the simple computation geometry presented in Fig. 7 did not reflect the side wire effects, which can cause deviations at low currents.

Plasma Region Thickness and Electric Field Strength on the Surface
Fig. 11 shows the relationships between the voltage and plasma region thickness, where the thickness was normalized to the wire radius.Overall, the thickness decreases with increasing applied voltage.This tendency becomes weak as the wire radius increases.Fig. 12 presents the relationships between the wire radius and plasma region thickness.The simulated results were compared with Takahashi et al. (1982)'s equation for the thickness.Since Takahashi et al. (1982)'s equation, t h = 1.56r w 0.65 (mm), was based on the positive discharge in the coaxial wire-to-cylinder geometry, it is believed that our simulation results predict a thicker plasma region than Takahashi et al. (1982)'s corona sheath.Most industrial ESPs are operated in as high voltage as possible to increase the precipitation efficiency and the wire is thick enough in order to protect it from irregular arcs.In this regard, we estimated the plasma region is approximately 1.5-2.5 times thicker than the wire radius in the wire radius range of 0.15 mm to 1.6 mm.The thinner plasma region was attributed to McDonald (1977)'s geometry having a lower wire-to-plate distance than other geometries.All the normalized thicknesses corresponding to q i,max ranged from 0.51 to 0.53, which coincides with the value of 0.5 = (120 Td -80 Td)/80Td.Fig. 13 shows the relationship between the supplied voltage and electric field strength on the wire surface, where the electric field strength was normalized to E o calculated by Eq. ( 4).The dashed line (E w /E o = 1.0) represents Kaptzov's hypothesis.The Kaptzov's hypothesis becomes valid as the applied voltage decreases and the wire radius increases.This coincides with Anagnostopoulos and Bergeles (2002)'s analysis in that a wire radius a few times larger than the actual one is essential for making a better prediction of the experimental data in the case of a small wire radius in the order of 0.1 mm.Also from Fig. 13, we can perceive that the wire radius should be greater than 0.5 mm to enhance E w /E o higher than 0.7.

Exponent α
The exponent α ranged from 1.5 to 2.2 for the compared experiments and it can be expressed as a function of the wire radius; α = 15 when r w > 0.8 mm, and α = 2.4 -1.1r w when r w ≤ 0.8 mm.

CONCLUSIONS
A computation model with uniform charge density around the wire electrode in the corona plasma region having a certain thickness was developed to solve the electric field and charge density distributions in a negative wire-to-duct corona discharge in air.The relationship between the applied voltage and charge density in the plasma region was estimated to be q i,e = k(V-V o ) α , while the maximum charge density in the layer was expressed as q i,max = k(V-V t ) α for a given electrode geometry.The exponent α ranged from 1.5 to 2.2 for the compared experiments.The simulation results from these relations and proposed computation procedure showed good agreement with the experimental data.
Kaptzov's hypothesis for the charge density conditions on the wire surface was valid only for the geometry of a wire radius > 0.5 mm.The plasma region thickness was approximately 1.5-2.5 times higher than that of the wire radius in the wire radius range of 0.15 mm to 1.6 mm.In addition, the calculated current density distribution in the single wire geometry followed the Warburg distribution, j/j c = cos m θ when m = 4.2.

Fig. 1 .
Fig. 1. 2-dimensional calculation domain for a single wire-to-duct geometry (not to scale).The upper schematic diagram shows the entire geometry.The shaded quarter represents the calculation domain below.Calculation domain: area abcde, wire surface: line ab, plate: line cd , Wire-to-plate gap: D (100 mm), Wire-to-edge gap : L (325 mm), f = 1.0, δ = 1.0.Voltage boundary conditions: Dirichlet for line ab and line cd, Neumann for line de, symmetry for line bc and line ae; Charge boundary conditions: constant in domain aa΄b΄b.

Fig. 2 .
Fig. 2. Normalized electric field strength on the plate for r w = 0.25 mm and V = 34.11kV.E c : electric field strength at c, x: distance from c to d along the line cd in Fig. 1.

Fig. 3 .Fig. 4 .
Fig. 3. Electric field variations around the wire electrode for r w = 0.25 mm.The intersecting points of the solid lines and dashed line represents r p .Unit for q i : C/m 3 .

Fig. 7 .Fig. 8 .
Fig. 7. 2 dimensional calculation domain for multiple wire-to-duct geometry (not to scale).The upper schematic depicts the entire geometry for the case of 5 wires.The shaded quarter represents the calculation domain below.Calculation domain: area abcde, wire surface: line ab, plate: line cd, wire-to-plate gap: D, half of wire-to-wire distance: L. Voltage boundary conditions: Dirichlet for line ab and line cd, symmetry for line bc, line de and line ae; Charge boundary conditions: constant in the domain aa΄b΄b.

Fig. 13 .
Fig. 13.Relationships between the applied voltage and normalized electric field strength on the wire surface with respect to the wire radius variations.