Calculating for surface electric field of converter valve shield system with fast multipole curved boundary element method

: The accurate and fast calculation of surface electric field of DC converter valve shield system have guiding significance in the designing process. However, the electrostatic field analysis is a large-scale problem with low efficiency using the conventional numerical algorithm. In order to solve this problem, fast multipole curved boundary element method (FMCBEM) was proposed. The FMCBEM is based on the indirect boundary integral equation of electrostatic field. Galerkin curved boundary element method based on coordinate transformation is used to improve the accuracy, then fast multipole method is used to improve the calculation speed and reduce the memory usage. Through simple examples, the accuracy and efficiency of the algorithm were verified. Compared with finite element method, it was found that the accuracy can meet the requirements of engineering. The efficiency was obviously better than the conventional boundary element method. Applying this algorithm, the surface electric fields of different shield systems of ±160 kV converter valve, which the degree of freedoms is more than 220 thousand, were analysed using desktop computer. This algorithm can calculate the large-scale electric field problem fast and accurately, which is an effective tool for the optimisation design.


Introduction
As the demand for the distributed generation to access the power grid increases, and the demand for interconnection of AC power grids increases, HVDC-flexible has attracted wide attention with the advantages that low requirements on the connected AC power grid and the reactive power and active power can be flexibly controlled [1]. Due to the converter valve which assembled in the limited space contain a large number of power electronic devices, the electromagnetic environment is complicated. The design of shield system is one of the key technologies for electromagnetic compatibility. Therefore, the shield system of converter valve must be designed reasonably to control the surface electric field, so as to ensure the safety margin between devices, and ensure no corona occupancy [2].
There is no unified standard for the shield structure design. Due to the structure of valve tower is complicated, the prototype test is much costly. Therefore, numerical simulation about the electric field could provide valuable references in the designing process [3]. The numerical analysis of shield system electric field is largescaled and complicated, the main methods are finite element method (FEM) and boundary element method (BEM). The FEM need to generate mesh in the whole computing domain, the degree of freedoms (DOFs) is large, so the requirement of the computer hardware is higher [4]. Compared with FEM, BEM reduces the DOFs and the spatial dimension by discretising computational domain boundary only. Besides, due to the analytic basic solution is used in BIE, the BEM has higher accuracy [5,6].
Since in general, the BIE is non-local, the matrix formed by BEM is fully populated. The numerical solution requires large amounts of time and memory. Thus, the application of BEM for large-scale problem is limited [7]. To overcome this obstruction, different fast algorithms have been proposed. The fast multipole method (FMM) introduces multipole expansion and local expansion using harmonic function [8], uses recursive operation in the hierarchical tree structure to complete the aggregation, transfer, and dispersion between source points and field points, so as to accelerate the computation of matrix vector product. Combining FMM and BEM, the FMBEM which the computational complexity is nearly linear with DOFs is formed [9].
The method of improving the computation accuracy in FMBEM is usually increasing the amount of element or using high-order element. The two methods will increase DOFs and not suitable for large-scale problems. Here, the curved boundary element method (CBEM) based on coordinate transformation was used, which can improve the calculation accuracy without increasing DOFs [10]. Combining with FMM, the fast multipole curved boundary element method (FMCBEM) was proposed. Applying this algorithm, the surface electric fields of different ±160 kV converter valve single bridge models were analysed using desktop computer.

Two typical models of ±160 kV valve tower
In the HVDC-flexible transmission system, the converter system contains three phase units which are composed of six bridges. The single bridge is composed of converter valve and reactor. The converter valve is composed of several valve sections, which consists of IGBT, diode, capacitor, and so on, connected in series electrically. Mechanically, the valve tier which composed several valve sections connected in series are interconnected with aluminium busbars. The coolant is distributed between valve sections by insulating manifold pipes [3]. In order to control the electric field, the shield system is necessary. In the valve hall of HVDC-flexible, the valve tower is supported on the ground by supporting insulators. Fig. 1 shows the calculation models of two typical valve tower used in the ±160 kV HVDC-flexible converter valve. The models are refined. Here, aiming at electric field distribution of shield system, the internal structure was less focused on. So the valve section with complex inner structure was simplified to cuboids, and the sheds on the insulators were also ignored. As shown in Fig. 1, the shield systems of the two models are different. In type I, the shield structure is set for two valve towers. Six corona rings are used for one level. In type II, the shield structure is set for one valve tower only. Every level contains eight corona rings, besides corona rings are also set near insulator and four corona spheres are set at the top of valve tower.

Voltage loading
The AC and DC mix in the converter valve which the highest harmonic frequency is reached at 4 KHz [4], so the induction electric field produced by the magnetic field changes can be ignored. Therefore, the electric field can be analysed according to the method of quasi-static electric field, as (1).
Two of them represent conduction current and displacement current. γ and ε is conductivity and permittivity, and ω is angular frequency. In order to simplify the calculation, the two items are compared. If the displacement current is far greater than the conduction current, the effect of conduction current can be ignored, therefore, it can be analysed based on the principle of electrostatic field. On the contrary, the influence of the displacement current can be ignored, and it can be solved based on the principle of current field. Here, the electric field in air and insulator is the attention. According to the material properties, the electric field is analysed by the method of electrostatic field. Electric fields of two single bridge models were analysed. The connection modes of the two single bridge models are shown in Fig. 2. The model of type I contains eight valve towers, and type II contains seven towers. The voltage distribution of the valve tower is the important input condition for the analysis of electrostatic field. The two simulation models of the corresponding converter systems are built. The voltage can be obtained. Here, the electric field is analysed when the AC input is −160 kV. The input is AC voltage and the output is DC voltage after the rectifier bridge. The voltage at a certain time is AC −160 kV at the AC input terminal, and the voltage is gradually increased, and the voltage at DC output terminal is 160 kV. The electric field at this moment can be analysed by the electrostatic field. The voltage of corona ring and bracket is same at the same position, and the voltage of module is the average of the left and right voltage.

Curved boundary element method
Let bounded domain Ω ⊂ ℝ 3 is the computational domain whose piecewise smooth surface Γ := ∂Ω is globally Lipschitz continuous. Γ can be subdivided into M smooth surfaces Γ = ⋃ i = 1 M Γ i . If the intersection of two different surfaces Γ i ∩ Γ i′ is non-void set, the intersection consists at most of a common vertex or a common edge. Let bounded parametric plane κ i ⊂ ℝ 2 , and the Then, there exists a smooth diffeomorphism γ i for each surface Generate smooth curved element for each surface Г i , the element amount is N i . Then, we i also consists at most of a common vertex or a common edge. Using the same diffeomorphism, we can obtain κ i = ⋃ j = 1 N i P j i , where P j i is 2D planar element. The diffeomorphism between curved element E and planar element P can be obtained by coordinate transformation.
For the typical structure in the power system, three different coordinate transformation methods are used, respectively, spherical, toroidal, and cylindrical transformation. In the Cartesian coordinate system, the coordinate of a point on surface satisfies the constraint equation, only two of (x,y,z) are independent. Through the corresponding coordinate transformation, one variable is constant and two independent variables can be found to form (a,b) plane parametric space [10]. The transformation process is shown in Fig. 3. The diffeomorphism γ 1 and γ 2 is achieved by spherical and cylindrical transformation. Through these transformations, the Cartesian coordinate of node a, which located at the common edge Γ 1 ∩ Γ 2 , can be transformed into spherical and cylindrical coordinates, respectively. The infinitesimal dS must be changed by Jacobian after transformation, where J T is Jacobian of transformation and is different in different coordinate transformation. The integral on the 3D surface is transformed into the integral on the parametric plane by the coordinate transformation. The process of generating curved element on 3D surface is difficult, however the plane element in surface which used in conventional BEM can be used to complete the calculation of CBEM. The mesh used in CBEM is shown in Fig. 4. Generate the plane mesh,   through the spherical transformation and the node coordinates of plane mesh, the planar mesh corresponding to curved mesh can be obtained. The only requirement is that the nodes are on the surface and the coordinates are known in advance. As the nodes of the curved mesh and the plane mesh are both on the surface of the model, and the coordinates are same. We do not need to know the shape of the actual curved element, only through the node coordinates can complete the transformation between curved mesh and planar mesh. The integral quantities, such as the coordinates of Gaussian points and the Jacobian are calculated in planar mesh, thus no discretisation error is introduced. Different coordinate transformations require different coordinate centres. When there are different surfaces in the model, it is necessary to transform the coordinate system to a specific position by translation and rotation [11]. Then, for complex models, the surfaces are transformed into a 2D parametric planes by different coordinate transformations. The electric field in the converter valve is excited by conductors under voltage and ceramic or composite insulators. In the electric field model presented here, these parts are taken into account. The indirect BEM is used to analyse this problem. The surface of metal element is covered with free surface charges, and the surface of dielectric is covered by surface charges representing dielectric polarisation [5]. At the boundary, let the field point be (x,y,z), the source point is (x′,y′,z′), and R is the distance between field point and source point. So the potential can be expressed as (3).
The potential of dielectric is unknown. According to the interface condition, the normal components of the electric displacement vector on both sides of interface are continuous, then (4) can be obtained.
where ε 0 is the vacuum permittivity. The unknown total surface charge density included free charge and polarisation charge is σ, and the known voltage is u. Let σ/ε 0 be the unknown variable, denoted by λ. According to the boundary condition on the conductor, the electric flux density on the conductor surface has only the normal component and is equal to the surface charge density, so the solution λ is the magnitude of the electric field strength on the conductor. For (4), ε a and ε b represent the relative dielectric constant of medium a and b, respectively. The normal direction is perpendicular to the interface and enters medium a from medium b. Generating mesh in 3D model, transforming the nodes into 2D parametric space (a,b), using linear interpolation of planar element, then applying Galerkin weighted residual method to (3) and (4), the CBEM equation can be obtained. The coefficient matrix formed by CBEM is fully populated, and the computation complexity is O N nod 2 . The complexity to solve using direct solver is O N nod 3 . The characteristic limits the solution scale of the conventional CBEM. However, the FMM can reduce the computation complexity near to O N nod .

Fast multipole curved boundary element method
The central idea of FMM is to expand 1/R by harmonic function [8]. As shown in Fig. 5, the spherical coordinate of p and q is (r 1 , θ 1 , ϕ 1 ) and (r 2 , θ 2 , ϕ 2 ), respectively, and satisfy r 1 > r 2 . Then, the expansion is: where o is the centre of expansion, and any point can be selected as the centre. S l,k and R l,k are the solid harmonics defined as P l k is the associated Legendre function and a superposed bar indicates the complex conjugate. Through this expansion, the coordinates of p and q in 1/R are separated, thus, the double layer integral in Galerkin CBEM can be separated into two independent single layer integrals. Thus, the multipole expansion and local expansion are introduced, and are recursive operated in the oc-tree structure [12] to complete the aggregation, transfer, and dispersion. The sketch map of FMBEM is shown in Fig. 6. The FMM can be used when the distance between source point and field point is far. First, Q2M represents aggregation, which is conducted to calculate the multipole expansion from source points in each leaf cell. Then, M2M represents the multipole expansion translation from a child cell to a parent cell until reaching level 2. The next step is the downward procedure, which includes M2L and L2L. M2L represents the transfer from multipole expansion to local expansion between interaction cells. L2L represents local expansion translation from a parent cell to a child cell, which starts from level 3 and ends once the highest level is reached. Finally, L2u represents dispersion, which is performed to calculate the voltage based on the local expansion in each leaf cell. The specific calculation formulas for FMM can be found in reference [7]. When source point and field point are close, the error becomes larger by FMM. Therefore, these interactions are calculated through (3) and (4) by direct numerical integration as the conventional CBEM.

Illustrative numerical results
The element used in Galerkin FMCBEM can be either triangular or quadrilateral. The method has been implemented in a code written in Fortran language and solved by single-core. All the computations were carried out on the same desktop computer with 3.3 GHz CPU and 8 GB memory. The computational accuracy of the algorithm was examined by comparing with FEM. The calculation of FEM was performed via COMSOL 5.2. Normalised root mean square error (RMSe) defined according to the reference [13]:  mesh, (b) planar mesh on 2D space (a, b), (c) curved mesh where the subscript i represents the global indices of the node, and the superscripts b and f denote the Galerkin FMCBEM result and the FEM result, respectively, x f max which is the normalised coefficient, represents the maximum value of the FEM result.
Here, the restarted generalised minimum residual method (GMRES(m)) was selected as the solver [14]. The restart parameter m was chosen as 8. The iteration stop conditions are the following three: (i) the relative residual norm of stopping, which is chosen to be 10 −6 , (ii) the maximum restart times, which is chosen to be 50, (iii) the relative residual norm drop rate is defined as the percentage difference between the relative residual norm for two successive restarts, which is chosen to be 4%. As long as one of the conditions is satisfied, the iteration is stopped. In addition, the sparse approximate inverse (SAI) preconditioner was used to accelerate the convergence of GMRES(m) [15]. In order to illustrate the validity of the algorithm, the following parameters are defined. t near is the time to calculate A near which is calculated by direct numerical integration, t sum is the total running time, and F m is the memory usage when the algorithm run.
A model of multiple coordinate transformations was used to verify this algorithm. The shield structure model of single valve tower of type II was used, which includes toroidal surface, cylindrical surface, and spherical surface. The valve tower is a three-layer structure including a lot of corona rings. The amount of element and node is 27,496 and 25,872 respectively, and the model and voltage are shown in Fig. 7.
The electric field strength of the nodes along line 1 is shown in Fig. 8. The result of FMCBEM is basically the same as 3D FEM. Artificial spherical boundary, which the voltage is zero, was set at 20 times of the model height in 3D FEM. The tetrahedron element amount 8.5 million and triangular element on the surfaces is 0.5 million.
The calculation time and efficiency between CBEM and FMCBEM was compared, as shown in Table 1. The CBEM is solved by Gaussian elimination method. The RMSes of the two methods is same. The calculation time and memory of FMCBEM are the 1.3% and 11.2% of CBEM, respectively. In the other word, FMCBEM can solve more large-scale problem quickly in the same hardware condition.

Application results
After calculation, the electric field of the two single bridges can be obtained. The surface electric fields of shield systems are shown in Figs. 9 and 10. The value represents the magnitude of the electric field strength. As the electric field strength has only the normal direction on the conductor surface, the positive value indicates that the direction is perpendicular to the surface outward, and the negative value indicates that the direction is perpendicular to the surface inward.
For type I, the amount of element and node is 111 thousand and 112 thousand, respectively. The running time is 7.1 h, and the memory is 1.3 GB. The maximum values of the surface electric field are 14.71 kV/cm, which are located at the position of max1 in Fig. 9, are much lower than 20 kV/cm, the average breakdown electric field strength of air [3]. The maximum position is located at the corner of the support fitting.
For type II, the amount of element and node is 223 thousand and 223 thousand, respectively. The running time is 26.9 h, and the memory is 5.2 GB. In conventional CBEM, if we use the single floating-point real number to store the coefficient matrix, it needs 185.3 GB, which is beyond the storage capacity of the desktop computer. The maximum values of the surface electric field are −8.47 kV/cm, which are located at the position of max2 in Fig. 10, are much lower than 20 kV/cm, the average breakdown electric  field strength of air. The maximum position is located at the corona ring near the support fitting.
Comparing the surface electric field of the two type, the corona ring near support fittings can reduce maximum values of the surface electric field. The maximum electric field strength on the shield structures of Type II was about 42% lower than those on Type I. The shield effect of type II is better.

Conclusion
The FMCBEM has been implemented for the electric field problem of converter valve shielded system. In this method, the conventional plane mesh is still used, and then the 2D parametric planar mesh corresponding to the surface curved mesh can be obtained by coordinate transformation. The boundary integration is constructed in the parametric space which the integral quantities, such as the coordinates of Gaussian points and the Jacobian, are calculated in planar mesh. Thus, no discretisation error is introduced. The FMM method is used to speed up the solution and enlarge the solution scale. Although FMCBEM is used to solve the electric field problem of converter valve shielded system, this method is also suitable for other large-scale electrostatic field problems.
The accuracy of FMCBEM can meet the requirements of engineering. The calculation speed and memory usage are much less than conventional CBEM. Besides, the electric field of two type of converter valve single bridge has been successfully analysed by this method. The maximum electric field strengths on the shield structures with Type I and Type II were, respectively, 14.71 and 8.47 kV/cm, all lower than the critical corona inception electric field strength of 20 kV/cm. The maximum electric field strength on the shield structures of Type II was about 42% lower than those on Type I.

Acknowledgments
This work was supported by the Fundamental Research Funds for the Central Universities (Grant No.2017XS006).