Calculation of AC resistance of single‐layer coils using boundary‐element method

Correspondence Yao Luo, School of Electrical Engineering and Automation, Wuhan University, Wuhan, China. Email: sturmjungling@gmail.com Abstract Boundary‐element method (BEM) is applied to the analysis of AC resistance of single‐ layer coils. For simplification of the analysis, the coil windings are replaced by a system of parallel straight conductors. Applying the method of images and closed‐form solutions of off‐diagonal matrix elements, and exploiting the symmetry of submatrices, high computational efficiency can be achieved by the proposed approach. As a numerical method, it is suitable for solenoids or pancake coils evenly wound by round or rectangular wires, or other wire shapes with proper symmetry. The proposed method is verified by the numerical results compared with those of finite‐element method (FEM) and measurements, from which it is confirmed that the proposed approach has the same accuracy as the FEM, and the former is much faster and more memory efficient.


| INTRODUCTION
Single-layer coils are widely used in the applications such as the wireless power transmission (WPT) [1][2][3][4], in which the coils are normally operated at the frequencies higher than 10 kHz. At these frequencies, the AC resistance is considerably higher than the dc case due to the skin and proximity effects of the AC currents flowing in the windings. Consequently, the coil AC resistance accounting for the energy dissipation is crucial for the efficiency of the WPT. In the pioneering work of Sommerfeld [5], an analytical method has been given for the AC resistance of an infinite long solenoid. In the analysis, the curvature of the coil turns is neglected and the calculation is simplified to a twodimensional (2D) problem of infinite long parallel conductors. According to Ref. [5], this simplification is always recommended for the calculation of coil AC resistance. Many other efforts have also been made to predict the frequency behaviour of the AC resistance [6][7][8][9][10][11][12][13][14]. In general, the analytical methods are quite efficient and easy to use, but most of them are only suitable for certain frequency ranges, specific winding compositions, or limited shapes of the wire cross section [14].
On the other side, numerical methods are alternative choices allowing for a general analysis of the coil AC resistance without the limitations of the analytical counterparts. Among them is the finite-element method (FEM), which has been applied in the relevant problems [15][16][17][18], but difficulties arise for very small skin depth (special techniques are required for the volumetric meshing), and the background enclosure must be large enough (especially for air-core coils) to give reasonable results, this may demand huge amounts of computing resources. Another available approach has been proposed by Manneback and Noether [19][20][21][22][23][24][25], which is based on the integral equations and the unknowns are the current densities over each conductor, hence a 2D discretization is still needed in the solving process. In this regard, a more preferable approach is the boundary-element method (BEM) [26][27][28][29][30][31][32][33], by which only the discretization on the contour lines of the wire cross section is required and the computing resources can be saved greatly. Another remarkable advantage of BEM is that the AC resistance can be given directly from the solution of the assembly matrix without the post processing of Poynting vector. Moreover, for applications of WPT normally a short solenoid or pancake coil with several turns are utilized, so a method such as BEM considering the end effect of the coil turns is desirable for the accurate evaluations of the AC resistance.
However, in most of the relevant works with BEM only two or three conductors are investigated [30][31][32][33], so the formulations should be optimized to acquire high efficiency in the circumstances of more conductors. In this work, the analysis will be given on the basis of Ref. [33], and with the proposed optimization, the computation time can be reduced drastically. The proposed method will be tested at frequencies ranging from 1 kHz to 1 MHz, and the FEM and experimental data will also be provided for the comparisons.
It is also possible to make a three-dimensional (3D) analysis with the axisymmetric BEM for circular coils. Nevertheless, the accuracy improvement of the AC resistance will be marginal with considerable higher computational cost. Certainly, an axisymmetric BEM is necessary if the inductance is also required, but the frequency dependence of the inductance is relatively weak [6] and mostly the inductance obtained by the Neumann formula [34,35] is acceptable for the practical applications. Finally, this work is dedicated to the calculations for solid-wire coils. For Litz-wire coils, the readers are referred to Ref. [36] for more details.

| BOUNDARY INTEGRAL EQUATIONS FOR THE VECTOR POTENTIAL
Neglecting the curvature of the coil turns, a set of parallel conductors of infinite length is shown in Figure 1. A 2D analysis will be carried out in the Cartesian coordinate system (x, y). Then the vector potential A only has a z-component A(x, y) along the wires. This analysis is applicable to both solenoids and pancake coils due to the simplification. The coils are supposed to be made of the non-magnetic uniform material, such as copper or aluminium. Due to the source inside the conductive region, the vector potential A e will satisfy an inhomogeneous Helmholtz equation (using the time convention e iωt ) where, J s is the source current density, and where, ω, σ, and μ 0 are the angular frequency of the currents, conductivity of the conductor material and the permeability of the vacuum, respectively. Noting that for a straight conductor with voltage source we have where, U s is the per-unit-length (p.u.l.) source voltage drop (a constant), it can be deduced from Equation (1) where, A s is an undetermined constant containing the information of the total current over the conductor cross section. For BEM formulation a homogeneous equation is preferred, and an auxiliary vector potential A can be introduced for Equation (4): and Equation (1) can be converted into a homogeneous Helmholtz equation with the help of Equation (5): For the conductor k with the total current I k , the p.u.l. AC resistance is found from the comparison of Equations (1), (3), and (4): According to the BEM procedure we apply the Green's second theorem to Equation (6) [28] for the conductor k and obtain where, n is the outward normal unit vector, i is the index of the nodes, c i is a constant determined by the solid angle of the observation point P i , for constant element c i ¼ half; A k and Γ k are the auxiliary vector potential and contour of the conductor k, respectively, and is the fundamental solution of Equation (6), where r is the distance between the observation and the integration points, and K 0 (x) is the modified Bessel function of the second kind.
In the other side, the vector potential A e of the surrounding air satisfies the Laplace equation For a coil with N turns the boundary integral equation for A e can be given as where, the normal unit vector has the same direction as that in Equation (8), and Then, with the relation it can be deduced from Equation (11) By means of the relation [29].
> > > > > : if P i ; P j belong to the same conductor 0: if P i ; P j belong to different conductors where, P i and P j are the observation and integration points, respectively, it follows from Equation (14) For a complete definition of this problem, the relations for the total currents over the conductor cross sections are necessary [31,33]:

| SIMPLIFICATION OF BOUNDARY INTEGRAL EQUATIONS
It is well known that the BEM matrix is not sparse. To obtain an efficient algorithm, the symmetry of the conductor configuration should be exploited and the method of images is suitable for the reduction of the matrix dimension. Here, we use the single-layer 6-turn coils wound by rectangular wires (see Figures 3 or 4) to describe the reduction procedure. In principle the result can be generalized to coils with arbitrary number of turns. Firstly, the coil windings are replaced by a set of six parallel straight conductors of infinite length (see Figure 1). With the symmetry of y-axis, it is allowed to merely consider the conductors 1, 2 and 3, for which the boundary integral equations corresponding to Equation (8) can be given as F I G U R E 2 Geometric descriptions for variables used in the solutions of off-diagonal matrix elements where, is the fundamental solution in the conductor k, and r k ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi denoting the distance between the observation point (x k,i , y k,i ) and integration point (x k , y k ) on the same contour. By the method of images, the remaining equations for conductors 4, 5, and 6 are superfluous.
In addition, for the observation point on Γ 1 we have where, with r 1k ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi À representing the distance between the observation point (x 1,i , y 1,i ) and integration point (x k , y k ) on the contours Γ 1 and Γ k , respectively. Due to the Neumann condition (∂A/∂n ¼ 0) on the y-axis the integration quantities on the contours Γ 4 , Γ 5 and Γ 6 should be transformed with the mirror pairs: where, the distance r 1′3 in G 1′3 should be interpreted as ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi À In similar manner the integral on Γ 5 and Γ 6 can be replaced with the image integrals on Γ 2 and Γ 1 . Consequently, it follows from Equation (21) where, with the subscript "1" indicating the image coordinate (À x 1i , y 1i ). The remaining two integral equations can be obtained in the same manner where, A further simplification can be achieved by exploiting the symmetry of x-axis (the Neumann condition ∂A/∂y ¼ 0). It is allowed to only consider the lower half (the contours in the third quadrant, see Figure 1) and replace the fundamental solutions in Equations (18a, b, and c) with where, r s k ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Similarly, Equations (26a, b, and c) should be modified as where, Γ k * indicating the lower half of the contour of conductor k, and with r mk ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi À x m;i À x k � 2 þ � y m;i À y k � 2 r ð32aÞ r m ′ k ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi À x m;i þ x k � 2 þ � y m;i À y k � 2 r ð32bÞ r s mk ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi À x m;i À x k � 2 þ � y m;i þ y k � 2 r ð32cÞ r s m ′ k ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi À Additionally, it follows from Equation (17) that Accordingly, a linear system can be established by the discretization of Equations (18), (26), and (33) with the Green's functions Equations (28) and (31) where, the elements of the submatrices are Equations (37a)-(37j) are valid for even N. If N is odd, N/2 in Equations (37a)-(37j) should be replaced by (Nþ1)/2, and the conductor (Nþ1)/2 will be divided symmetrically by x ¼ 0 and y ¼ 0, then Γ (N þ 1)/2 * will become a quarter of the entire (N þ 1)/2-th contour, hence the final element in Equation (37j) should be replaced with I (Nþ1)/2 /2. Normally, all elements in Equation (37j) should be set to the same value. Equation (36) can be reduced to the one with lower dimension, as well, that is When m ¼ k and i ¼ j, analytical solutions are recommended to evaluate Equations (35a)-(35d) [37]: where, L n (x) is the modified Struve function of order n, and l j is the length of the jth element on the conductor k. Contrary to most literature about the calculation of constant elements, in fact, even for m≠k or i≠j the analytical solutions still exist for Equations (35c) and (35d). For example, the result relevant to the term r mk in Equation (35c) is where, ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi À x m;i À x k1 � 2 þ � y m;i À y k1 � 2 r ρ 2 ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi À x m;i À x k2 � 2 þ � y m;i À y k2 � 2 r and the meaning of the other variables in Equation (40) are shown in Figure 2. In addition, for the term relevant to r mk in Equation (35d) we have where, x f À x k;1 cos θ ; for x k;2 < x k;1 and θ ≠ π=2 y k;1 À y f ; for y k;2 > y k;1 and θ ¼ π=2 y f À y k;1 ; for y k;2 < y k;1 and θ ¼ π=2 and x k;2 À x f cos θ ; for x k;2 > x k;1 and θ ≠ π=2 x f À x k;2 cos θ ; for x k;2 < x k;1 and θ ≠ π=2 y k;2 À y f ; for y k;2 > y k;1 and θ ¼ π=2 y f À y k;2 ; for y k;2 < y k;1 and θ ¼ π=2 where, θ is the inclination angle of the line element to the x-axis. Considering the rows of H 0 and G 0 are densely occupied in the matrix Equation (36), the analytical matrix elements will be quite beneficial to the computational efficiency. The calculations can be accelerated further by using the symmetry of the submatrices in Equations (37a)-(37d). For a coil of uniform turn spacing, we have H 1 ¼ H 2 ¼⋯H N/2 , and G 1 ¼ G 2 ¼ ⋯G N/2 . More complicated relations also exist for the elements in Equations (37c) and (37d). For example, with N ¼ 6, denoting the components of H ij related to Equations (32a)-(32d) as H ij α , H ij β , H ij

| Numerical and experimental results for the solenoid and pancake coil wound by rectangular wires
The presented method has been verified by the experiment of impedance measurement, for which a solenoid and a pancake coil have been constructed. The inner radii of the solenoid and pancake coil are 122.5 and 80.25 mm, respectively. Both of them are six turns and made by the same aluminium alloy wires (σ ¼ 3.18 � 10 7 S/m) of rectangular cross section (a ¼ 4 mm, b ¼ 2 mm), and the turn spacing of the solenoid and pancake coil are c ¼ 2.24 mm and c ¼ 3.5 mm, respectively.An impedance analyser (GW Instek LCR-8220) has been used for the measurement of AC resistance in a frequency range spanning from 1 kHz to 1 MHz. The experimental sites are shown in Figures 3 and 4. The coil impedance is measured with the "L S -R S mode" using the frequency scanning function of the impedance analyser. The 4-wire probes (LCR-06B) used in the measurement are reliable at frequencies up to 1 MHz according to the information given by the manufacturer. It is notable that measurement errors will be included in the measured data due to the stray capacitance C P [38]: Hence, the obtained values from the measurements are R S , rather than R itself. In fact, parallel resonance will occur at the frequency ω r ¼ (LC P ) -1/2 . This will lead to significant overestimation of the AC resistance at high frequencies.
The measured data of the solenoid and pancake coil are illustrated in Figures 5 and 6 as red dots, and the BEM results are shown with black lines, accompanying with those of FEM denoted by the black circles obtained by the software package CST EM Studio. It is shown clearly that good agreement between the measurement and BEM results is obtained for both coils, and the BEM and FEM coincide very well over the full range of calculated frequency. The discrepancy is observed at the frequencies higher than 300 kHz. For pancake coil, the deviations occurring at lower frequencies may be due to the coil frame plate, which is made of phenolic resins with relative permittivity ε r ¼ 8 and comparatively higher stay capacitance could be anticipated.
In the FEM simulations, the coil windings are replaced with parallel straight conductors, as well. Thin slices of 0.1 mm thickness are used for the modelling of the cross sections of coil windings. About 8 � 10 5 4-node tetrahedral elements are used to obtain a result with the relative estimated error of about 5 � 10 À 4 , and the boundary condition E t ¼ 0 is imposed on the end faces of the slices to simulate the infinite long conductors. No remarkable improvement of accuracy is observed with a thicker model. On the other hand, the 3D fullscale simulation cannot be achieved due to the memory overflow and unrealistic simulation time.
The FEM simulations were implemented on a personal computer of a 3.5 GHz processor and 128 GB RAM memory. About 13 min are required to finish a simulation, making a strong contrast with the elapsed time of only about 3 s for an evaluation by the proposed method, with 30 constant elements for every winding. The memory consumption is about 40 GB for the FEM simulation and ignorable for our method.To demonstrate the skin and proximity effects more clearly, the AC resistance of every single wire is given for the solenoid and pancake coils. The results are shown in Tables 1 and 2, respectively, in which p.u.l. R AC1 , p.u.l. R AC2 , and p.u.l. R AC3 are AC resistance per unit length of the wires 1, 2, 3, and R Total is the AC resistance of the entire coil. Due to the symmetry only the values of half the coil wires are given. DC values of resistance are also provided here, which are obtained by the well-known equation R ¼ ρl/S. It should be noted that they are not deducible by the proposed method due to the kernel function K 0 (αr) included in the solving process (which becomes infinite at f ¼ 0). It is shown from Table 1 that the p.u.l. AC resistance of conductor 3 (the nearest one to the symmetric plane x ¼ 0) can be about two times higher than that of F I G U R E 5 AC resistance of the tested solenoid versus frequency F I G U R E 6 AC resistance of the tested pancake coil versus frequency LUO ET AL. conductor 1 (the outmost one) at high frequencies. This demonstrates a strong proximity effect between the turns of the experimental coils.

| Numerical results for the pancake coils wound by round wires
A further investigation of the proposed method can be given by calculations of three pancake coils of N ¼ 12, 14, and 16 turns, and of the inner radii 10, 20, and 30 mm, respectively. They are all wound by the same round copper wires (σ ¼ 5.6 � 10 7 S/m, wire radius ¼ 0.5 mm) with the same turn spacing of 0.5 mm (c ¼ 0.5 mm). In Figure 7, the AC resistance curves are plotted for these pancake coils obtained by our method compared with the FEM results represented by the dots. Very good agreement is achieved over the full range of the calculated frequencies. Additional comparisons are given for the AC resistance considering solely the skin effect (omitting the proximity losses), which are shown by the dashed  (44) and (45) 10lines. For the round wires considered here, the skin loss can be found by the analytical method. In fact, the AC current density in a round wire of radius a is where, I is the total current over the wire cross section, α is given by (2), and J n (x) is the Bessel function of the first kind and order n. Therefore, the p.u.l. AC resistance due to the skin effect in a single winding is where, the asterisk denotes the complex conjugate. By the comparisons with the dashed lines in Figure 7, it is shown clearly that the proximity loss becomes significant at the frequencies higher than 100 kHz, and its proportion increases gradually with frequency. The average elapsed time for an evaluation of the coils of N ¼ 12, 14, and 16 is about 5, 6, and 7 s, which is also minimal compared with that of FEM (ca., 3 min). The shorter elapsed time of simulations is ascribed to thinner wire considered here in comparison with the aforementioned rectangular wire. About 3 � 10 5 four-node tetrahedral elements are required to obtain a result with the relative estimated error of about 1 � 10 À 4 . The memory consumption is about 15 GB for the FEM and minimal for our method.

| CONCLUSION
The approximation of parallel straight conductors is suitable for the evaluations of the AC resistance of single-layer coils. Using the boundary element method, the AC resistance of the windings with round or rectangular cross sections can be easily obtained. High computational efficiency has been attained by virtue of a series of improvements such as the method of images, the closed-form solutions of the off-diagonal matrix elements, and the symmetry of the submatrices. This approach coincides well with the measured data below 300 kHz. At higher frequencies, discrepancies arise mainly due to the influence of stray capacitance in apparatus. Furthermore, the proposed approach proves to be highly consistent with the FEM in the entire examined frequency region, while the former is much faster than the latter. Analysis of coils wound by more realistic Litz wires is the subject of our future work.