Analytical-Numerical Approach to the Skin and Proximity E ﬀ ect in Lines with Round Parallel Wires

: Power and communication lines with round wires are often used in electrical engineering. The skin and proximity e ﬀ ects a ﬀ ect the current density distribution and increase resistances and energy losses. Many approaches were proposed to calculate the e ﬀ ects and related quantities. One of the simplest approximate closed solutions neglects the dimensions of neighboring wires. In this paper, a solution to this problem is proposed based on the method of successive reactions. In this context, the solution with substitutive ﬁlaments is considered as the ﬁrst approximation of the true solution. Several typical arrangements of wires in single-phase communication lines or three-phase bus ducts are considered to detect the limits of applicability of the ﬁrst approximation. The error of the ﬁrst approximation grows with wire radius to skin depth ratio and wire radius to wire spacing ratio. When the wire radius to skin depth ratio is up to 1, and the gap between the wires is above the wire radius, the error is at a level of 1%. However, lowering the distance and / or skin depth leads to a much larger error in the ﬁrst approximation.


Introduction
Power and communication lines with wires of circular cross-section (often called round wires) are very often used in electrical engineering. Two wires are used in single-phase cables or in communication technology. Three-phase cables with and without neutral wire are often arranged in flat or trefoil geometry. When an alternating current passes in a wire, it generates an alternating magnetic field, which induces eddy currents in the wire so that the current is pushed towards the conductor surface (skin effect) [1]. The current also generates eddy currents in neighboring conductors, so the current density in them is also disturbed (proximity effect). Therefore, the current density is considerably disturbed compared to the uniform density occurring during the direct current passage. As a result, the effective cross-section area drops and the resistances of the bus duct elements enlarge, which lowers bus duct ampacity and power transfer efficiency. Hence, the knowledge of current distribution in power and communication lines is important to properly assess or optimize their properties.
Despite the fact that the phenomenon is well known, and much research has been done in this field, the topic is still vivid, and new methods are still being developed, e.g., [2]. When analyzing known approaches, it seems they can be categorized into differential, integral and circuit-based. In the first approach, the subject is reduced to solving the Helmholtz equation in the conductors and the Laplace equation in the non-conducting regions. However, except for coaxial conductors, the equations are hard to solve. For example, when considering a twin line with round wires, one would like to use the method of separation of variables in bipolar coordinates. Unfortunately, the variables in the Helmholtz equation do not separate in that coordinate system. The situation is even worse for systems with multiple wires. Therefore, the differential approach requires either some simplification or numerical support. One of the first analytical solutions was found by Carson [3], who considered wave propagation along two parallel round wires. The solution was expressed as infinite series with coefficients, which have to be found via successive approximations. Later Arnold [4] used this result as well as some others to obtain formulae for the proximity factor defined as AC resistance to DC resistance for round parallel conductors in single and three-phase cases. There were many approaches with semi-analytical solutions, where fields are expressed as infinite series, yet the coefficients in the series have to be determined numerically by solving a system of linear algebraic equations. For example, Hussain and Bringer [5] obtained E and H fields for a system of parallel conductors, Ahmed et al. [6] used the method of fundamental solutions to express E and H fields as a linear combination of such solutions, and Machado et al. [7] used infinite series together with multipole expansion in the dielectric. Nevertheless, the final calculations always require the determination of coefficients in the series, which is mostly done via series truncation and solving a linear system of equations (but successive approximations are possible, too). Therefore, the most general approach is just using a purely numerical approach, like the finite element method (FEM) (e.g., [8][9][10][11]) or boundary element method (BEM) (e.g., [12]). Furthermore, combined models were used, like FEM + BEM [13] or Bubnov-Galerkin + variable separation method [14,15]. A pure numerical approach is quite general but does not allow us to obtain closed formulae, which makes it difficult to formulate general conclusions.
The second approach uses an integral formulation. One of the earliest solutions was found by Manneback [16], who solved an integral equation for eddy currents in a round conductor due to nearby current filament. An interesting approach useful for high frequencies, when currents are constrained to a thin surface layer, was considered by Smith [17]. Dlabač and Filipović [18] solved the integral equation for two round parallel wires by expressing the solution in the form of power series; of course, this required series truncation and solving a linear system of algebraic equations. The most general approach leads to the Fredholm integral equation of the second kind. Again, it can be solved numerically, e.g., via the collocation method, like in [19].
Finally, the circuit-based approach uses the idea of splitting the conductors into a certain set of parallel subconductors. Then each subconductor is assigned DC resistance and inductance based on assumed current density, which is usually considered constant (although unknown) if the subconductor's cross-section area is small enough. Then equations known from circuit theory are used to formulate a system of algebraic equations [20]. In a way, a similar approach was proposed by Coufal [2,21], but with voltage excitation, which is more realistic. Nevertheless, all the approaches require some kind of numerical support. The Manneback solution was later used to model a set of round wires by treating the neighboring wires as filaments, e.g., [22,23]. This is one of the simplest approximate approaches. However, the main drawback is not taking into account the finite dimensions of neighboring wires. The method of including the dimensions into the solution was presented independently in [24,25]. In this work, the method is generalized for lines with multiple wires. In this approach, the solution with substitutive current filaments is regarded as the first approximation of the true solution. Several typical lines with round wires are considered as examples to check the performance of the method as well as the limits of applicability of the first approximation.

Problem Description
Let us consider a multi-wire bus duct with round conductors (Figure 1). The bus duct consists of K long wires placed parallel to each other in a nonconductive and non-magnetic environment. The length is assumed much larger than the distances between the wires. All the wires are non-magnetic, have constant conductivity σ k and invariable circular cross-section of radius R k , where k = 1, 2, . . . , K. It is assumed that the wires carry sinusoidal currents with angular frequency ω and complex root-mean-square (RMS) value I k . The frequency is assumed low enough to neglect the displacement currents (electromagnetic wave length is much larger than circuit size).
Energies 2020, 13, x FOR PEER REVIEW 3 of 21 where = 1,2, … , . It is assumed that the wires carry sinusoidal currents with angular frequency and complex root-mean-square (RMS) value ℐ . The frequency is assumed low enough to neglect the displacement currents (electromagnetic wave length is much larger than circuit size).

Figure 1.
Multi-wire bus duct with round parallel wires.
The above assumptions allow us to use 2D calculations in the cross-section of the configuration. Due to the default assumption of linearity, the current density in wire can be expressed as follows: where are functions of size, wire position, frequency, and material properties, but independent of currents in the wires. It is convenient to separate the part dependent on current in wire . This part is related to the so-called skin effect (superscript "s"). The remaining part is attributed to eddy currents induced in wire due to current passage in the neighboring wires, which is called the proximity effect (superscript "p"). Finding closed exact analytical forms for is possible only in special cases with coaxial symmetry. Another geometry requires either a kind of numerical approach or introducing additional simplifications.

Approximate Analytical Solution
In the case of a system with multi-wire round conductors, the approximate analytical solution can be found, among others, in [22]. The approach is based on introducing simplifications and treating the neighboring wires as infinitely thin filaments. This leads to a fully analytical solution to a substitutive problem and is often sufficient in practice. This approximate solution will be temporarily marked with a tilde. It is convenient to use polar coordinates ( , ) associated with wire -see Figure 2. If ( , ) are coordinates of the axis of wire then it follows that: where is the modified Bessel function of the first kind of order , and: is the propagation constant, in which = √2 ( 0 ) ⁄ is the skin depth for wire and j is the imaginary unit. Function Λ equals: and is interpreted as eddy current density induced at point of polar coordinates ( , ) in a unit radius round wire aligned with axis and having propagation constant due to a parallel current filament with the current π crossing point ( , 0)-see Figure 3. The above assumptions allow us to use 2D calculations in the cross-section of the configuration. Due to the default assumption of linearity, the current density in wire i can be expressed as follows: where W ik are functions of size, wire position, frequency, and material properties, but independent of currents in the wires. It is convenient to separate the part dependent on current in wire i. This part is related to the so-called skin effect (superscript "s"). The remaining part is attributed to eddy currents induced in wire i due to current passage in the neighboring wires, which is called the proximity effect (superscript "p"). Finding closed exact analytical forms for W ik is possible only in special cases with coaxial symmetry. Another geometry requires either a kind of numerical approach or introducing additional simplifications.

Approximate Analytical Solution
In the case of a system with multi-wire round conductors, the approximate analytical solution can be found, among others, in [22]. The approach is based on introducing simplifications and treating the neighboring wires as infinitely thin filaments. This leads to a fully analytical solution to a substitutive problem and is often sufficient in practice. This approximate solution will be temporarily marked with a tilde. It is convenient to use polar coordinates (r, θ) associated with wire i-see Figure 2. If (r k , θ k ) are coordinates of the axis of wire k then it follows that: where I n is the modified Bessel function of the first kind of order n, and: is the propagation constant, in which δ i = 2/(ωµ 0 σ i ) is the skin depth for wire i and j is the imaginary unit. Function Λ equals: and is interpreted as eddy current density induced at point of polar coordinates (ρ, ϕ) in a unit radius round wire aligned with z axis and having propagation constant γ due to a parallel current filament with the current π crossing point (ξ, 0)-see Figure 3.  Such an approach neglects the finite dimensions of the cross-section of the adjacent conductors. The question is, what error causes such neglect. To estimate this, a comparison with the exact solution is needed. However, there is no such method. One way is to use a solution obtained via other methods (e.g., finite elements) assumed as more accurate. Another approach is to construct a more precise approximation of the original problem solution. Both methods were used in the research.

Improving the Solution
The true solution differs from that given by Equations (1)-(5) because the parts of current in a neighboring wire which are closer to the considered wire interacts stronger than it is expressed by Equation (3). In these terms, the above solution can be regarded as a first approximation of the true solution. A more accurate approximation can be obtained via the method of successive reactions. The eddy currents induced in a wire due to the neighboring substitutive filaments are now a new source of the magnetic field, which affects the neighboring wires themselves and induces additional eddy currents. The additional currents are another source of a magnetic field and induce new eddy currents and so forth. Of course, this image is just an auxiliary construction of our minds. In fact, only the total of the reactions can be observed. However, by calculating the successive reactions, we approach the true solution.
To formalize this idea, let us introduce superscripts ( ) and [ ] to denote -th reaction and  Such an approach neglects the finite dimensions of the cross-section of the adjacent conductors. The question is, what error causes such neglect. To estimate this, a comparison with the exact solution is needed. However, there is no such method. One way is to use a solution obtained via other methods (e.g., finite elements) assumed as more accurate. Another approach is to construct a more precise approximation of the original problem solution. Both methods were used in the research.

Improving the Solution
The true solution differs from that given by Equations (1)-(5) because the parts of current in a neighboring wire which are closer to the considered wire interacts stronger than it is expressed by Equation (3). In these terms, the above solution can be regarded as a first approximation of the true solution. A more accurate approximation can be obtained via the method of successive reactions. The eddy currents induced in a wire due to the neighboring substitutive filaments are now a new source of the magnetic field, which affects the neighboring wires themselves and induces additional eddy currents. The additional currents are another source of a magnetic field and induce new eddy currents and so forth. Of course, this image is just an auxiliary construction of our minds. In fact, only the total of the reactions can be observed. However, by calculating the successive reactions, we approach the true solution.
To formalize this idea, let us introduce superscripts ( ) and [ ] to denote -th reaction and Such an approach neglects the finite dimensions of the cross-section of the adjacent conductors. The question is, what error causes such neglect. To estimate this, a comparison with the exact solution is needed. However, there is no such method. One way is to use a solution obtained via other methods (e.g., finite elements) assumed as more accurate. Another approach is to construct a more precise approximation of the original problem solution. Both methods were used in the research.

Improving the Solution
The true solution differs from that given by Equations (1)-(5) because the parts of current in a neighboring wire which are closer to the considered wire interacts stronger than it is expressed by Equation (3). In these terms, the above solution can be regarded as a first approximation of the true solution. A more accurate approximation can be obtained via the method of successive reactions. The eddy currents induced in a wire due to the neighboring substitutive filaments are now a new source of the magnetic field, which affects the neighboring wires themselves and induces additional eddy currents. The additional currents are another source of a magnetic field and induce new eddy currents and so forth. Of course, this image is just an auxiliary construction of our minds. In fact, only the total of the reactions can be observed. However, by calculating the successive reactions, we approach the true solution.
Function J i represents the approximate solution presented in Section 2.2, whereas ∆J [M] i is the correction due to the finite size of the wires after considering M higher-order reactions.
Reaction m ≥ 2 in wire i can be calculated by solving the equation: and then calculating: The point is that Equation (9) is not easy to solve. If it were, the whole procedure would not be necessary because the true solution could be obtained in one-step. The proposition is to represent J (m−1) j at point X as follows: where δ(X − Y) is a 2D Dirac delta function placed at point Y. Physically, this is equivalent to think of wire j as an infinite set of infinitesimal filaments with currents J (m−1) j (Y)dY. Solution for a round wire and one filament is, however, analogical to that which is given by Equation (3). Hence, by the method of superposition reaction m in current density can be expressed as follows ( Figure 4):

Numerical Implementation
Integral (12) is hard to evaluate analytically; therefore, some numerical support will be required. In its simplest form, wire, j can be divided into Q "small" fragments of areas S j,q and then: where Y j,q is the formal center of the fragment S j,q . Let us denote J Y j,q and introduce designation: Energies 2020, 13, 6716 Then a numerical version of Equation (12) can be rewritten as follows: It is convenient to arrange the discrete current densities J as well as introduce matrices Λ ij with elements Λ i,j;p,q . Then it follows that: with J (0) i and J (1) i evaluated from Equations (7) and (8). Such an approach is equivalent to representing each wire as a set of the finite number of current filaments distributed throughout the wire cross-sections. Hence, wire j, when acting on the neighboring wires, can be modeled as a set of Q current filaments placed at points Y j,q and carrying currents J and then calculating: The point is that Equation (9) is not easy to solve. If it were, the whole procedure would not be necessary because the true solution could be obtained in one-step. The proposition is to represent ( −1) at point as follows: where ( − ) is a 2D Dirac delta function placed at point . Physically, this is equivalent to think of wire as an infinite set of infinitesimal filaments with currents ( −1) ( ) . Solution for a round wire and one filament is, however, analogical to that which is given by Equation (3). Hence, by the method of superposition reaction in current density can be expressed as follows ( Figure 4):

Numerical Implementation
Integral (12) is hard to evaluate analytically; therefore, some numerical support will be required. In its simplest form, wire, can be divided into "small" fragments of areas , and then: where , is the formal center of the fragment , . Let us denote , ( −1) = ( −1) ( , ) and introduce designation: Energies 2020, 13, x FOR PEER REVIEW 6 of 21 Then a numerical version of Equation (12) can be rewritten as follows: It is convenient to arrange the discrete current densities , ( ) in vectors ( ) as well as introduce matrices with elements Λ , ; , . Then it follows that: with (0) and (1) evaluated from Equations (7) and (8). Such an approach is equivalent to representing each wire as a set of the finite number of current filaments distributed throughout the wire cross-sections. Hence, wire , when acting on the neighboring wires, can be modeled as a set of current filaments placed at points , and carrying currents , , ( Figure 5). Having found reaction , we can use it to construct reaction + 1 and so forth. Let us observe that this procedure is semi-analytical-it gives an analytical solution but requires a bit of a numerical approach. However, in contrast to many numerical methods, it does not require solving a system of algebraic linear equations but only performing a summation. Moreover, matrices , are calculated only once, and then they can be used to find corrections of any order. The successive corrections are weaker and weaker so that the calculations can be finished when correction is weak enough. The algorithm can be summarized as follows: 1. Select a set of filaments representing the wires and calculate matrices , ; (1) Having found reaction m, we can use it to construct reaction m + 1 and so forth. Let us observe that this procedure is semi-analytical-it gives an analytical solution but requires a bit of a numerical approach. However, in contrast to many numerical methods, it does not require solving a system of algebraic linear equations but only performing a summation. Moreover, matrices Λ i,j are calculated only once, and then they can be used to find corrections of any order. The successive corrections are weaker and weaker so that the calculations can be finished when correction m is weak enough. The algorithm can be summarized as follows: For m = 2, 3, . . . to obtain presumed accuracy or to consider the presumed number of corrections M;

4.
Use vectors ∆J to represent the correction of current density in wire i as follows:

5.
Repeat steps 2-4 for other currents in wires if necessary.

Power Losses and Resistance
The power losses in wire i per unit of length can be calculated using the standard integral over wire cross-sections as follows: In the simplest approach, this integral can be calculated using the decomposition of wire cross-section into sectors related to the filaments. A more sophisticated but also more accurate is using Equation (6) in the above formula and perform analytical integration. Brief information on this is presented in Appendix A. Then power can be used to find resistance, e.g., in a twin line with opposing currents (P 1 + P 2 )/I 2 .

Analyzed Configurations
Several typical arrangements of current bus ducts with round wires are considered below. They are as follows: ) and zero sequences (I L3 = I L2 = I L1 ) are analyzed. In the case when the total of currents is not zero, it is assumed the current returns via additional wire placed far from the considered ones.
In each case, it is assumed that the wires are identical, and their radius equals R. The value of Γ is established as (1 + j)/δ, where δ is the skin depth. The spacing between the closest wires in each arrangement is equal to s (the distance between wire axes is d = s + 2R). Analysis of Equation (12) indicates that current density depends on R/δ and d/R ratios.  In the case of the single-phase line, two versions are analyzed: opposing and same currents of equal RMS value. In the case of three-phase lines, versions of equal RMS currents and positive (ℐ 2 = ℐ 1 exp(−j120°), ℐ 3 = ℐ 2 exp(−j120°)), negative (ℐ 2 = ℐ 1 exp(+j120°), ℐ 3 = ℐ 2 exp(+j120°)) and zero sequences (ℐ 3 = ℐ 2 = ℐ 1 ) are analyzed. In the case when the total of currents is not zero, it is assumed the current returns via additional wire placed far from the considered ones.
In each case, it is assumed that the wires are identical, and their radius equals . The value of Γ is established as (1 + j)⁄ , where is the skin depth. The spacing between the closest wires in each arrangement is equal to (the distance between wire axes is = + 2 ). Analysis of Equation (12) indicates that current density depends on ⁄ and ⁄ ratios.

Comparison with Literature
As a test example, a twin line with opposing currents was selected, and the results were compared with those obtained in [2], Problem 2. The parameters of the line were as follows: wire radius = 10 mm, distance between wire axes = 40 mm, wires of hypothetical conductivity ( Al = 36.9/c MS/m), where is a constant, frequency = 60 Hz. Figure 7a shows the current density RMS value on the left wire surface for = 1, which corresponds to / = 0.935. The values are expressed in relation to DC current density when the total voltage drop on the line equals 1 V/m. The reference values, which were obtained by analyzing [2], were marked with dots. Figure 7b shows the current density RMS value on the diameter of the

Comparison with Literature
As a test example, a twin line with opposing currents was selected, and the results were compared with those obtained in [2], Problem 2. The parameters of the line were as follows: wire radius R = 10 mm, distance between wire axes d = 40 mm, wires of hypothetical conductivity (σ Al = 36.9/c MS/m), where c is a constant, frequency f = 60 Hz. Figure 7a shows the current density RMS value on the left wire surface for c = 1, which corresponds to R/δ = 0.935. The values are expressed in relation to DC current density when the total voltage drop on the line equals 1 V/m. The reference values, which were obtained by analyzing [2], were marked with dots. Figure 7b shows the current density RMS value on the diameter of the left wire for c = 1, 0.05, 0.02 and 0.01 as in work [2]. The values correspond to R/δ ratios equal to 0.935, 4.18, 6.61 and 9.35, respectively. The current density was expressed in absolute units of A/mm 2 . In both cases, the full agreement was obtained. left wire for = 1, 0.05, 0.02 and 0.01 as in work [2]. The values correspond to / ratios equal to 0.935, 4.18, 6.61 and 9.35, respectively. The current density was expressed in absolute units of A/mm 2 . In both cases, the full agreement was obtained.    Table 1 shows the resistance per length of the unit for selected values of c. It was calculated via Equation (18) and compared with that given in work [2], Table 1, as well as with that obtained via finite element method (FEMM software, version 4.2, 64-bit 21 April 2019 by David Meeker, MA, USA, http://www.femm.info/wiki/HomePage) with very fine mesh. When the skin depth is comparable with wire radius, the results agree very well. For smaller skin depths, some discrepancies arise. This is probably the result of too rough discretization (each wire was divided into 100 filaments).  Figures 8 and 9 show current density distribution on wire surfaces for configurations 6a and 6c, respectively, for R/δ = 3 and d/R = 2.5 with selected variants of currents. The solid lines represent the improved solution with higher-order reactions up to M = 6. Each wire was divided into 100 filaments to ensure high accuracy. The solutions were compared with the results obtained by means of FEM. The FEM solution is represented by dots in Figures 7 and 8. It follows that the improved semi-analytical solution agrees very well with FEM results. The dashed lines represent the first approximation. The differences between the first approximation and the improved solution are noticeable. Further analysis reveals that the differences grow with R/δ ratio and drops with d/R ratio. Figure 10 shows the plots of current density correction given by Equation (17) in wire cross-sections for selected configurations and supply variants. The values are expressed in percentage of uniform current density. It follows that the largest differences appear on wire surfaces. Analysis of analogous plots for the remaining configurations (not attached) leads to similar conclusions.

The Effect of Wire Discretization
The improved method requires using filaments to represent the neighboring wires. The more filaments used, the higher accuracy obtained, but also the more time required in computations. In this section, the effect of filament number is tested. The wires are divided into n sectors of similar areas. Figure 11 shows the tested discretization variants.  Figure 10 shows the plots of current density correction given by Equation (17) Figure 10 shows the plots of current density correction given by Equation (17) Figure 10 shows the plots of current density correction given by Equation (17)

The Effect of Wire Discretization
The improved method requires using filaments to represent the neighboring wires. The more filaments used, the higher accuracy obtained, but also the more time required in computations. In this section, the effect of filament number is tested. The wires are divided into sectors of similar areas. Figure 11 shows the tested discretization variants.  Figure 12 shows the relative error in current density modulus on the wire surface in a twin line with opposing and similar currents. The left part of each plot corresponds to opposing currents, whereas the right part is for similar currents. The four plots correspond to weak and strong skin effects as well as closely or loosely placed wires. It follows that when the skin effect is weak ( ≥ ) and the wires are not too close ( ≥ 3 ), the error is below 1% even for the first approximation ( Figure 12a). However, for closely placed wires ( = 2.2 ) to keep the error at a level of 1%, it is necessary to take at least four filaments. In case of stronger skin effect ( ≤ 3 ⁄ ) and not too closely placed wires, the number of filaments should be around 16. When the wires are placed closely, even more filaments are required. It is probably possible to propose a more efficient discretization. However, the results indicate that in low-frequency cases with not too close wires, it is usually enough to use the first approximation, whereas this approximation is not enough for higher frequency cases. As the frequency rises, the skin depth decreases and the eddy currents due to neighboring wires concentrate mainly near the wire surface. Therefore, eddy current density is expected larger near the wire surface, whereas very small in deeper layers. In such cases, it would be more efficient to use a denser mesh of filaments near wire surfaces with fewer filaments near the wire axis. Of course, it should be kept in mind that this model is valid for frequencies at which displacement currents can be neglected, as stated in Section 2.1. (d) Figure 10. The values of correction given by Equation (17)

The Effect of Wire Discretization
The improved method requires using filaments to represent the neighboring wires. The more filaments used, the higher accuracy obtained, but also the more time required in computations. In this section, the effect of filament number is tested. The wires are divided into sectors of similar areas. Figure 11 shows the tested discretization variants.  Figure 12 shows the relative error in current density modulus on the wire surface in a twin line with opposing and similar currents. The left part of each plot corresponds to opposing currents, whereas the right part is for similar currents. The four plots correspond to weak and strong skin effects as well as closely or loosely placed wires. It follows that when the skin effect is weak ( ≥ ) and the wires are not too close ( ≥ 3 ), the error is below 1% even for the first approximation ( Figure 12a). However, for closely placed wires ( = 2.2 ) to keep the error at a level of 1%, it is necessary to take at least four filaments. In case of stronger skin effect ( ≤ 3 ⁄ ) and not too closely placed wires, the number of filaments should be around 16. When the wires are placed closely, even more filaments are required. It is probably possible to propose a more efficient discretization. However, the results indicate that in low-frequency cases with not too close wires, it is usually enough to use the first approximation, whereas this approximation is not enough for higher frequency cases. As the frequency rises, the skin depth decreases and the eddy currents due to neighboring wires concentrate mainly near the wire surface. Therefore, eddy current density is expected larger near the wire surface, whereas very small in deeper layers. In such cases, it would be more efficient to use a denser mesh of filaments near wire surfaces with fewer filaments near the wire axis. Of course, it should be kept in mind that this model is valid for frequencies at which displacement currents can be neglected, as stated in Section 2.1.  Figure 12 shows the relative error in current density modulus on the wire surface in a twin line with opposing and similar currents. The left part of each plot corresponds to opposing currents, whereas the right part is for similar currents. The four plots correspond to weak and strong skin effects as well as closely or loosely placed wires. It follows that when the skin effect is weak (δ ≥ R) and the wires are not too close (d ≥ 3R), the error is below 1% even for the first approximation (Figure 12a). However, for closely placed wires (d = 2.2R) to keep the error at a level of 1%, it is necessary to take at least four filaments. In case of stronger skin effect (δ ≤ R/3) and not too closely placed wires, the number of filaments should be around 16. When the wires are placed closely, even more filaments are required. It is probably possible to propose a more efficient discretization. However, the results indicate that in low-frequency cases with not too close wires, it is usually enough to use the first approximation, whereas this approximation is not enough for higher frequency cases. As the frequency rises, the skin depth decreases and the eddy currents due to neighboring wires concentrate mainly near the wire surface. Therefore, eddy current density is expected larger near the wire surface, whereas very small in deeper layers. In such cases, it would be more efficient to use a denser mesh of filaments near wire surfaces with fewer filaments near the wire axis. Of course, it should be kept in mind that this model is valid for frequencies at which displacement currents can be neglected, as stated in Section 2.1.

The Effect of Skin Depth
To assess the effect of skin depth, all the considered cases were analyzed for constant distance d = 3R and R/δ ratio varying in range 0.1 to 4. The correction given by Equation (17) related to uniform current density was calculated at some characteristic points of all the configurations at selected variants of supply currents. The results are plotted in Figures 13 and 14. It follows that for δ > R the correction modulus rises approximately proportionally to R/δ ratio. The differences depend not only on skin depth to radius ratio but also on phases of currents in the wires. In the case of three-phase lines, the phase sequence is important. When the current distortion is high, there are significant harmonics with high R/δ ratio and various phase sequence. In such a case, the first approximation may introduce significant errors. selected variants of supply currents. The results are plotted in Figures 13 and 14. It follows that for > the correction modulus rises approximately proportionally to ⁄ ratio. The differences depend not only on skin depth to radius ratio but also on phases of currents in the wires. In the case of three-phase lines, the phase sequence is important. When the current distortion is high, there are significant harmonics with high ⁄ ratio and various phase sequence. In such a case, the first approximation may introduce significant errors.

The Effect of Spacing between Wires
The spacing between the wires also greatly affects the current distribution. Figures 15 and 16 show the modulus of current density correction given by Equation (17) versus distance between the wires at weak skin effect ( = ). It follows that the correction modulus is approximately inversely proportional to the distance between axes of the wires. For weak skin effect and distance between the wires not smaller than wire radius, the correction seems negligible. However, as shown in Figures 13 and 14, the correction grows with the skin depth, so for higher harmonics, the correction may be significant.

The Effect of Spacing between Wires
The spacing between the wires also greatly affects the current distribution. Figures 15 and 16 show the modulus of current density correction given by Equation (17) versus distance between the wires at weak skin effect (δ = R). It follows that the correction modulus is approximately inversely proportional to the distance between axes of the wires. For weak skin effect and distance between the wires not smaller than wire radius, the correction seems negligible. However, as shown in Figures 13 and 14, the correction grows with the skin depth, so for higher harmonics, the correction may be significant.

Power Losses
The results in Sections 3.5 and 3.6 indicate that using the approximate analytical solution presented in Section 2.2 leads to errors in current density, which are the bigger, the closer the wires and the smaller the skin depth. Current density is directly related to power losses via Equation (18). Therefore, power losses were calculated both with higher-order reactions taken into account ( ) and without them as resulting from the first reaction ( [1] ), and the difference was then calculated. Figure  17 shows the difference presented as % = ( − [1] ) [1] ⁄ × 100%. Power was used instead of resistance because the latter leads to ambiguity for three-phase bus ducts, which would require additional explanations. In cases such as twin lines with opposing currents, the error in resistance

Power Losses
The results in Sections 3.5 and 3.6 indicate that using the approximate analytical solution presented in Section 2.2 leads to errors in current density, which are the bigger, the closer the wires and the smaller the skin depth. Current density is directly related to power losses via Equation (18). Therefore, power losses were calculated both with higher-order reactions taken into account (P) and without them as resulting from the first reaction (P [1] ), and the difference was then calculated. Figure 17 shows the difference presented as δP% = P − P [1] /P [1] × 100%. Power was used instead of resistance because the latter leads to ambiguity for three-phase bus ducts, which would require additional explanations. In cases such as twin lines with opposing currents, the error in resistance will be the same as δP%, because power losses are proportional to resistance. The calculations were performed for 4 sets of (R/δ, d/R) ratios: (1,2), (1,4), (4,2), (4,4). In all the cases, typical variants of currents in the bus duct were taken into account. Figure 17 shows cases for which δP% values were at least 1%. As expected, the most significant errors occur for the highest tested value of R/δ and closely placed wires, i.e., for parameters (4,2). In such cases, the errors are nearly 30% in some variants, but typically up to around 15%. The errors for a twin line are relatively small compared to those for three-phase lines. Of course, in practice, the wires are not placed so tightly. Therefore errors are then smaller. For d/R = 4 the errors are usually below 1%. (a-f as in Figure 6); numbers in parentheses show and ratios and the preceding number is the number of the wire, where 4 is attributed to N wire. Figure 17. Percentage error δP% for selected variants of supply in the considered configurations (a-f as in Figure 6); numbers in parentheses show R δ and d R ratios and the preceding number is the number of the wire, where 4 is attributed to N wire.

Conclusions
A semi-analytical method for a set of parallel round wires was proposed in the paper. It can be regarded as an extension of the analytical approach in which the neighboring wires are treated as current filaments, which is called here the first approximation. The first approximation leads to the fully analytical formula with neglecting the dimensions of the neighboring wires. The proposed method also leads to analytical formula, but the dimensions of the neighboring wires are taken into account. This, however, usually requires some numerical computations. In contrast to many numerical methods, it does not require solving an algebraic system of equations.
In the further stage of research, the method was used for several typical configurations of singleand three-phase lines. The results were compared with those obtained by means of finite elements, and a full agreement was obtained. The results were also compared with those from the first approximation as well as to results found in the literature. The first approximation offers accuracy in current density about 1% when the skin depth is not less than the wire radius and when the gap between the wires is not less than the wire radius. The error is the bigger, the closer the wires are and the smaller the skin depth. The analysis revealed that if the skin depth is smaller than the wire radius, the error increases almost proportionally to the skin depth. This may introduce errors in power losses and resistances of bus duct elements, especially for a higher frequency, and in the case of three-phase systems-for currents with considerable levels of higher harmonics.
Like all methods, also the considered one has its pros and cons. It offers a semi-analytical solution, which sometimes can be useful. There is no necessity to solve an algebraic system of equations because the successive reactions quickly tend to zero. However, the numerical effort in this method is quite noticeable because function Λ must be evaluated many times. When the skin depth is very small compared to the wire radius, the required discretization may be too large to effectively perform calculations with reasonable accuracy and in an acceptable time. Therefore, it seems that the most suitable application is for skin depths above R/5 and closely placed wires (as for distant wires, the first approximation offers quite an acceptable accuracy).

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Using Equation (6) in Equation (18) The first term corresponds to the first approximation. The two other are corrections due to higher-order reactions. After putting Equations (7), (8) and (17) into the above formula, the integrals can be evaluated analytically because cosine in Equation (5) makes the mixed terms in the series products integrate to zero. The final formula for the power can be expressed as follows: I k I * l C i,k,l + 2 K k=1,k i w i 0 (Γ i R i ) for k = l = i, Θ(ξ ki ξ li , θ ki − θ li , Γ i R i ) for k i and l i, 0 otherwise, (A3) In the above formulas, w and v stand for any of additional filaments representing the wires except wire i. Parameter ξ ki = |X k − X i |/R i is the relative distance of filament k from the axis of wire i, and θ ki is the angular coordinate of filament k in the local coordinate system connected with wire i. Functions i n and Θ are defined as follows: i n (γ) = γI n (γ) Θ(ζ, θ, γ) = ∞ n=1 ζ −n i n (γ) cos nθ. (A6)