3.2 Thermal Model II
3.2.5 Thermal anisotropic properties of windings
The full 3-D numerical model of electric motor presented in the following sections con-tains homogeneous solid structures representing coils wounded on each stator core tooth.
In reality, windings consist of a copper conductor, its insulation and air gaps between neigh-bouring wires. Such an approach applied to the windings forces the implementation of ef-fective thermal conductivity in all the directions. However, the crucial values need to be determined in the two specific directions across the wires: transverse and axial directions with respect to the motor shaft. Due to the complex shape of the windings and three compo-nents materials, i.e. copper conductor, wire insulation and air between coils, characterised by different thermal conductivity, the effective thermal conductivity in the specific direction
Bobbin
Air gap area Winding
Diameter wire
Equivalent air gap thickness
Isolation
Figure 3.11: Scheme of the wire alignment to the bobbin surface
was calculated using separate 2-D CFD models and analytic equations.
To calculate the effective thermal conductivity of the windings in the transverse direc-tion, 2-D CFD model was prepared. In this secdirec-tion, the model description with the obtained results is presented.
In the winding cross-section model, the numerical domain was limited by four surfaces.
In this part of the study, two domain of two different windings arrangement (crossed and straight patterns) were investigated. These two arrangements are presented in Fig. 3.12, where the temperature-dependent thermal conductivity of the copper and the air is also pre-sented. Temperature difference in this case between boundary conditions was set as 5 K and thus the presented thermal conductivity values do not vary significantly. This temperature difference between two sides of the windings is expected in the real motor work conditions and that was confirmed during the experimental part. Within the presented domains of the winding cross-section, a number of 65 wires was used.
The 2-D numerical meshes for each winding arrangement used for calculation consisted
of approx. 3 million quadrilateral elements chosen after the mesh sensitivity analyses. The mesh sensitivity analysis was conducted for the selected case for the overestimated differ-ence of temperature between two domain surfaces at the level of 100 K. The presented results of the total heat transfer rate and the calculated effective transverse thermal conductivity of the winding are presented in Fig. 3.13. The finer mesh tested for the selected case achieved 20 million elements and allowed to obtain similar results as those presented in that figure.
· ·
Figure 3.12: Crossed and straight pattern of the winding wire alignment
In the numerical investigation, the effective thermal conductivity of the windings in the transverse direction mainly depended on the thermal conductivity of the coils insulation and the contact area between specific coils. The two neighbouring wires within the wind-ings are schematically presented in Fig. 3.14, where the conductor was pointed out using orange colour, the conductor insulation using brown colour and blue colour was used to de-note the air gaps. In Fig. 3.14, the conductor diameter was indicated as dimension A which was approx. 0.0007 m, while a dimension of the contact area between coils was depicted as the dimension B and was varied in a model as an input data. The contact area between coils could also be expressed as a ratio between this contact length and the constant conductor diameter. This ratio was treated as the input parameter in the following results presentation.
Therefore, the ratio of the contact area to conductor diameter can be estimated as the re-lation of B to A in Fig. 3.14. The material data about the thermal conductivity of the wire insulation were not available. For this reason, it was decided to investigate its range from 0.14 W · m−1·K−1to 1.2 W · m−1·K−1.
1.061 1.071 1.081 1.091 1.101 1.111 1.121
310 312 314 316 318 320 322 324 326 328
50 000 500 000 5 000 000
Total heat transfer rate Effective thermal conductivity Total heat transfer rate, W
25 000 000 mesh elements Effective thermal conductivity, W/(m·K)
Figure 3.13: Mesh sensitivity analysis of the 2-D model for the windings thermal conductiv-ity estimation
A B
Figure 3.14: Characteristic dimensions of the windings: conductor diameter (dimension A) and contact area between coils (dimension B)
The 2-D models of the winding cross-section were built on the basis of the winding di-mensions and the winding component materials that were discussed previously in the cur-rent section. One of the numerical meshes, built on the basis of crossed pattern and the medium value of the contact area between coils, is presented in Fig. 3.15. On the left-hand
0.0035 m
0.0102 m
(a) (b)
(c) (d)
Figure 3.15: Exemplary numerical mesh used in 2-D model of the crossed pattern for the windings thermal conductivity estimation: (a) full mesh view, (b) zoomed vew of three neigh-bouring wires, (c) zoomed view of the air gap between three wires and (d) zoomed view of the insulation and contact surface between neighbouring wires
side of the figure, the whole domain mesh is presented with its global dimensions. The mesh in this view is not visible due to the high amount of the elements. Therefore, in the pre-sented figure, the zoomed views are introduced on the left-hand side. The following zooms are shown in the coloured frames pointing each zoom position in Fig.3.15(b) to (d).
The model scheme including the computational domain and boundary conditions for the determination of the effective thermal conductivity in the horizontal direction is pre-sented in Fig. 3.16. In this model configuration, the assumed boundary conditions were based on the Dirichlet type, i.e. constant temperature that was defined on the left and right sides of the domain, while the zero heat flux was prescribed to the top and bottom sides of the domain. For the rated operating point of the motor work, the temperature difference reached approx. 5 K, while the higher temperature was at the level of 95◦C. Therefore, these
temperature values were defined as the boundary conditions in the discussed model of the winding cross-section. In the fluid zone, both convection and radiation phenomena were taken into account. As a result of the presented numerical model, the heat flux was esti-mated in a particular direction. Then the effective thermal conductivity was calculated from the analytical heat conduction equation on the basis of the area-averaged heat flux, the dis-tance between the surfaces where the temperature boundary conditions were defined and the temperature difference at the boundary conditions.
T
LOWT
HIGHq
q=0
q=0
Figure 3.16: Computational domain and boundary condition types for the numerical model of the winding cross-section for the determination of the effective winding thermal conduc-tivity in the horizontal direction
The results for the effective thermal conductivity of the windings in the transverse direc-tion for two coils arrangements and in the funcdirec-tion of the discussed factors are presented in Fig. 3.17. Assuming that the thermal conductivity of the coil insulation is at the level of 0.22 W · m−1·K−1and proportion of contact area to conductor diameter achieved 10.7%, the resulting effective transverse conductivity is approx. 1.36 W · m−1·K−1for the windings wounded in the crossed pattern. In the case of the same parameters for the straight pat-tern, the model results in the windings effective transverse conductivity at the level of 0.95 W · m−1·K−1.
Taking into account that in the real geometry the winding is wounded on the iron in a mixed way, the winding effective thermal conductivity for the assumed parameters can be estimated as an average value of the results obtained from the crossed and straight pattern.
Therefore, the effective thermal conductivity value in transverse directions transferred to
the winding region in the 3-D model simulations was set to 1.15 W · m−1·K−1. The calcu-lated thermal conductivity is consistent with the values reported in the experimental and numerical studies [57,90,99].
3.2.6 Boundary conditions
In the previous section, the governing equations with an accompanying set of bound-ary conditions were introduced. These boundbound-ary conditions included three different types of given data on the domain boundary. On the external surfaces of the air block surround-ing motor, the constant pressure inlet with a value of 0 Pa was employed. At the domain boundary above the motor, the pressure outlet with value 0 Pa was defined. At the described constant pressure boundary conditions, the temperature of the room surfaces was defined in the numerical model. The emissivity value of the room walls was assumed as 1 due to their higher dimensions comparing to the motor size. The emissivity of all specific motor sur-faces was set individually at the same level as presented in Section3.1.6. The only difference was present in variants above Variants A and was related to the emissivity of the aluminium brackets on the level of 0.95. Moreover, the same emissivity value was set for all the radiators.
The surface of the aluminium base was treated as a wall boundary condition with zero heat flux. Between the two motors, the one investigated and the other working as a generator, a symmetry boundary condition was used on a surface normal to the axis of both motors.
More details on the role of the generator were already given in Chapter2.
0.75
Figure 3.17: Effective thermal conductivity of windings in transverse direction as a function of the insulation thermal conductivity and the ratio of the contact area to the conductor diameter: (a) for the crossed pattern and (b) for the straight pattern
53
3.3 2-D Electromagnetic model
The separate electromagnetic model was developed to simulate electromagnetic phe-nomena that occur within the electric motor. As it was mentioned, the main aim of this model was the power loss estimation independently from the experiment tests. The detailed 2-D model is discussed in Sections from3.3.1to3.3.4. To formulate the model, the Maxwell software was used for the EMAG simulations [100]. This solver is based on the FEM tech-nique. Moreover, the domain preparation and mesh construction for the EMAG analysis were conducted within the frame of the EMAG solver software, not separately as it was car-ried out in the case of the CFD model.
3.3.1 2-D EMAG domain
In the EMAG analysis, the numerical domain was simplified comparing to the one pre-sented in thermal models. The main components in the domain were as follows: stator core, windings, magnets, rotor core, shaft and separated surfaces of air. In the EMAG analysis, the PCB material, plastic elements, housing, aluminium fixing and other additional elements were omitted because they are not relevant from the electromagnetic point of view. The EMAG domain was depicted with the accompanying mesh in Fig. 3.2. The numerical mesh is described in the next section.
3.3.2 2-D EMAG mesh
Numerical mesh used in 2-D EMAG analysis consisted of over 28 thousand nodes. As it was mentioned, the EMAG analysis was based on the FEM technique. Therefore, the in-terpolation process between the EMAG and CFD solvers was performed as it is described in Section3.5. In EMAG model, the air gap between stator and rotor was one of the most crucial regions where very fine mesh was generated. Therefore, a number of 7 thousand elements were implemented in the air gap between stator teeth and rotor, named as band, which was treated as a rotating part. The mesh was built automatically preserving defined size of the specific elements. The detailed information about the numerical mesh used in 2-D EMAG analysis were presented in Tab. 3.2. Therefore, a crucial aspect of the presented mesh qual-ity is the distance between neighbouring nodes. In Tab. 3.2, besides the elements number in the specific zones also the other mesh parameters are listed. In the presented table, the
Table 3.2: Numerical mesh parameters used in the 2-D EMAG analysis
Element Num. of Min. Max. Min. elem. Max. elem. Mean elem.
elements length, mm length, mm area, mm2 area, mm2 area, mm2
iron rotor 615 0.3828 1.997 0.19345 1.664 0.692
iron stator 9214 0.1445 0.999 0.02687 0.414 0.175
magnet 1 431 0.0226 1.000 0.00083 0.345 0.071
magnet 2 416 0.0226 0.983 0.00108 0.370 0.074
magnet 3 431 0.0226 1.000 0.00083 0.345 0.071
magnet 4 421 0.0226 1.000 0.00083 0.345 0.073
magnet 5 441 0.0226 1.000 0.00083 0.345 0.070
magnet 6 421 0.0226 1.000 0.00083 0.345 0.073
magnet 7 431 0.0226 1.000 0.00083 0.345 0.071
magnet 8 441 0.0226 1.000 0.00083 0.345 0.070
shaft 284 0.3775 1.993 0.18443 1.447 0.575
band 7216 0.0293 0.180 0.00102 0.011 0.008
inner Region 2848 0.0077 1.075 0.00009 0.412 0.017
outer Region 3644 0.0358 2.892 0.00252 2.756 0.213
windings 1 320 0.2981 0.795 0.06413 0.227 0.134
windings 2 341 0.3274 0.790 0.06815 0.258 0.126
windings 3 320 0.2981 0.795 0.06413 0.227 0.134
windings 4 341 0.3274 0.790 0.06815 0.258 0.126
windings 5 320 0.2981 0.795 0.06413 0.227 0.134
windings 6 346 0.3274 0.799 0.06719 0.236 0.124
windings 7 320 0.2981 0.795 0.06413 0.227 0.134
windings 8 320 0.2981 0.795 0.06413 0.227 0.134
windings 9 341 0.3274 0.790 0.06815 0.258 0.126
windings 10 320 0.2981 0.795 0.06413 0.227 0.134
windings 11 346 0.3274 0.799 0.06719 0.236 0.124
windings 12 341 0.3274 0.790 0.06815 0.258 0.126
minimum and maximum distances between the nodes are pointed out for specific elements.
Moreover, in the mentioned table, the maximum, minimum and mean areas between neigh-bouring nodes are presented for specific regions. The scheme of the numerical mesh is also presented in Fig.3.18. As one can see in the presented table and figure describing numerical mesh, the finest mesh was built in the air gap region between magnets and stator core. This region is crucial for the induction field being the calculation results. Therefore, in Fig. 3.18 (b) the air gap area is present where at least 5 nodes between magnets and stator core were introduced.
(a)
(b)
windings
stator magnets core
rotor core
air gap
Figure 3.18: Numerical mesh in the 2-D EMAG analysis: (a) the general mesh overview (b) zoom to the air gap
3.3.3 2-D EMAG governing equations
In the transient EMAG analysis, the governing equations derived from on the three Maxwell equations (Eqs (21)–(24)) according to the solver documentation [100] and according to the standard literature [101].
∇∇∇ × H =σE (21)
∇∇∇ × E = −∂B
∂t (22)
∇∇
∇ · B = 0 (23)
B = ∇∇∇ × A (24)
where H denotes the magnetic field,σ represents the conductivity, E is the electric field, B represents the magnetic flux density and A is the magnetic vector potential. The relation between the magnetic field and the magnetic flux inside the domain is given by Eq. (25)
B = 1
ϕH (25)
whereϕ is the reluctivity of medium. The documentation of the applied solver [100] pointed out that the time-dependent magnetic equation for 2-D problem can be written as in Eq.
(26).
∇∇∇ ×υ∇∇∇ × A = Js−σ∂A
∂t −σ∇∇∇ϕ + ∇∇∇ × HC+σv × ∇∇∇ × A (26) where: υ is the reluctivity, Js is the source current density, σ is the electric conductivity of material,ϕ is the electric scalar potential, HC is the coercivity of the permanent magnet, v is the linear velocity vector.
The motion equation is changed comparing to Eq. (26) because the moving components have been fixed to their own coordinate system and the partial time derivative is the total time derivative of A. Therefore, the motion equation can be expressed by Eq. (27) and is solved at each time step at every node in the finite element model. The motion is contained
implicitly in the total derivative of A.
∇∇
∇ ×υ∇∇∇ × A = Js−σ∂A
∂t −σ∇∇∇ϕ + ∇∇∇ × HC (27) In the electromagnetic analysis, similarly as in the thermal models, the motor winding at each tooth is represented geometrically as uniform solid material. The stranded conductor characteristic is assigned to this component. Therefore, the current density is averaged over the area of the analysed stranded conductor. The filaments of the stranded conductor can be connected in series or parallel. Furthermore, when analysis eddy current effects are not being considered as well as electric potentials, Eq. (27) simplifies to Eq. (28).
∇∇∇ ×υ∇∇∇ × A = Js (28)
The current in Eq. (28) was calculated from uniformly distributed current density from Eq. (29).
jS=df Nf ·if
Sf ·a · p (29)
where if is the total current flowing within a filament coil group, Nf is the conductor number in the windings, a is the number of parallel branches in the windings, df is the polarity to represent the current direction, Sf is the total area of the cross-section of the solid region representing windings, p is the ratio of the original full model to the field domain to be solved.
The essential aspect for the thermal analysis is the calculation of losses occurring dur-ing the considered motor work. The Joule heatdur-ing losses caused by the current distribution in the conductor can be expressed by Eq. (30). The important aspect of the EMAG analy-sis is estimation of the windings conductivity. The material of the windings was defined as annealed copper corrected by the filling factor. This conductivity was implemented as the value dependent on temperature.
∆PCu= |Js|21
σ (30)
The main losses occurring in the core material can be divided into the hysteresis (∆Ph), eddy current (∆Pc) and excess losses (∆Pe) [102]. The calculation of the core losses, also
named as iron losses, was proposed by Bertotti in [102, 103] and can be expressed by Eq.
(31) under the sinusoidal flux according to the software documentation [100] and standard literature [104]. The estimation of K coefficients in Eq. (31) should be conducted for the specific core material. The influence of the Kecoefficient on the core losses was investigated experimentally and numerically in the work of [105]. This coefficient takes into account the complex phenomena nonlinear electromagnetic field diffusion in the lamination [106,107].
It also depends on the type of the core production process. Zhang et al. [108] described two approaches while the first one negligibly the excess losses and the second approach taking it into account. In the presented model, the excess losses were analysed.
∆Pi r on= ∆Ph+ ∆Pc+ ∆Pe=Khf Bm2 +Kcf2Bm2 +Kef1.5Bm1.5 (31) where ∆Ph is the hysteresis losses, ∆Pc represents eddy current losses and ∆Pe depicts the excess losses, Bm is the amplitude of the AC flux component, f is the frequency, Kh is the hysteresis loss coefficient with a value of 320, Kcis the eddy-current loss coefficient reaching 1.3 and Keis the excess loss coefficient at the level of 1.72.
3.3.4 Material properties
In the EMAG simulations, the air was treated as the default with the same electromag-netic properties as a vacuum. Therefore, it was treated as non-conductive material with rel-ative magnetic permeability at the same level as vacuum reaching a value of approx. 1.256
·10−6H · m−1. The permeability of a core was implemented according to B-H curve for the estimated core material as presented in Fig. 3.19. The applied permanent magnets made of neodymium material were characterised by magnetic coercivity at the level 870 kA·m−1. The conductor material was characterised as annealed copper whose bulk electric conductivity was corrected by the ratio of filling factor reaching approx. 0.75. Moreover, the conductor conductivity was implemented as temperature-dependent value according to the function presented in Fig. 3.20. This function was converted on the basis of the resistivity function which was temperature-dependent and it has been already expressed by Eq. (13). The dis-cussed material properties were used to solve the governing equations presented in Section 3.3.3.
. . .
B, T
H, A/m
Figure 3.19: Simplified scheme of the auxiliary model for a simulation of the current wave-form control system
Conductivity, MS/m
Temperature, °C Figure 3.20: Conductivity of annealed copper dependent on temperature
3.3.5 0-D auxiliary model for current regulation system
Modelling of the current waveforms that occur in the PM BLDC motor required formula-tion of a separate 0-D auxiliary model of the control system to simulate the work of the elec-tronic commutator. The auxiliary model was defined in the software based on the Simplorer
scheme that is a part of the EMAG solver. The scheme of the auxiliary system is presented in Fig. 3.21. In this figure, the voltage supply is represented by DC source on the left-hand side and was defined for each analysed load separately. The diodes presented on the scheme defined the current flow direction and voltage switches mimicked the work of transistors applied in the commutator system. The voltage switches were controlled by voltmeters pre-sented on the left-bottom part of the scheme where voltage pulse signals were initiated after achieving a proper angle of the rotor position. The voltmeters simulated the work of hall effect sensors used in the system of the rotor position control. In the EMAG model, the rotor position is calculated on the basis of the actual rotor angle and the starting rotor position. Ac-cording to Fig.3.21, three-phase wires supplied the phase windings representing the triangle connection. A challenging task during EMAG simulations was the definition of proper val-ues of inductances and resistances for the presented phases. It was achieved as a part of the experimental work using advanced electromagnetic analysers and their results calibration.
The windings connection including the number of parallel branches within this connection was identified by conducting the destructive examination of the investigated motor. In Fig.
The windings connection including the number of parallel branches within this connection was identified by conducting the destructive examination of the investigated motor. In Fig.