Domain decomposition with subdomain pre-processing for finite element modelling of transformers with stranded conductors

Finite element analysis of transformers with stranded conductors can be computationally heavy due to the fine mesh required for winding domains. In this paper, a domain decomposition method with subdomain pre-processing is proposed to overcome with this problem. In order to solve the field-circuit problem with finite element analysis, the winding domain is modelled initially. Next, the reduced matrix system is obtained from the winding domain formulation. The solution of the entire problem is obtained by solving the winding domain formulation and the reduced matrix system separately. The method is tested on a transformer model in no load operation. Accurate solutions are obtained in both frequency and time domains 7-9 times faster than when using brute-force method. © 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5130076., s


I. INTRODUCTION
In order to reduce AC-losses, the windings of transformers and electrical machines are sometimes composed of several stranded conductors. Using stranded conductors instead of big solid conductors reduces the eddy current losses. Furthermore, it makes the manufacturing process easier and selection of conductors more flexible. On the other hand, it causes circulating current losses, which are the result of the uneven distribution of the current among the parallel strands. For accurate prediction of winding losses, a fine mesh inside the winding domains is needed. However, using a fine mesh with finite element (FE) analysis may result in very long computation time. Therefore, homogenization techniques and models that reduce computation time are required.
In this paper, a FE modelling method is proposed by decomposing domains with subdomain pre-processing to analyze the proximity effect and circulating current losses in transformers with stranded conductors with reduced computation time. The method is heavily inspired by the similar approach proposed earlier for harmonic analysis of electrical machines with multi-conductor windings. 23 The proposed method in this paper is applicable for both frequency and time domain analyses of the transformer. It allows analyzing different connections of winding domains.
A transformer composed of a linear core and 120 stranded conductors in winding domains have been studied. The method is simulated in both frequency and time domains with different connections of winding domains. Sensitivity analysis of the method for different mesh refinements and different numbers of boundary nodes have been made. Also, the simulation results of the proposed method and brute-force method have been compared. According to the results and comparisons, high accuracy has been obtained with considerably shorter computation time.

II. REDUCED BASIS MODELLING
In this section, the basic theory for solving a 2D field-circuit FE problem will be presented first. Afterwards, the proposed method will be explained in detail.

A. Basic theory
In this paper, the well-known A-V formulation is used to solve the field-circuit problem. 25 The problem consists of nodal vector potentials a, voltages over the conductors u, and a set of linearly independent loop currents i. In time domain, the matrix system with A-V formulation can be expressed as where S and M are well-known stiffness and mass matrices respectively. C J and C E matrices describe the current density in the conductors created by the voltage term and induced back-emfs on the conductors respectively. The entries of C J and C E are where φi is the shape function of the node i in domain Ωj c , σ is the conductivity, and le is the axial length. Moreover, in the matrix system given by (1), R is the diagonal matrix that consists of DC-resistances of each conductor, Z is the matrix related to end-winding impedances, which is taken as 0 in this paper, and L is the loop matrix denoting the connection of the windings, with the entries Finally, us is sum of supply voltages for each current loop.

B. Description of the proposed method
In brute-force method, the matrix system in (1) is solved to find the solutions of a, u, and i. As the size of the matrix system is large, long computation time is required.
In proposed method, two different domains are defined as nonwinding and winding domains. Also, boundary is defined between the domains. By separating the nodes of domains, the matrix system in (1) can be expanded as Instead of solving the A-V formulation directly, the vector potential of the nodes of winding domains aw and the voltages of the winding conductors u are expressed with a matrix system, winding domain formulation, in terms of the vector potential of the nodes of non-winding domain an, those of its boundary a b , and the currents flowing through individual conductors i. Then, the expression obtained for the winding domains are inserted into the other rows of A-V formulation to form reduced matrix system. From the reduced matrix system, an, a b , and i are solved first. Afterwards, aw and u are solved from the winding domain formulation. In the following parts, the procedure mentioned above will be explained in detail.

C. Winding domain formulation
From the third and fourth row of the matrix system in (5), aw and u can be expressed as where aw and u correspond to entire winding domains. If the winding domains have identical geometry and mesh, vector potential and conductor voltages of one winding domain can be used since the expressions will be the same for other winding domains. The inverse matrix in (6) will be defined as Q matrix, and Q11, Q12, Q21, and Q22 entries will be used for the remainder of the paper. Therefore, aw and u can be expressed as

D. Reduced matrix derivation
Next, the winding domain formulation is inserted into the remaining rows of the matrix system in (5) in order to form the reduced matrix system. The reduced matrix system can be expressed as Therefore, an, a b , and i can be solved from the reduced matrix system.

E. Solution procedure for harmonic analysis
In harmonic analysis, d dt terms in (5) and (6) are replaced with jω, where ω is angular frequency. The winding domain formulation and reduced matrix system are obtained in pre-computation stage and stored. Since linear core case is studied in this paper, stiffness AIP Advances 10, 015133 (2020); doi: 10.1063/1.5130076 10, 015133-2 ARTICLE scitation.org/journal/adv matrix terms can also be stored as they are constant. For non-linear problems, the Newton-Raphson method can be used, which does not change the main logic of the method. After obtaining and storing winding domain formulation and the reduced matrix system in pre-computation stage, an, a b , and i terms are solved from the reduced matrix system. Next, aw and u terms are solved from the winding domain formulation by using the solutions of an, a b , and i. The flow diagram of the harmonic analysis is shown in Fig. 1.

F. Solution procedure for time-stepping analysis
In time-stepping analysis, the d dt terms in (5) and (6) are timediscretized according to the preferred approach. In the simulations of this paper, Implicit Euler approach has been used.
The winding domain formulation and reduced matrix system are obtained and stored in pre-computation stage. After pre-computation, the time-stepping analysis is initiated. In timestepping, firstly an, a b , and i terms are solved from reduced matrix system. Next, aw and u are solved from winding domain formulation by using the solutions of an, a b , and i. Afterwards, the load vector is updated and next time-stepping can be performed similarly. The flow diagram of the time-stepping analysis is shown in Fig. 2.

G. Non-conforming boundary case
The proposed method is applicable also when the meshes of the non-winding and winding domains are non-conforming. In that case, two different sets of boundary nodes a b1 and a b2 are defined. These boundary nodes are interpolated linearly, and an interpolation matrix P is obtained such that a b2 = Pa b1 . After obtaining the continuity between the boundaries, the vector potentials can be expressed as which can be simplified to a = Ptotanew.
Afterwards, a b2 can be eliminated by inserting (12) into Sa = F expression and multiplying both sides with P T tot as Then, the problem becomes suitable for the application of the solution procedures explained in the previous parts.

III. SIMULATIONS
The method has been simulated for a transformer in no load operation consisting of 120 identical stranded conductors. in parallel with each other. However, the method is applicable with any type of connections. A mesh consisting of 9177 nodes in total; 411, 30, and 8736 of which belongs to the non-winding domain, boundary, and winding domains respectively, has been constructed and this case has been chosen as reference case. The mesh of the transformer and one winding domain is shown in Fig. 3. Simulations have been performed by using both proposed method and brute-force method. The field-circuit problem has been solved with harmonic and time-stepping analyses.
After solving the field-circuit problem, flux density and current density plots have been obtained. The accuracy of the proposed method has been validated by comparing the results with brute-force method. Error norms for the vector potentials have been found. Computation times of the methods have been compared for both analyses. Finally, sensitivity analyses have been performed with 6 different mesh structures. The effect of mesh refinement in the efficiency of the proposed method has been investigated.

Harmonic analysis
Firstly, harmonic analysis has been performed to solve the fieldcircuit problem. 140 V, 50 Hz has been applied as a supply voltage. The flux density distribution of the transformer and the current density distribution of the conductors obtained from the simulations are shown in Fig. 4a and Fig. 4b. Also, comparison of the computation times of the methods in harmonic analysis can be found in Table I. In the proposed method, pre-computation takes 0.36 s. The pre-computation stage can be stored once and re-used for different harmonic analyses. The time spent in harmonic analysis of the proposed method is 9 times smaller than the that of brute-force method. If the harmonic analysis is repeated several times, the precomputation time becomes negligible and the speed of the proposed method becomes more significant.

Time-stepping analysis
Secondly, time-stepping analysis has been performed to solve the field-circuit problem. 17.5 V has been applied as a sinusoidal voltage. 2 periods have been analyzed with 100 time-stepping per period. Comparison of the computation times of the methods in time-stepping analysis can be found in Table II. According to Table II, the time-stepping analysis of the proposed method is 9 times faster than the that of brute-force method. The pre-computation is negligible in this case. Error norms of vector potentials in last period can be found in Fig. 5.

B. Sensitivity analysis
In order to investigate the effect of the mesh refinement of the winding domains on the computation time, different mesh refinements have been simulated. Initially, one coarser and one finer mesh compared to the reference case have been investigated. Then, the number of boundary nodes in the winding domains has been changed to 18 and its effect have been analyzed with 3 different mesh refinements inside the winding domains. The computation times of each mesh refinements for 10 and 18 boundary nodes cases can be found in Table III and Table IV respectively.
According to the results of each case, it is observed that when the mesh of the winding domains is coarser, the proposed method becomes faster. The results of different mesh refinements are consistent with each other.
When the cases are compared with each other, it can be observed that as the number of boundary nodes increases from 10 to 18 for each winding domain, the time spent on pre-computation

IV. CONCLUSIONS AND DISCUSSIONS
In this paper, a domain decomposition method with subdomain pre-processing has been proposed to observe winding losses from finite element analysis while decreasing the computation time. Firstly, the voltage across the conductors u and vector potentials of winding domains aw have been expressed in terms of the vector potentials of non-winding domain an and boundary a b , and the currents of the individual conductors i in order to obtain reduced matrix system form. Reduced matrix system and winding domain formulation have been solved separately to obtain the entire solution.
The proposed method has been tested with a transformer model in no load operation consisting of a linear core in nonwinding domain and 120 stranded conductors in winding domains. The results of the harmonic and time-stepping analyses show that the proposed method is 7-9 times faster than the brute-force method. As the mesh refinement of winding domains is coarser, the speed increases relatively. Also, if the number of boundary nodes of a winding domain is less, the proposed method becomes more time efficient.
According to the results of the simulations, the proposed method shows good accuracy compared with the brute-force method in 2-D finite element analysis. As a future work, this method should be implemented for a more realistic transformer, and the effect of eddy and circulating currents should be observed for different frequencies.