2D electrostatic modeling of twisted pairs

Abstract

Purpose

This paper aims to model a three-dimensional twisted geometry of a twisted pair studied in an electrostatic approximation using only two-dimensional (2D) finite elements.

Design/methodology/approach

The proposed method is based on the reformulation of the weak formulation of the electrostatics problem to deal with twisted geometries only in 2D.

Findings

The method is based on a change of coordinates and enables a faster computational time as well as a high accuracy.

Originality/value

The effectiveness of the adopted approach is demonstrated by studying different configurations related to the IEC 60851-5 standard defined for the measurement of the electrical properties of the insulation of the winding wires used in electrical machines.

Citation

Hazim, K., Parent, G., Duchesne, S., Nicolet, A. and Geuzaine, C. (2022), "2D electrostatic modeling of twisted pairs", COMPEL - The international journal for computation and mathematics in electrical and electronic engineering, Vol. 41 No. 1, pp. 48-63. https://doi.org/10.1108/COMPEL-04-2021-0139

1. Introduction

Using twisted wires in electrical engineering goes back to the late 19th century, when Bell introduced them to mitigate the crosstalk between the first telephone and a telegraph (Bell, 1876). They soon became commonplace in many electric and electronic applications and are in particular used as a means to test the aging of the insulation of electrical machines (Guastavino and Dardano, 2012) and to extract the electrical properties of their insulation (IEC, 2019).

Several studies have investigated the effect of twisted wires on the surrounding environment (Moser and Spencer, 1968), a nearby single wire (Paul and Jolly, 1982) or the ground (Pignari and Spadacini, 2011), by using predictive models. However, the accurate calculation of the distributed capacitances or the local electromagnetic fields, e.g. the local electric field strength for realistic configurations, requires costly three-dimensional (3D) simulations (Acero et al., 2014; Lyly et al., 2012). While a method based on several two-dimensional (2D) simulations (on slices of the 3D geometry) was proposed in Gustavsen et al. (2009), it introduces an intrinsic modeling error because of the performed averaging and cannot recover the true 3D local fields.

In this paper, following the approach originally derived in Nicolet et al. (2006, 2007a), we propose to reformulate the 3D problem in helicoidal instead of Cartesian coordinates (Waldron, 1958), leading to a 2D formulation with an adapted metric, easy to implement in existing 2D finite element codes. For purely helicoidal and infinitely long 3D geometries, the proposed 2D formulation is exact. The gains are twofold: the geometrical modeling and meshing of the potentially complex multi-layer conductors usually with large aspect ratios can be performed on a 2D slice instead of in 3D; and the finite element solution is greatly reduced.

The paper is organized as follows. After describing the geometrical configuration of twisted pairs in Section 2, the 3D electrostatic problem is reformulated in helicoidal coordinates in Section 3, leading to a 2D finite element formulation with an adapted metric – the “2D twisted” model. This 2D twisted model is verified against the full 3D model in Section 4 by computing the conductor charges in function of the number of turns in a twisted pair of enameled wires. Next, the improvement in accuracy of the 2D twisted model compared to the usual 2D (straight) model is analyzed in Section 5 in function of the number of turns, the conductor diameter and the insulation thickness. The effectiveness of the approach is finally demonstrated by studying different configurations related to the IEC 60851-5 standard for the measurement of the electrical properties of winding wires.

2. Geometry of twisted wires

The twisted geometry is defined by a number of turns N_turns over a swirl length ℓ_s (Figure 1). One spatial period of the twisted pair is depicted in Figure 2, where the extrusion length ℓ_e (i.e. the linear length of the spatial period) and the length of the helix ℓ_h (i.e. the length of the mean fiber of each twisted conductor) are defined as follows:

(1) {ℓe=ℓsNturnsℓh=(π(d+2ℓi))2+ℓe2

with d and ℓ_i being the diameter of the conductors and the insulation thickness, respectively. Finally, the torsion of the twisted pair is defined as α=2πℓe.

3. Twisted electrostatic model

To study the capacitive coupling between the twisted wires using the finite element method, we first derive a weak form of Maxwell’s equations in the electrostatic approximation.

3.1 Electrostatic model

Assuming zero free charge density, and denoting the (static) electric field and the dielectric permittivity in Cartesian coordinates (x₁, x₂, x₃), respectively, by e(x₁, x₂, x₃) and ε(x₁, x₂, x₃), Maxwell’s equations reduce to:

(2) curlx(e)=0, divx(ε e)=0,

where curl_x and div_x are the classical curl and divergence operators in the Cartesian coordinates system, respectively. Defining the scalar electric potential v(x₁, x₂, x₃) such that e = −grad_x(v), (2) leads to the second-order equation as follows:

(3) −divx(ε gradx(v))=0,

where grad_x denotes the gradient operator in Cartesian coordinates.

Our aim is to solve (3) in a domain Ω_x enclosing the twisted pair, with a prescribed scalar potential on the (boundary of the) conductors. To use the finite element method, Ω_x is bounded and an appropriate geometrical transformation is defined in an annular region to handle the natural extension of fields to infinity (Henrotte et al., 1999; Remacle et al., 1994; Bossavit, 1998). The annular region (of inner and outer radii R_int and R_out, respectively) is depicted in Figure 3. The boundary ∂ Ω_x of Ω_x is thus composed of the boundary of the conductors in the twisted pair Γ_c and the outside boundary of the annular region Γ_out.

3.2 Weak formulation in Cartesian coordinates

Let us assume that the potential v is fixed to zero on Γ_out and to a prescribed value v_c(x₁, x₂, x₃) on Γ_c. The weak formulation of (3) writes (Ern and Guermond, 2013): find v ∈ H¹(Ω_x) with v = v_c on Γ_c and v = 0 on Γ_out, such that:

(4) ∫Ωxε gradx(v)·gradx(v′)dΩx=0

holds for all test functions v′∈H01(Ωx). Here, H¹(Ω_x) denotes the classical Sobolev space on Ω_x, and H01(Ωx) its subspace with functions vanishing on the boundary ∂ Ω_x. For the 3D geometry, periodic boundary conditions along the x₃ axis are used. After discretizing Ω_x with a finite element mesh, and choosing suitable finite dimensional subspaces (e.g. with piecewise linear or quadratic shape functions), (4) becomes:

(5) ∫Ωx(ε gradx(vh))T gradx(v′h) dΩx=0,

where the superscript “T” denotes the transpose, whereas v_h and v′h denote the finite dimensional approximation of v and v′, respectively, leading thus to the classical matrix system.

Following the approach described in Nicolet et al. (2006, 2007a), we now reformulate the weak formulation (5) in helicoidal coordinates.

3.3 Helicoidal coordinates

Helicoidal coordinates (ζ₁, ζ₂, ζ₃) are related to the Cartesian coordinates (x₁, x₂, x₃) through the following relations (Figure 4) (Waldron, 1958):

(6) {x1=ζ1cos⁡(αζ3)+ζ2sin⁡(αζ3)x2=−ζ1sin⁡(αζ3)+ζ2cos⁡(αζ3)x3=ζ3,

or, conversely:

(7) {ζ1=x1cos⁡(αx3)−x2sin⁡(αx3)ζ2=x1sin⁡(αx3)+x2cos⁡(αx3)ζ3=x3.

The Jacobian of the transformation from helicoidal to Cartesian coordinates is defined as follows:

(8) J=∂(x1,x2,x3)∂(ζ1,ζ2,ζ3)=(∂x1∂ζ1∂x1∂ζ2∂x1∂ζ3∂x2∂ζ1∂x2∂ζ2∂x2∂ζ3∂x3∂ζ1∂x3∂ζ2∂x3∂ζ3),

which, using (6), writes:

(9) J=(cos⁡(αζ3)sin⁡(αζ3)αζ2cos⁡(αζ3)−αζ1sin⁡(αζ3)−sin⁡(αζ3)cos⁡(αζ3)−αζ1cos⁡(αζ3)−αζ2sin⁡(αζ3)001).

If grad_ζ denotes the gradient in helicoidal coordinates and J^−T the inverse of the transpose of J, we have:

(10) gradx(vh)=J−Tgradζ(vh),

and the weak form (5) in helicoidal coordinates becomes (because dΩ_x = det(J)dΩ_ζ: find v_h such that:

(11) ∫Ωζ(ε J−T gradζ(vh))T J−Tgradζ(v′h) det⁡(J) dΩζ=0⇒∫Ωζε gradζ(vh)TJ−1 J−T gradζ(v′h)det⁡(J) dΩζ=0

where Ω_ζ refers to the domain of study in helicoidal coordinates. Equation (11) can be further simplified as follows:

(12) ⇒∫Ωζε gradζ(vh)T T−1 gradζ(v′h) dΩζ=0,

with T⁻¹ the symmetric matrix:

(13) T−1(ζ1,ζ2)=T-T(ζ1,ζ2)=J−1J-Tdet⁡(J)=(1+α2ζ22−α2ζ1ζ2−αζ2−α2ζ1ζ21+α2ζ12αζ1−αζ2αζ11),

thus making (12) become:

(14) ∫Ωζ(εT−1 gradζ(vh))T gradζ(v′h) dΩζ=0.

Formally, comparing (5) to (14), the dielectric material represented by the scalar permittivity ε in Cartesian coordinates is thus replaced in helicoidal coordinates by an equivalent anisotropic and inhomogeneous permittivity tensor density εT⁻¹ (Nicolet et al., 2007a, 2007b). The matrix T⁻¹(ζ₁,ζ₂) being independent of ζ₃, (14) can be reduced to a 2D problem in terms of coordinates ζ₁ and ζ₂, on a “slice” of the 3D geometry with constant ζ₃ = x₃, without information loss. When expressed in polar coordinates (ρ,ϕ) defined as follows:

(15) {ζ1=ρcos⁡(ϕ)ζ2=ρsin⁡(ϕ),

with ϕ=2arctan⁡(ζ2ζ1+ζ12+ζ22), T−1 can be rewritten as follows:

(16) T−1(ζ1,ζ2)=R(ϕ)(10001+α2ρ2αρ0αρ1)R(−ϕ),

where R(ϕ) is the rotation matrix:

(17) R(ϕ)=(cos⁡(ϕ)−sin⁡(ϕ)0sin⁡(ϕ)cos⁡(ϕ)0001).

As mentioned in Section 3.1, the domain of study is bounded by an annular region of inner and outer radii R_int and R_out, respectively, to handle the natural extension of fields to infinity. The matrix is thus piecewise defined, with (16) being valid for ρ ≤ R_int. The expression of the matrix for ρ ∈ ]R_int, R_out] is presented next.

3.4 Infinite transformation

In the annular region ρ ∈ ]R_int, R_out], a geometrical transformation is applied to handle the natural extension of the fields to infinity (Henrotte et al., 1999; Remacle et al., 1994), and the matrix T⁻¹ writes:

T−1(ζ1,ζ2)=R(ϕ)(rdrρ000(1+α2r2)drρrαdrr0αdrrdrrρ)R(−ϕ).

The parameters r and d_r are the transformed radial cylindrical coordinates and the derivative of r, respectively:

(18) {r=r(ρ)=(Rint−Rout)ρ−Routρdr=dr(ρ)=∂r∂ρ=−RoutRint−Rout(ρ−Rout)2.

The infinite transformation can be seen as a change in the physical properties of the “infinite air,” thus an anisotropic and piecewise defined tensorial permittivity is used in the 2D finite element implementation.

4. Verification of the two-dimensional twisted model

Given the spatial periodicity of the geometry, only one spatial period referred to as the 3D geometry, with its cross section defining the 2D geometry is studied. For all geometries, created and meshed using Gmsh (Geuzaine and Remacle, 2009), an “infinite” transformation is applied in an annular region as described in Section 3.4. The 3D geometry was used to solve (5) and the 2D to solve (14). The finite element computations were performed using the open source finite element solver GetDP (Dular et al., 1998a).

In real applications, enameled conductors are coated in successive thin layers to ensure correct polymerization and cross-linking. This defines the thermal and electrical performance of the wire. A classical arrangement is polyester-imide (PEI) as the first layer and polyamide-imide (PAI) on top. Additionally, for this case study the layers were given the same thickness. Thus, ε is defined as a piecewise function and its values were chosen for the two layers in 3D to be ε_PEI = 4.5ε₀ and ε_PAI = 2.4ε₀, where ε0=136π 10−9F m⁻¹ is the vacuum permittivity, 4.5 and 2.4 being the relative permittivities of PEI and PAI, respectively. The equivalent anisotropic and inhomogeneous dielectric permittivity for the 2D twisted study is hence easily deduced (Section 3.3). This is portrayed in Figure 5, where a capacitor is used to illustrate the capacitive coupling between the two conductors.

To verify the 2D twisted model, simulations were run where the quantity of interest is the charge per meter Q of the system (Dular et al., 1998b). Using finite elements, and for both geometries, the computation of the charges per meter is done as follows:

(19) {3D  case:  Q3D=1ℓe∫Γcε gradx(vh)·n dΓc2D  case:  Q2Dα=∫ΓcεT−1 gradζ(vh)·n dΓc,

where n is the outside normal to Γ_c, Q_3D is the charge computed in 3D converted to a 2D value by dividing by ℓ_e and Q2Dα is the charge computed by the 2D twisted model.

Finite element computations were performed for a predefined range of N_turns per meter varying from 64 to 768, for a fixed value of ℓ_s. In all computations the potential difference between the two conductors was set to 1 V, and the potential at infinity was set to 0 V. The results are presented in Figure 6, where the corresponding 3D geometries were added to illustrate the amount of twist. The relative error used for comparison is defined as follows:

error(%)=100.|Q3D−Q2Dα|Q3D.

As depicted in Figure 6(b), the relative error does not exceed 0.06%, which validates the twisted model. (The remaining small error is because of the small discrepancy between the 3D mesh and the exact twisted geometry.) The computational efficiency of the 3D and 2D models is reported in Table 1. The number of nodes in the 3D mesh is about 21 times the number of nodes in 2D mesh, which explains the large difference in the computational times: 16 min for the 3D model and only 8 s for the 2D model.

As noticed in Figure 6(a), the increase of N_turns yields an increase in the computed charge, which is because of the increase of the capacitive coupling between the different parts of the studied conductor. This extra capacity, referred to as C_s, is illustrated in Figure 7. The higher N_turns, the less air there is between the two electrodes of the equivalent capacitance C_s, the higher is the computed charge. It is also worth mentioning that for very high and non-physical values of N_turns, the relative error starts to increase, because of the increased deformation of the 3D mesh. Finally, it should also be noted that in addition to global quantities like the charge, the 2D model also allows to retrieve local quantities like the electric field strength.

5. Impact of the twisted geometry on the computed values

To emphasize the importance of using the proposed 2D twisted model instead of simplifying the study to a 2D straight model, we study the influence of the geometrical parameters N_turns, d and, ℓ_i over the computed charges Q2Dα and Q2D0. The latter is the computed 2D charge for a straight geometry, i.e. with α = 0.

5.1 Impact of the number of turns

To study the impact of the number of turns N_turns over the computed 2D charges, d was fixed to 0.8 mm, and ℓ_i to 28 µm. A new variable Q2Dαc is introduced corresponding to a corrected 2D twisted 2D charge deduced from a straight charge analytically, as follows:

(20) Q2Dαc=Q2D0ℓhℓe=Q2D0(π(d+2ℓi)ℓe)2+1,

where ℓ_e and ℓ_h are defined in Figure 2. The results are plotted in Figure 8, where error₁ and error₂ are computed as in (21):

(21) {error1(%)=100.|Q2Dα−Q2D0|Q2Dαerror2(%)=100.|Q2Dα−Q2Dαc|Q2Dα.

The higher the value of N_turns, the higher Q2Dα, whereas Q2D0 remains constant because no twist is geometrically introduced. For N_turns per meter equal to 64, the relative error is roughly 1.0549% and increases with N_turns because of the emergence of new capacitive coupling as elaborated in the previous section, increasing thus the discrepancy between the two.

The engineering approach adopted to correct the 2D straight values Q2D0 to twisted ones Q2Dαc thus gives approximate, but not precise, values of the twisted charges. Therefore, the ratio of lengths is not enough to compute the value of the twisted charge. This proves the added value of the 2D twisted model and the importance of its usage when it comes to a twisted geometry. Indeed, the proposed 2D twisted model has the same computational efficiency as the 2D straight model, does not need a post processing (i.e. correcting the straight values with the length ratio) and fully considers the impact of the twist over the computed values.

5.2 Impact of the diameter of the conductor

To study the impact of the diameter of the conductor, N_turns per meter was set to 320 turns and ℓ_i to 28 µm. The diameter was varied from 0.4 to 5.16 mm. The results are seen in Figure 9.

As expected, the charges, straight and twisted, increase with the diameter of the conductor because the surface of the copper increases as well. This is illustrated in Figure. 5 and even though the two conductors are not planar, the formula for plane capacitors can be used to understand the impact:

(22) Cpp=εAℓi,

where A is the surface of the parallel copper plates, which increases with d. Moreover, the two charges Q2Dα and Q2D0 do not increase in the same manner. This is seen in Figure 7, where the twisted geometry introduces capacitive coupling between the different parts of the twisted conductors (C_s): the higher α, the lesser the thickness of the equivalent insulation of C_s, the higher is the value of the computed charge. The values of Q2Dαc were also computed and plotted in Figure 9, where the discrepancy between Q2Dαc and Q2Dα is because of the increase of the values of d, which increases the value of the length ratio (ℓ_h increases with d, whereas ℓ_e is constant because N_turns is constant).

5.3 Impact of the insulation thickness

Here, d was set to 0.8 mm and N_turns per meter to 320. The insulation thickness ℓ_i was varied from 10 to 50 µm. The results are seen in Figure 10 where a 2D geometry is presented to visualize the impact.

Following the same explanatory approach used in Section 5.2, and by making use of the parallel plate capacitive coupling analogy, the results can be easily interpreted as well. The capacitance (or the charge because ΔV = 1 V) is inversely proportional to the thickness of the dielectric. Thus, the smaller ℓ_i, the easier the storage of the charges. For this study also, the values of Q2Dαc were computed and seem to converge to Q2Dα for very high non-physical values of ℓ_i.

5.4 Study of the standard IEC 60851-5

As an application to the developed model, we propose to investigate the twist configurations in the IEC 60581-5 standard (IEC, 2019), which are summarized in Table 2.

For each diameter range d, a value of N_turns (ultimately α) is defined for ℓ_s = 0.125 m, which is also fixed by the standard. The higher the value of d, the lesser the N_turns and the less twisted the geometry is. The performed study consists in computing for each d varying from 0.1 to 2.4 mm, ultimately the upper and lower bounds of the specified range, the 2D charges: Q2Dα and Q2D0. Their variations according to d are plotted in Figure 11. The relative error was computed as in (21) (i.e. error₁). Here, the maximum number of turns N_turns per meter does not exceed 264 turn per meter (value computed for N_turns = 33) which is inferior to the values chosen in Sections 5.1–5.3.

In Figure 11, two parameters vary: the diameter d and the number of turns N_turns. To understand the variations, each sub-range where only one parameter varies will be studied separately. Taking, for instance, the sub-range of d where N_turns = 4, corresponding to a diameter varying between 1.4 and 2 mm, increasing d increases simultaneously the charges ( Q2Dα, Q2D0) as well as the relative error between the two. This agrees well with the results of Section 5.2. Moreover, the variation of the slope of the relative error is because of the decrease of N_turns with the increase of the value of d. The higher the d is, the less twisted the geometry is, the smaller the value of C_s, the slighter is the difference between the twisted and straight charges until it eventually nullifies for a N_turns = 0. It should also be mentioned that in real applications, the insulation thickness ℓ_i also varies with the diameter d (IEC, 2013) which was considered constant for the study of the standard.

6. Conclusion

This paper presents an efficient method to solve a 3D twisted electrostatic problem using 2D finite elements based on a helicoidal change of coordinates. The object of study is a twisted pair widely used in different electrical engineering applications. The low time complexity of the 2D model is highlighted, as well as its high accuracy, when compared to the 3D model. The importance of using a 2D twisted geometry instead of a 2D straight geometry to study a 3D twisted problem is emphasized.

The proposed model, by inserting the correct permittivities of the dielectric materials used in real application, could efficiently provide the value of the capacitance between twisted wires. The model may even be used as an alternative to measurements, enabling the analysis of twisted wires only by simulations. The model could also be adopted for various other case studies of twisted pairs: a Debye model could, for example, be used to describe the frequency-dependent behavior of the used insulation; or the impact of the temperature over the permittivity that could also be taken into account by using, e.g. an Arrhenius model.

Figures

Figure 1.

An example of a twisted pair with N_turns turns over a swirl length ℓ_s

Figure 2.

(a) One spatial period of the twisted pair. (b) The geometrical parameters of the cross-section

Figure 3.

Bounded domain with annular region with geometrical transformation to handle the natural extension of fields to infinity

Figure 4.

Definition of Cartesian and helicoidal coordinates

Figure 5.

Illustration of the capacitive coupling between the two conductors

Figure 6.

(a) Comparison of the computed charges per meter for the 2D and 3D geometries. (b) The relative error between Q_3D and Q2Dα

Figure 7.

Figure portraying the new capacitive couplings C_s a twisted geometry presents

Figure 8.

(a) Study of the impact of the number of turns per meter on Q2Dα, Q2D0 and Q2Dαc. (b) The computed relative error between Q2Dα and Q2D0. (c) The computed relative error between Q2Dα and Q2Dαc

Figure 9.

Study of the impact of the diameter of the conductor on the computed charges

Figure 10.

Study of the impact of the insulation thickness on the computed charges

Figure 11.

Study of the standard

Table 1.

Comparison of the computational efficiencies of the 2D and 3D models

	2D	3D
No. of nodes	47,755	1,002,945
No. of elements	101,188	2,226,182
Computational time	8 s	16 min

Table 2.

IEC 60581-5 Standard

Range of the diameter (mm)	Nturns	Nturns per meter
0.100 ≤ d < 0.250	33	246
0.250 ≤ d < 0.355	23	184
0.355 ≤ d < 0.500	16	128
0.500 ≤ d < 0.710	12	96
0.710 ≤ d < 1.060	8	64
1.060 ≤ d < 1.400	6	48
1.400 ≤ d < 2.000	4	32
2.000 ≤ d < 2.500	2	16

Source: IEC (2019)