Investigation of corona discharges in cylindrical geometries: application to wind turbines and aircraft

BACHELOR FINAL THESIS Investigation of corona discharges in cylindrical geometries. Application to wind turbines and aircraft Author: Asensio Garcia, Cristian Director /Co-director: Montanyà Puig, Joan / Sousa Arcanjo, Marcelo Augusto Degree: Bachelor's degree in Aerospace Technology Engineering Examination session: Spring, 2021 Document: Report UPC Lightning Research Group Investigation of corona discharges in cylindrical geometries. Application to wind turbines and aircraft Report Degree: Bachelor’s Degree in Aerospace Technology Engineering Course: Bachelor final thesis Delivery date: 22-06-21 Student: Asensio Garcia, Cristian Director: Montanyà Puig, Joan Co-Director: Sousa Arcanjo, Marcelo Augusto Abstract This thesis develops the university-level investigation of corona discharges in cylindrical geometries with application to wind turbines and aircraft. The main objective of this investigation is to create a code that simulates the conditions and development of electrical discharges, that studies the corona discharges in these conditions and that can be adapted to different cylindrical geometries. It is also taken into account the effect of wind and altitude, with the aim of ultimately simulating wind turbines blades and aircraft static dischargers under high electric fields. The methodology consist on performing an initial research on atmospheric electricity, aircraft electrostatic and corona discharges. Then creating a solver by first solving the electrostatics problem with open boundaries and secondly solving the continuity equation for time-dependent problems. The last step of the solver is to include the corona equations and the influence of wind and altitude. Once the solver is complete, simulations with different geometries are performed to simulate a vertical wire model that includes the corona effect, a static discharger model under the influence of the wind, and the evolution over time of positive ions next to the vertical wire. Regarding the wind turbines, the solver designed doesn’t simulate actual wind turbines but similar geometries that enhance the corona discharge effect. In this case, this similar geometry is the vertical wire which has been already used for this purpose by the UPC-LRG research group[1]. The report concludes with a discussion of the results and recommendations to implement this solver on more realistic geometries and complement the corona discharge effect with additional equations. Resum Aquesta tesis desenvolupa la investigació a nivell universitari de descàrregues de corona en geometries cilíndriques, amb aplicació a aerogeneradors i avions. L’objectiu principal d’aquesta investigació és crear un codi que simuli les condicions i el desenvolupament de les descàrregues elèctriques, que estudiï les descàrregues de corona en aquestes condicions i que es pugui adaptar a diferents geometries cilíndriques. També es té en compte l’efecte del vent i l’altitud, amb la finalitat de simular en última instància les pales dels aerogeneradors i els descarregadors estàtics d’aeronaus sota camps elèctrics elevats. La metodologia consisteix a realitzar una investigació inicial sobre electricitat atmosfèrica, electrostàtica d’avions i descarregues de corona. A continuació, es crea un solver resolent primer el problemes d’electrostàtica amb límits oberts i, en segon lloc, resolent l’equació de continuïtat per a problemes dependents del temps. L’últim pas del solver és incloure les equacions de la corona i la influència del vent i l’altitud. Un cop completat el solver, es realitzen simulacions amb diferents geometries per simular un model de fil vertical que inclou l’efecte corona, un model de descarregador estàtic sota la influència del vent i l’evolució al llarg del temps d’ions positius a prop del fil vertical. En quant als aerogeneradors, el solver creat en aquest treball no simula pròpiament aerogeneradors però si que simula geometries similars que realcin els efectes de la descàrrega de corona. En aquest cas, la geometria similar es el fil vertical que ja s’ha utilitzat prèviament amb aquest objectiu pel grup de recerca UPC-LRG [1]. La memoria conclou amb una discussió sobre els resultats i recomanacions per implementar aquest solver en geometries més realistes i complementar l’efecte de descàrrega de corona amb equacions addicionals. Declaration on Honour I declare that, the work in this Degree Thesis is completely my own work, no part of this Degree Thesis is taken from other people’s work without giving them credit, all references have been clearly cited, I’m authorised to make use of the research group UPC-LRG related information I’m providing in this document. I understand that an infringement of this declaration leaves me subject to the foreseen disciplinary actions by The Polytechnical University of Catalonia - BarcelonaTECH. Cristian Asensio García 22-06-21 Name of the Student Signature Date Thesis title: Investigation of corona discharges in cylindrical geometries. Application to wind turbines and aircraft Contents List of Figures 1 List of Tables 3 Acronyms 4 List of Symbols 4 1 Introduction 6 1.1 Objective ..................................... 6 1.2 Scope ....................................... 6 1.3 Requirements ................................... 7 1.4 Justification.................................... 7 2 Development 8 2.1 Background .................................... 8 2.2 Stateoftheart .................................. 8 2.3 Methodology ................................... 9 2.4 Approach and selection of alternatives . . . . . . . . . . . . . . . . . . . . . 10 3 Vertical wire model theory 12 3.1 Potential and Electric field . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.2 Coronainceptionzone .............................. 13 3.3 Continuityequation................................ 16 4 Solver 17 4.1 Poissonsolver................................... 17 4.1.1 Discretized equation . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.1.2 Openboundaries ............................. 19 4.1.3 Neumann boundary condition . . . . . . . . . . . . . . . . . . . . . . 19 4.2 Successive Over Relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.3 Continuity equation solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.3.1 Discretization............................... 20 4.3.2 Numerical stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.4 Iondriftsolver .................................. 22 5 Vertical wire model 24 5.1 Geometry ..................................... 24 5.2 Results....................................... 24 5.2.1 Corona inception zone . . . . . . . . . . . . . . . . . . . . . . . . . . 24 5.2.2 Ion trajectories with wind and fair weather . . . . . . . . . . . . . . . 27 5.2.3 Ion trajectories with wind and severe weather . . . . . . . . . . . . . 32 5.3 Similarcasegeometry............................... 35 5.4 Similarcaseresults ................................ 36 6 Static discharger model 37 6.1 Geometry ..................................... 37 6.2 Results....................................... 38 7 Ion drift model 42 7.1 Geometry ..................................... 42 7.2 Results....................................... 43 8 Discussion of results 47 8.1 Verticalwiremodel................................ 47 8.2 Static discharger model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 8.3 Iondriftmodel .................................. 48 8.4 Windturbinesmodel............................... 48 9 Conclusions 50 Bibliography 51 A Code 53 B Poisson solver validation 56 C Relaxation factor 58 Investigation of corona discharges in cylindrical geometries List of Figures 1 Work Breakdown Structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2 Coefficient values for different values of electric field. . . . . . . . . . . . . . 15 3 Effective ionization coefficient value for different values of electric field and altitude. ...................................... 16 4 Five point stencil representation. . . . . . . . . . . . . . . . . . . . . . . . . 17 5 SORscheme. ................................... 20 6 Iondriftsolver. .................................. 23 7 Corona inception altitude for different ambient electric fields (in logarithmic scale) ranging from 120 V/m to 5000 V/m, with two corona inception indicators. 25 8 Corona inception zone size for different ambient electric fields ranging from 120 V/m to 5000 V/m, for a 100 mhighwire. ................. 26 9 Corona inception zone size for different ambient electric fields ranging from 120 V/m to 5000 V/m, for a 200 mhighwire. ................. 27 10 Ion speed module for Eambient = 120 V/m, zoomed view of the vertical wire tip, with normalized vectors of 0.1mand a logarithmic color-map. . . . . . . 29 11 Ion speed module for Eambient = 120 V/m, zoomed view of the vertical wire tip, with normalized vectors and a logarithmic color-map. . . . . . . . . . . . 30 12 Ion speed module for Eambient = 120 V/m, zoomed view of the vertical wire tip, with normalized vectors of 0.08 mand a logarithmic color-map. . . . . . 31 13 Ion speed module for Eambient = 1204 V/m, zoomed view of the vertical wire tip, with normalized vectors of 0.1mand a logarithmic color-map. . . . . . . 32 14 Ion speed module for Eambient = 1204 V/m, zoomed view of the vertical wire tip, with normalized vectors and a logarithmic color-map. . . . . . . . . . . . 33 15 Ion speed module for Eambient = 1204 V/m, zoomed view of the vertical wire tip, with normalized vectors of 0.08 mand a logarithmic color-map. . . . . . 34 16 Real TARS geometry. Source: [19] ....................... 35 17 Corona inception zone size for different wire heights ranging from 200 mto 5000 m, for an ambient electric field of 120 V/m. ............... 36 18 Static discharger geometry and model. . . . . . . . . . . . . . . . . . . . . . 37 19 Ion speed module, for a static wick, with normalized vectors of 0.08 mand a logarithmiccolor-map. .............................. 39 20 Ion speed module, for a static wick, with normalized vectors of 0.08 mand a logarithmiccolor-map. .............................. 40 21 Ion speed module, for a static wick, with normalized vectors of 0.08 mand a logarithmiccolor-map. .............................. 41 1 Investigation of corona discharges in cylindrical geometries 22 Ion drift model domain representation: wire in black, final domain in blue and Gaussian distribution above the wire. . . . . . . . . . . . . . . . . . . . 42 23 2D distribution of positive ions density with interpolation to smooth the result: initial distribution at left and final distribution at right. . . . . . . . . . . . . 44 24 1D distribution of positive ions density of the vertical nodes on the zaxis, for differentpointsintime............................... 45 25 Speed verification of the ion drift model. Real speed computed numerically, and Mobility speed computed as |~ Wp|=µp~ E. ................. 46 26 Code scheme. On the left the functions used for the ion drift model, on the right the functions used for the vertical wire and on the bottom the ones used for the static discharger. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 27 Validations results extracted from Jesús Alberto López Trujillo [11] except for the ’New solver’ curve, which is the model developed in this thesis. . . . . . 57 28 Optimum relaxation factor found experimentally. . . . . . . . . . . . . . . . 58 2 Investigation of corona discharges in cylindrical geometries that in recent years has been heightened further due to the use of composite materials, that have more issues with dissipating the electric charges on their surfaces. It is also know that corona discharges are greatly affected by wind [7], but deeper studies need to be done to understand its behavior completely and control them. That is why research is being done by the UPC Lightning Research Group and the Department of AeroAstro at MIT as a means of predicting the corona discharges, with simulations in wind tunnels and measurements done by drones. In that context, this thesis aims to investigate the conditions and development of corona discharges from cylindrical geometries with the aim of determining the occurrence in tips of wind turbine blades and aircraft static dischargers. 2.3 Methodology In this thesis an academic-based approach is proposed to study the corona discharges. The first step in any investigation is to perform an initial research on the topic, in this case information about atmospheric electricity, aircraft electrostatic and corona discharges was gathered. The next step consists on creating a solver that would meet the requirements. Since lightning parameters are not easily measurable, predictive modeling [8] is often applied to model the pre-lightning strike and aircraft electrostatics, therefore the same methodology is followed in this project. The solver consists on a MATLAB program that requires the additional use of the FEMM software. The solver starts with solving the electrostatics problem with open boundaries and then solving the equations for time-dependent problems. The last step of the solver is to include the corona equations and the influence of wind and altitude. The process followed with the different steps and time dedication can be seen in the following Work Breakdown Structure in figure 1: 9 Investigation of corona discharges in cylindrical geometries Figure 1: Work Breakdown Structure. 2.4 Approach and selection of alternatives After an initial research on the topic, in some aspects of the project a decision on how to proceed was needed. Therefore in this section the different alternatives and the final approaches are discussed. 10 Investigation of corona discharges in cylindrical geometries Initially, when designing a solver the main decision to be made is programming language, in this case two options are considered: MATLAB and C++. •MATLAB: It is programming and numeric computing platform with a programming language that expresses matrix and array mathematics directly. Its main advantages are the simplicity of the language, the pre-programmed functions and its capabilities to post-process and plot data. •C++: It has object-oriented and generic programming features, and provides facilities for low level memory manipulation. Since it compiles directly to the machine’s native code, it is one of the fastest languages when it is optimized, and by far faster than MATLAB. Due to the lack of expertise on C++ and the broader experience on MATLAB, in addition to the great versatility of MATLAB to post-process data, it is finally decided to code the solver in MATLAB language. Another decision that will be discussed later, on section 3.1, is the method of computing the induced charge on a vertical wire due to the potential gradient of the atmosphere. Two alternatives are being considered: the Method of moments or the use of a commercial software such as FEMM. It is finally decided to compute the induced charge with the FEMM software due to its faster implementation. Regarding the different types of models simulated, it is decided to simulate a static discharger and a vertical wire, which is used as an approximated geometry of a wind turbine. The reason for this is that, a complete wind turbine with its blades cannot be simulated with the axisymmetric cylindrical geometry, because the blades break the axisymmetric geometry. Therefore a thin vertical wire is simulated with a radius equivalent to the radius of curvature of the trailing edge of a wind turbine blade, as previously done in Pol Fontanes Molina et al. [1]. The final decision to be made is whether to perform experimental tests or not. However, due to a limitation of time it is decided not to include any experimental result in this report, as it would be required to conduct many complex experimental tests to validate the simulation model. 11 Investigation of corona discharges in cylindrical geometries 3 Vertical wire model theory 3.1 Potential and Electric field Corona discharges depend on the values of the electric field at each point, so the first step in this model is obtain the potential field and electrical field for a vertical wire ideally suspended on air with fair weather conditions. In order to simulate this, a potential that increases with altitude can be estimated along the boundary of the domain as V(h) = |~ Eambient| · h, where the ~ Eambient is the ambient electric field with a value between 100 V/m and 120 V/m [9] in fair weather conditions. The vertical wire is connected to the ground or has a fixed potential so it remains with a V=ct at all times. Given these conditions, the electric potential is found for the whole domain by solving Laplace equation of electrostatics. ∇2V= 0 Where Vis the potential field at each point and ∇2is the Laplacian operator that for cylindrical coordinates with axisymmetric geometry, where V(θ) = ct, results in: ∂2V ∂r2+1 r ∂V ∂r +∂2V ∂z2= 0 Moreover, this potential field created by the atmosphere induces a negative surface charge density in the wire that is calculated with the FEMM software. This negative induced charge generates at the same time a potential and electric field that must be solved with Poisson’s equation. ∇2V=−ρ 0 Where ρis the space charge density of the wire in [C m−3]and 0is the free space permittivity in [F m−1]. Again developing the Laplacian operator, it is obtained: ∂2V ∂r2+1 r ∂V ∂r +∂2V ∂z2=−ρ 0 Finally the total potential of the domain is calculated with the superposition of the potential created by the electric ambient field, and the potential created by the surface charge on the wire, in other words: Vtotal =Vambient +Vcharge. Therefore the total electric field is: ~ Etotal =−∇Vtotal. 12 Investigation of corona discharges in cylindrical geometries 3.2 Corona inception zone The first step in this corona discharge study is to determine the conditions of inception of corona discharges. In other words, it must be calculated if there is a chance to obtain electrical breakdown at the tip of the wire, which would ignite a corona discharge. To do so, the Townsend criterion for spark formation or Townsend breakdown criterion is applied. In non-uniform fields like this one, and at low pressures the Townsend criterion for sparks is [10]: γñexp ÇZd 0 αdxå−1ô= 1 Where γare the number of electrons released from the cathode per incident positive ion, d is the set gap length, and αis the effective ionisation coefficient. The integration must be taken along the line dx of the highest field strength. Taking into account a non-uniform distribution of α, the criterion condition for breakdown results in: exp ÇZxc<d 0 αdxå=Ncr Where Ncr is the critical electron concentration in an avalanche that gives rise to the initiation of a streamer1which is shown to approx. 108m−3[10], and xcis the path of avalanche to reach the Ncr. This equation can be written with the following formulation, which is the one that will be applied in the code in order to detect the region of inception of corona in a cylindrical wire with axisymmetric geometry [7]: Zs0 a αds = lnNcr ≈18 −20 Where ais the radius of the wire and sis the distance from the wire center along the local electric field line. Thus, the integration is done along the electric field line as the avalanche starts at the distance where the field breakdowns (s=s0) and ends on the wire (s=a). Since it is a non-uniform field, the field strength and hence the effective ionisation coefficient vary across the domain. The effective ionisation coefficient is computed as α=α−η, where αis the ionisation coefficient in air and ηis the electron attachment coefficient in air. Below are the equations used for the coefficients, which are extracted from R Morrow and JJ Lowke. [12]. In these equations |E|is the vector module of the electric field and Nis the density of the background gas (in this case air) at standard sea-level conditions (ISA), which is defined as N= 2.45 ×1019 cm−3according to Lipeng Liu and Marley Becerra [13]. The equations for the coefficients are: 1Streamers electrical discharges at atmospheric temperature, in other words plasma channels, that are generated in regions where the electric field reaches the breakdown criterion, such as electrified clouds, that rise and propagate creating filaments or ramifications. [11] 13 Investigation of corona discharges in cylindrical geometries Ionization coefficient αin air: α N=2 ×10−16exp Ç−7.248 ×10−15 |~ E|/N åcm2 if |~ E|/N > 1.5×10−15 V cm2 = 6.619 ×10−17exp Ç−5.593 ×10−15 |~ E|/N åcm2 if |~ E|/N ≤1.5×10−15 V cm2 Electron attachment coefficient in air: η=η2+η3 η2 N=8.889 ×10−5Ç|~ E| Nå+ 2.567 ×10−19 cm2 if |~ E|/N > 1.5×10−15 V cm2 = 6.089 ×10−4Ç|~ E| Nå−2.893 ×10−19 cm2 if |~ E|/N ≤1.5×10−15 V cm2 η3 N2=4.7778 ×10−59 Ç|~ E| Nå−1.2749 cm2 (1) Since the coefficients depend only on the variable electric field, they can be plotted: 14 Investigation of corona discharges in cylindrical geometries Figure 2: Coefficient values for different values of electric field. It is of great interest to study the corona effect at different altitudes, therefore a correction for the αand ηcoefficients is introduced by modifying the number density Nof the air at each altitude. This correction follows the International Standard Atmosphere (ISA) model, where for h < 11000m[14] provides the following equations for finding the temperature T, the pressure P, the air density ρand the number density Nat each altitude. T= 15.04 −0.00648h[◦C] P= 101.29 ïT+ 273.1 288.08 ò5.256 [kPa] ρ=P 0.2869(T+ 273.1) [kg/m3] In order to use more useful units the density is converted to [g/cm3]. Then it can finally be obtained the number density of air: N=NA MM ρ Where NAis Avogadro’s number and MM the molar mass of air. A new plot can be obtained for the effective coefficient against the electric field value for different altitudes: 15 Investigation of corona discharges in cylindrical geometries Figure 3: Effective ionization coefficient value for different values of electric field and altitude. 3.3 Continuity equation When studying the behavior of electrons and positive ions on the air Poisson’s equation is coupled with the continuity equation that takes into account the movement of ions in air due to a electric field. Many streamer models contain continuity equations with multiple species and reactions, but in this investigation the continuity equation are simplified to the simplest case where only the positive ion mobility is taken into account. From the continuity equations provided by Lipeng Liu and Marley Becerra [15] and taking into account only the positive ions, the continuity equation ends up being: ∂Np ∂t =−∇ · (Npµp~ E) Applying cylindrical coordinates with axisymmetric geometry: ∂Np ∂t =−µp r d(r·NpEr) dr −µp d(NpEz) dz Where Npis the number density of positive ions particles, that is related with the space charge density ρ=e·Np, being ethe electron charge of 1.6022 ×10−19 C. The parameter µpis the ion mobility for positive ions, which relates the electric field with the drift velocity of ions ~ Wp. For positive ions in air, the ion mobility is µp= 2.34 cm2/V s [15] and constant in rand z. ~ Wp=µp~ E 16 Investigation of corona discharges in cylindrical geometries 4 Solver In order to study the effect of corona discharges on cylindrical geometries a solver has been created with the following characteristics: •2D solver in cylindrical coordinates with axisymmetric geometry, f(θ) = ct. •Node-centered discretization and a quadrilateral mesh. •Successive Over-Relaxation (SOR) method applied on the Poisson equation. •Open boundaries in order to reduce the size of the domain. Since the aim of this solver is to simulate real problems, such as wind turbines and static dischargers, the axisymmetric geometry is chosen because it is a geometry that allows to solve 3D simple geometries by reducing them to a 2D problem. In this section, the different parts of the solver are explained, with specific emphasis on the numerical methods applied. The solver can be found on the appendix A. 4.1 Poisson solver Recalling Poisson’s equation in cylindrical coordinates with axisymmetric geometry, where V(θ) = ct: ∇2V=−ρ 0 ∂2V ∂r2+1 r ∂V ∂r +∂2V ∂z2=−ρ 0 (2) A discretization of the domain in 2 dimensions is performed by using the five stencil points on a rectangular mesh, as can be seen on figure 4. Figure 4: Five point stencil representation. 17 Investigation of corona discharges in cylindrical geometries Therefore the second degree Taylor Expansion is applied to compute the derivatives, the following process is extracted and modified from Jesús Alberto López Trujillo [11]: V1=V0+∂V ∂z hz+1 2 ∂2V ∂z2h2 z V2=V0−∂V ∂z hz+1 2 ∂2V ∂z2h2 z V3=V0+∂V ∂r hr+1 2 ∂2V ∂r2h2 r V4=V0−∂V ∂r hr+1 2 ∂2V ∂r2h2 r                        (3) 4.1.1 Discretized equation By adding the equations of V1and V2from (3) it is obtained: V1+V2= 2V0+∂2V ∂z2h2 z−→ ∂2V ∂z2=1 h2 z (V1+V2−2V0)(4) Repeating the procedure with V3and V4from (3), it is obtained: ∂2V ∂r2=1 h2 r (V3+V4−2V0)(5) Subtracting V4from V3the remaining derivative can be obtained: ∂V ∂r =1 2hr (V3−V4)(6) By replacing the three derivatives on the Poisson equation (2) it can be found the standard discretized equation for the new value of V0as a function of its adjacent points: V0=1 4rh2 z+ 4rh2 rï2rh2 rV1+V2)+2rh2 z(V3+V4+hrh2 z(V3−V4) + 2rh2 rh2 zρ 0ò For the nodes on the axis, it shall be used the Poisson equation taking into account the singularity on r= 0, where L’Hôpital’s rule is applied: lim r→0 1 r ∂V ∂r =∂2V ∂r2: 2∂2V ∂r2+∂2V ∂z2=−ρ 0 By replacing the derivatives on this modified Poisson equation and knowing that for the nodes on the axis V3=V4, it can be obtained: V0=1 4rh2 z+ 4rh2 rïh2 rV1+V2)+2h3 z(V3+V4+h2 rh2 zρ 0ò This discretized equation can be also applied to solve the Laplace equation imposing that ρ= 0. 18 Investigation of corona discharges in cylindrical geometries one cell (0.75 mm). This initial radius is not realistic because it can happen that the radius is smaller than one cell, but due to discretization limitations it can only be represented as one cell. Therefore other indicators for quantifying when the corona inception zone appears must be adopted. The thresholds are arbitrarily defined as the moment when the base of the corona inception zone is equal or greater than the radius of three cells and the radius of four cells. In this case, three cells correspond to a radius of 2.3mm and four cells to 3mm. The following graph 7 shows the height of the wire where the corona zone is initiated for different conditions of ambient electric field (in logarithmic scale). The electrical field ranges from conditions of fair weather of 120 V/m up to severe weather conditions of 5000 V/m. The corona inception altitude is represented for the two indicators previously explained: the 2.3mm and 3mm corona radius. Figure 7: Corona inception altitude for different ambient electric fields (in logarithmic scale) ranging from 120 V/m to 5000 V/m, with two corona inception indicators. As it can be seen, the first indicator (2.4mm) and the second indicator (3mm) follow a stable decrease trend, that suggest that as the ambient electric field increases the size of the corona increases too in the vertical direction. It should be mentioned, that after many simulations it has been observed that the corona model has a slight irregularity for small values of electric field. That’s because due to how 25 Investigation of corona discharges in cylindrical geometries the breakdown criterion is implemented,there is always a residual corona region for small values of electric field. For this reason, the first values for the 2.4mm indicator are not represented. This previous result can be contrasted with simulations that study how the size of the corona changes with the ambient electric field. Since the corona doesn’t have a simple and constant structure, the size of the corona inception zone will be defined by its maximum radius and its 3D volume, which will be computed numerically. The following graph represents the corona radius and volume vs different ambient electrical fields ranging from 120 V/m to 5000 V/m, for a fixed wire height. Figure 8 shows the results of a 100 mwire, while figure 9 shows the results from a 200 mwire. Figure 8: Corona inception zone size for different ambient electric fields ranging from 120 V/m to 5000 V/m, for a 100 mhigh wire. 26 Investigation of corona discharges in cylindrical geometries Figure 9: Corona inception zone size for different ambient electric fields ranging from 120 V/m to 5000 V/m, for a 200 mhigh wire. It can be observed that for both cases of 100 mand 200 mof wire height, the graphs follow the same pattern: the corona radius grows with the electric field, and the corona volume decreases initially and then gradually increases with the electric field. However, the difference between wire heights lies on the values of both the corona radius and corona volume. For the wire of 100 mof height the maximum corona radius increases at a rate of 1.8mm every 1000 V/m, while for the wire of 200 mof height the maximum corona radius increases at a rate of 2.3mm every 1000 V/m. And for severe weather conditions of 5000 V/m the maximum corona radius is 50% higher for the 200 mhigh wire compared to the 100 mhigh wire. Moreover, regarding the 3D corona volume the rate of increase is 4 times higher for the 200 mhigh wire. Again it should be mentioned that, the initial decrease of the 3D corona volume may be as a result of the irregularity of the model for small values of electric field. 5.2.2 Ion trajectories with wind and fair weather Another important effect that can be simulated is the effect of wind to the ion trajectories of the ions around the vertical wire. In this case, one of the geometries simulated in the previous section 5.1, the one with a 200 mhigh vertical wire and 120 V/m of ambient electric 27 Investigation of corona discharges in cylindrical geometries field, is simulated with an additional wind in the radial positive direction ~ Vwind =Vwind ·ˆur. The geometry can be summarized with the parameters of the following table 3. Table 3: Geometry parameters for the ion trajectory simulation. Wire height 200 m Nr: cells in r 1000 Wire radius 0.75 mm Nz: cells in z 1200 Domain radius 0.75 m hr: cell size in r 0.75 mm Domain height 600 m hz: cell size in z 0.5 m Eambient: ambient electric field 120 V/m In the following figure 10, the positive ion trajectories are plotted along with the ion speed module, in a zoomed view of the vertical wire tip, when no wind is applied. It can be seen that above the wire, when the ions are not disturbed by the presence of the wire, the trajectories are vertical as they follow the direction of the ambient electric field. However, as the ions get closer to the wire they get attracted to it due to the electric field it generates. 28 Investigation of corona discharges in cylindrical geometries Figure 10: Ion speed module for Eambient = 120 V/m, zoomed view of the vertical wire tip, with normalized vectors of 0.1mand a logarithmic color-map. In the following figures 11 and 12, the trajectories are plotted for different wind speeds from 10 m/s up to 50 m/s. It can be seen that as the wind appears it drags with himself most of the ions, and only the ions closer to the wire remain attracted. As the wind speed increases, less ions are attracted to the wire up until 50 m/s where the ions attracted are the ones from the cells next to the wire. 29 Investigation of corona discharges in cylindrical geometries (a) Zoomed view with 10 m/s wind and normalized vectors of 0.1m. (b) More zoomed view with 10 m/s wind and normalized vectors of 0.08 m. Figure 11: Ion speed module for Eambient = 120 V/m, zoomed view of the vertical wire tip, with normalized vectors and a logarithmic color-map. 30 Investigation of corona discharges in cylindrical geometries (a) Zoomed view with Wind of 30 m/s. (b) More zoomed view with Wind of 30 m/s. (c) Zoomed view with wind of 50 m/s. (d) More zoomed view with wind of 50 m/s. Figure 12: Ion speed module for Eambient = 120 V/m, zoomed view of the vertical wire tip, with normalized vectors of 0.08 mand a logarithmic color-map. 31 Investigation of corona discharges in cylindrical geometries 5.2.3 Ion trajectories with wind and severe weather In this subsection, the previous geometry of the 200 mwire, is simulated for conditions of approaching severe weather, which is approximated with an ambient electric field of 1204 V/m. In the following figure 13, the positive ion trajectories are plotted along with the ion speed module, in a zoomed view of the vertical wire tip, when no wind is applied. It can be seen that the results are similar to the ones obtained with fair weather, in section 5.2.2, but in this case the trajectories are even more influenced by the electric field of the wire since they get attracted with higher speeds. This occurs because a higher electric ambient field induces a bigger charge in the wire and hence, a bigger electric field is generated. Figure 13: Ion speed module for Eambient = 1204 V/m, zoomed view of the vertical wire tip, with normalized vectors of 0.1mand a logarithmic color-map. In the following figures 14 and 15, the trajectories are plotted for different wind speeds from 10 m/s up to 50 m/s. It can be seen that the graphs are nearly identical as the case of fair 32 Investigation of corona discharges in cylindrical geometries weather, section 5.2.2, but in this case more ions are attracted to the wire since more arrows are pointing towards it. Moreover, higher speeds are achieved in the nodes closer to the wire due to the higher electric field. As the wind speed increases, again less ions are attracted to the wire. (a) Zoomed view with 10 m/s wind and normalized vectors of 0.1m. (b) More zoomed view with 10 m/s wind and normalized vectors of 0.08 m. Figure 14: Ion speed module for Eambient = 1204 V/m, zoomed view of the vertical wire tip, with normalized vectors and a logarithmic color-map. 33 Investigation of corona discharges in cylindrical geometries (a) Zoomed view with Wind of 30 m/s. (b) More zoomed view with Wind of 30 m/s. (c) Zoomed view with wind of 50 m/s. (d) More zoomed view with wind of 50 m/s. Figure 15: Ion speed module for Eambient = 1204 V/m, zoomed view of the vertical wire tip, with normalized vectors of 0.08 mand a logarithmic color-map. 34 Investigation of corona discharges in cylindrical geometries (a) Wind of 100 m/s. (b) Wind of 150 m/s. Figure 21: Ion speed module, for a static wick, with normalized vectors of 0.08 mand a logarithmic color-map. 41 Investigation of corona discharges in cylindrical geometries 7 Ion drift model The final model to be implemented is the ion drift model of a 2D Gaussian distribution of positive ions under the influence of an electric field generated by a vertical wire. 7.1 Geometry In this model the geometry used is one of the geometries already used in the vertical wire model from section 5.1. A thin wire of 0.75 mm of radius and a height of 100 mand an ambient electric field of 120 V/m, where the domain is only the region of the tip of the wire as can be seen in figure 22. The geometry can be summarized with the parameters of the following table 6 Table 6: Geometry parameters for the ion drift model. Wire height 100 m Nr: cells in r 1000 Wire radius 0.75 mm Nz: cells in z 2500 Domain radius 0.5 m hr: cell size in r 0.5 mm Domain height 1 m hz: cell size in z 0.4 mm Eambient: ambient electric field 120 V/m Figure 22: Ion drift model domain representation: wire in black, final domain in blue and Gaussian distribution above the wire. 42 Investigation of corona discharges in cylindrical geometries The seed population of positive ions Npfollows a 2D Gaussian distribution that is placed on the zaxis a little above the wire and that is only half Gaussian due to the symmetry, as can be seen in the previous figure 22. The equation that follows the 2D Gaussian is: Np=N0·exp ï−Å(x−µx)2+ (y−µy)2 2σ2ãò With values taken from TN Tran et al. [22] and adapted to this case: Table 7: Parameters for the 2D Gaussian distribution. N0: Gaussian peak 1×107cm−3 σ: deviation 0.25 mm µr: offset in r 0 µz: offset in z 2.5mm above the wire 7.2 Results Due to the long duration of ion drift model simulations, only one simulation is performed simulation a total period of time of 161 µs, and a time step of dt = 1 µs. This time step is calculated in each time iteration to ensure that it follows the CFL condition (previously mentioned), but since the conditions of the simulation don’t vary significantly between each time step, the time step is rounded in each iteration to the same value. Table 8: Ion drift simulation, time parameters. dt: Time step: 1 µs Number of time steps 161 Simulation time 159 µs The first result that can be obtained is a 2D representation of the movement of positive ions by comparing the initial seed Gaussian and the final distribution, as can be seen in figure 23. It can be seen that the distribution of positive ions moves downwards towards the wire, the wire ends at 100 m. The distribution seems to keep a sort of Gaussian shape, but deformed in the vertical direction, as if physically the Gaussian was pulled down. This effect is confirmed by the maximum value of the distribution, that seems to decrease slightly. 43 Investigation of corona discharges in cylindrical geometries Figure 23: 2D distribution of positive ions density with interpolation to smooth the result: initial distribution at left and final distribution at right. 44 Investigation of corona discharges in cylindrical geometries A better representation of the movement can be seen in figure 24, where the 1D representation of the vertical nodes on the zaxis is plotted. It is observed that the distribution of particles flattens out and skews to the left as time goes by, and as a result the peak value of the distribution is reduced. Due to computation limitations, this simulation can only be performed for short periods of time, but the results suggest that for higher times of simulations the distribution would skew to left up until it reached the wire. Figure 24: 1D distribution of positive ions density of the vertical nodes on the zaxis, for different points in time. 45 Investigation of corona discharges in cylindrical geometries Finally, the speed of movement of the Gaussian can be computed numerically and can be compared with the drift velocity of ions |~ Wp|=µp~ Edue to the mobility µp. This comparison is shown in figure 25. Where the speed of the Gaussian is computed by finding the centroid of the 2D Gaussian at each time and calculating the difference in position with respect to the previous step. It can be seen that the speed of drift due to the electric field matches the speed of movement of the Gaussian. Moreover, both speeds increase as time goes by, this is because as the time increases the Gaussian gets closer to the wire tip, where the electric field is higher and thus it becomes more attracted to the wire. Figure 25: Speed verification of the ion drift model. Real speed computed numerically, and Mobility speed computed as |~ Wp|=µp~ E. 46 Investigation of corona discharges in cylindrical geometries 8 Discussion of results 8.1 Vertical wire model From the results of the corona inception size obtained in section 5.2.1 it can be extracted that: •The corona inception altitude decreases as the ambient electric field increases, if the irregularity of the model for small electric values is prevented. That is because a higher electric field induces more charge on the wire, which at the same time generates a higher electric field near the wire. •The size of the corona increases as the ambient electric field increases, since both the corona radius and the 3D corona volume increase linearly with the ambient electric field. As explained previously, the ratio of increase of the corona radius is around 20% higher for the wire of 200 mof height compared with the 100 mhigh wire, and for severe weather conditions of 5000 V/m the corona radius is bigger for the high wire. Which means that for higher wires the corona effect becomes more important as the weather conditions become more severe. Regarding the 3D corona volume, the results show an irregularity for fair weather conditions of 120 V/m due to the fact that the model is not accurate for small values of electric field. That is because the way the breakdown criterion is implemented, it generates always a residual corona corona for small values of electric field. It can be seen that the 3D corona volume increases 4 times faster for the wire of 200 mof height, which again verifies that for higher wires the corona effect is heightened. From the results of the ion trajectories in a vertical wire from section 5.2.2 with and without wind in the radial positive direction ~ Vwind =Vwind ·ˆur. It can be seen that: •When no wind is applied the ions follow the trajectory of the ambient electric field, but as they closer to the wire they get attracted to it and bend their trajectory until they are perpendicular to the conductor. As the ambient electric field is increased, the ions are more strongly attracted to wire, due to a bigger charge density on the wire. •When wind is applied, it drags with himself most of the ions of the domain, and only the ions that are really close to the wire remain attracted. As the wind increases, this effect is heightened up to a point where only the ions on the cells next to the wire are attracted. The effect of increasing the ambient electric field is that the ions are more attracted to the wire even with high wind speeds. 47 Investigation of corona discharges in cylindrical geometries Regarding the results from the similar case geometry, the TARS, in section 5.4, it can be stated that: As the wire height increases, both the maximum corona radius and 3D corona volume increase linearly. 8.2 Static discharger model From the results of the influence of wind on the vertical positive direction ~ Vwind =Vwind ·ˆuz on a static discharger, obtained in section 6.2, it can be said that: •When no wind is applied, all the ions are attracted to the static discharger as they "bend" their trajectory to follow the direction of the electric field generated by the conductor. In other words, a direction perpendicular to the surface of the static discharger. •When wind is applied, the ions start to follow the trajectory of the wind and they feel less attracted to the static discharger, and as wind increases this behavior increases too. For high winds of 150 m/s or more, the ions are little attracted to the conductor, and they fly out with the wind. With all of this results, it can be said that the aircraft static discharger model already mimics the real shape of a static wick. Future improvement of this model would be the proposal of new electrode shapes to increase the wind speed. 8.3 Ion drift model From the results of the simulation of the ion drift on the tip of a vertical wire in section 7.2, it can be said that positive ions particles are attracted to the charged conductors as it has already been proven with the other two previous models. From a 1D point of view, the Gaussian distribution skews towards the electrode tip, this behavior matches the results of previous 1D models such as Floris Wulf [23]. From the comparison of the Gaussian speed of movement in figure 25 it can be observed that the speed of drift due to the electric field matches the simulated Gaussian speed of movement. This result, can be considered as a first validation of the ion drift model. 8.4 Wind turbines model Regarding the implementation of this models on wind turbines, the vertical wire is a preliminary step of simulating wind turbine blades. The similarity between these two models lays on the fact that the radius of the wire is equivalent to the radius of curvature of the trailing edge of a wind turbine blade. The trailing edge of a wind turbine blade is the sharpest point of a wind turbine, thus the place where a corona discharge could appear, therefore by only 48 Investigation of corona discharges in cylindrical geometries simulating the trailing edge, the parts more exposed to discharges are already included. As a consequence, the model can be reduced to a simple 2D cylindrical axisymmetric geometry. 49 Investigation of corona discharges in cylindrical geometries 9 Conclusions The goal of this thesis was to investigate and simulate the corona discharges in cylindrical geometries, having the ultimate intention of simulating real applications. Based on the previous results, it can be concluded that despite the difficulty of predictive models and corona discharge behavior, the objectives of the thesis are fulfilled. However, the simulation model and the thesis are still far from perfect, and conclusions can be established regarding the model and the solver. The code works properly and is successfully validated for the range of potential values needed for the corona discharge simulations. Moreover, despite not using the ideal programming language for this types of solver, which should be C++, the solver seems to converge for all the tests that have been made. But, on the other hand, the code is relatively slow, since some simulations require up to 16 hours. Therefore, for future improvements the solver could be optimized by implementing the multigrid method [24] or solving the main loops on C++. Finally, it can be said that the model would have been more developed and completed if validation tests had been conducted, but due to time limitations this was not able. 50 Investigation of corona discharges in cylindrical geometries Figure 27: Validations results extracted from Jesús Alberto López Trujillo [11] except for the ’New solver’ curve, which is the model developed in this thesis. 57 Investigation of corona discharges in cylindrical geometries C Relaxation factor The equation used for the optimum relaxation factor ωas a function of the number of cells Nis a logarithmic regression of the following type : ω=a·ln(N)4+b·ln(N)3+c·ln(N)2+d·ln(N) + e Table 9: Parameters for the ωequation. a b c d e -2.3614×10−50.001462 -0.03406 0.3549 0.5967 Figure 28: Optimum relaxation factor found experimentally. 58