A Coupled , Semi-Numerical Model for Thermal Analysis of Medium Frequency Transformer

Medium-frequency (MF) transformer has gained much popularity in power conversion systems. Temperature control is a paramount concern, as the unexpected high temperature declines the safety and life expectancy of transformer. The scrutiny of losses and thermal-fluid behavior are thereby critical for the design of MF transformers. This paper proposes a coupled, semi-numerical model for electromagnetic and thermal-fluid analysis of MF oil natural air natural (ONAN) transformer. An analytical model that is based on spatial distribution of flux density and AC factor is exploited to calculate the system losses, while the thermal-hydraulic behavior is modelled numerically leveraging the computational fluid dynamics (CFD) method. A close-loop iterative framework is formulated by coupling the analytical model-based electromagnetic analysis and CFD-based thermal-fluid analysis to address the temperature dependence. Experiments are performed on two transformer prototypes with different conductor types and physical geometries for validation purpose. Results suggest that the proposed model can accurately model the AC effects, losses, and the temperature rises at different system components. The proposed model is computationally more efficient than the full numerical method but it reserves accurate thermal-hydraulic characterization, thus it is promising for engineering utilization.


Introduction
Power transformer is a critical component in electricity distribution systems.Nowadays, the medium-frequency (MF) transformer has become an important part of many power conversion systems.The high frequency reduces the volume of magnetics but risks increasing the core and winding loss densities, which cause deteriorated thermal performance and potentially lead to challenges, such as insulation damage, shortened lifespan, and even malfunctions or explosions [1].In this regard, the thermal limit of temperature rise on structural parts is the major concern restricting the design of MF transformer [2].
Due to the generated heat, temperature will increase largely if the transformer is left thermally isolated.Depending on the amount of heat to be dissipated, different cooling systems are used to maintain the temperature at an expected level.Amongst others, the oil natural air natural (ONAN) approach cools the device internally and externally with the natural movement of oil and air driven by complexity and occupies too much computing resources, thus it is not favorable in a design program that commonly takes several iterations to achieve the final solution.Hence, a trade-off has to be made to achieve sufficient modelling accuracy while keeping the computing cost at an expected level as well.An attempt was shown in [35], where the dimensionless least squares and upwind FEM were combined to simulate the thermal-fluid field in an oil-immersed transformer to improve the computing efficiency of coupling analysis.
In this paper, a coupled, semi-numerical model is proposed for the electromagnetic-thermal-fluid analysis of MF ONAN transformers.An analytical model that is based on spatial distribution of flux density and AC resistance factor is exploited for electromagnetic analysis to calculate the core and winding losses.A 3D numerical model that is based on CFD technique is then used in conjunction with the analytical model to predict the temperature rises of different system components.The electromagnetic and the thermal-fluid behaviors are solved iteratively in a close-loop manner to address the temperature-dependent properties and achieve a converged solution.Load experiments on two transformer prototypes with different conductor types and physical geometries are performed to validate the proposed method.The proposed method contributes to improving the modeling of the MF transformer in the following points.Firstly, the refined core loss modeling well addresses the flux inhomogeneity stemming from the nonlinearity of magnetic materials and the difference of magnetic paths, thus the core loss can be estimated with higher fidelity.Secondly, the AC effects are modeled in a multi-layer manner in seeking to determine the total winding loss more prudently.Thirdly, the proposed semi-numerical method avoids the time-consuming 3D electromagnetic numerical calculation, while it keeps a detailed modeling of both fluid flow and temperature distribution, thus it can be expected to better manage the trade-off between accuracy and complexity.Due to the improved electro-magnetic modeling, the proposed method also has good potential to be further simplified by replacing the numerical thermal model with the analytical model, so as to ease the real-time application.
The rest of paper is organized, as follows.Section 2 introduces the fundamentals of electromagnetic analysis for core and winding loss determination.Section 3 presents the 3D CFD model and the overall iterative framework of the proposed method.Section 4 describes the experimental setup and tests.The validation and discussion are presented in Section 5, while Section 6 draws the key conclusions.

Electromagnetic Modeling
The accurate estimation of losses, which depends on high-fidelity electromagnetic modelling, is critical, especially in the thermal point of view.Instead of using the complicated numerical method, this section describes an analytical method for the determination of core loss and winding loss, which are respectively based on the spatial distribution of flux density and AC factor.

Core Loss
The total core loss is composed of hysteresis loss, eddy current loss, and residual loss.At medium/ high frequency operations, amorphous and nanocrystalline materials exhibit lower hysteresis loss and they are commonly employed in medium/high frequency transformer applications.While separating the core loss into three sources needs extensive efforts on experiments and coefficient extraction, the empirical method requires much less measurements, thus it is more attractive in real applications.
Depending on the flux density, the core loss in per unit volume occurred could be approximated by the well-known empirical Original Steinmetz Equation (OSE) for purely sinusoidal excitation [36]: where K, α, and β are Steinmetz constants that depends on the material properties, f is the operating frequency in Hz, and B is the flux density.The sinusoidal excitation is applied on the low-voltage (LV) side in this study, which well fits to the assumption in OSE model.In the case of arbitrarily shaped excitations, several models have been developed for core loss modeling, including the modified It should be noted the flux density in Equation ( 1) is assumed to be the average of flux inside the core.Due to the nonlinearity of magnetic materials and different magnetic paths, however, the flux is not distributed homogenously depending on the material characteristics and core configuration.
E-E core is used in this paper, and half of the core geometry is shown in Figure 1.The amorphous alloy core is assembled by lacing thin ribbon layers together in seeking to attenuate the eddy current effect.Despite the better flux distribution of this structure, a different magnitude of flux density can still be observed in each ribbon due to the different magnetic paths of flux.To capture the inhomogeneous spatial distribution, the non-uniform variation of the flux density is taken into account by dividing the core into n equal segments with the thickness of dx.The flux density for the i-th segment and the total flux can be expressed as: where R t , and ϕ t denotes the total reluctance and flux of the core, while R i and ϕ i represents the corresponding parameter of the i-th magnetic segment.The reluctance of the i-th magnetic segment can be calculated as: where l i is the length of the magnetic path, ε is the material permeability, and A i is the cross-sectional area of the i-th magnetic path.For the E-E core in use, the following equation can be drawn: Combining Equations ( 2)-( 5), the following relationship in terms of the reluctance and flux density can be drawn: (6) where B i is the flux density of the i-th magnetic segment.It is explicit from Equations ( 1) and ( 6) that the actual core loss varies spatially according to this non-uniform distribution of flux density, i.e. higher core loss allocation occurs in inner locations.To this end, the proposed model scrutinizes the spatial distribution of flux density thus contributes to improving the modeling accuracy of core loss.

AC Effect and Winding Loss
The ohmic losses in the windings are calculated by: where J is current density and σ is electric conductivity.
As eddy current effects become more prominent at higher frequency, the effective AC resistance of winding increases as a consequence of the skin and proximity effects.To take eddy current effect into consideration, two well-known analytical models that are applicable to foil and round conductors, respectively, are used in this paper.The classical Dowell's equation [38] gives the AC resistance of the foil-type winding, assuming that foil conductors occupy the entire height of the core window.The AC factor of the j-th layer (Fr n ) can be calculated by: (8) where ∆ is the penetration ratio, d foil the foil thickness, δ the skin depth, and ξ 1 and ξ 2 the skin and proximity effect respectively.Dowell's equation for foil conductors can also be applied on round-type conductors by introducing the porosity factor (η), which is the ratio of the height occupied by conductors to the window height.In addition, in 1990, Ferreira proposed a new solution for round-type wire using the Kelvin-Bessel functions by considering the orthogonality between skin and the proximity effect.However, this formula overlooked the porosity factor.To take Dowell's porosity coefficient into account and to give a more accurate solution, Bartoli modified Ferreira's formula and proposed a new model as [39]: The multilayer configuration of windings are taken into account using Equations ( 8) and ( 9), so that the AC resistance for each layer can be calculated more prudently.With the conductor geometry and winding arrangement, the total AC resistance can be derived by accumulating all individual layer resistances: Dowell's equation is applicable when the transformer is subject to sinusoidal excitation.In the case of non-sinusoidal excitation, the current can be decomposed into harmonics with Fourier transform, while the loss for each harmonic is calculated by multiplying the harmonic current with the AC resistance under corresponding frequency.The total winding loss can then be determined by summing all of the harmonic losses.
It should be noted that the resistivity of the conductor copper is also temperature-dependent, suggesting that the winding loss is sensitive to the result of thermal modelling.The resistivity of copper can be calculated empirically as: where ϕ 0 is the resistivity at T 0 and τ is the temperature coefficient of copper, which is equal to 0.004041.

Thermal-Fluid Simulation
This section goes further to describe the CFD-based thermal-fluid modelling, leveraging the losses that were obtained from the electromagnetic analysis in Section 2.

Governing Equations
The heat balance for each component of the transformer is expressed as: Multiple heat transfer processes occurs for the ONAN transformer.The heat is transferred by conduction among the core, winding, and insulation layer.Meanwhile, the heat generated is dissipated from the transformer solid surfaces to the oil via convection.The oil carries the heat that is generated by convection to the tank and finally the heat is dissipated by radiation and convection to the ambient.
The heat conduction, convection, and radiation can be described, respectively, by Fourier's law, Newton's law, and Stefan-Boltzmann's law.The equations are: where Q is heat generation rate, k is the thermal conductivity, h is the convection coefficient, A is the heat-transfer area, κ is the surface emissivity, δ is the Stefan-Boltzmann's constant, h r is the equivalent radiation coefficient, and T s and T ∞ is the surface and ambient temperature, respectively.The mathematical model consists of governing equations simplified from the general expressions for the conservation of mass, momentum, and energy, as follows: where ρ is the material density, v is velocity, F is body force per unit mass, p is static pressure, µ is dynamic viscosity, Cp is specific heat capacity, and λ is heat conductivity.

Boundary Conditions
No-slip boundary condition is defined at the fluid-solid interface.The temperature and heat transfer are continuous at all interfaces.In the meantime, equal momentum and heat transfer between liquid-gas interfaces and free surface are defined.Convection heat transfer occurs when the fluid flows along the solid surface between the fluid and the solid.The convection heat transfer coefficient can be calculated with: where Nu is the Nusselt number, k f is the thermal conductivity of fluid, and L is the characteristic length.The Nusselt number for vertical and horizontal enclosures and for vertical flat plates are calculated as: 0.42Ra 1/4 Pr 0.012 (H/L) −0.3 0.069Ra 1/3 Pr 0.074 0.825 + 0.387Ra 1/6 / 1 + (0.492/Pr) 9/16 8/27 (16) where Ra, Gr, and Pr are the Rayleigh, Grashof, and Prandtl numbers of the fluid, respectively, which are defined as: where g is acceleration of gravity, λ is the fluid volumetric expansion coefficient, and ∆T is the temperature difference between surface and fluid.

Determination of Thermal Parameters
The thermal-fluid system can be solved with Equations ( 14)- (17), where the parameters are temperature-dependent.Therefore, these parameters should be included into the simulation in the form of parametric equations.The dependence of thermal parameters of the core and winding materials to the temperature is considered by incorporating the equations in Table 1.The temperature-dependent physical properties of mineral oil can be calculated by equations in Table 2.

Parameter Specification
Table 2. Thermal Parameters of Mineral Oil.

Parameter Specification
The core is assumed to be a solid compound consisting of flat steel sheets interspersed with oil.The schematic illustration of the core structure is shown in Figure 1, where θ oil and θ core represents the layer thickness of oil layer and steel sheet, respectively.The conduction process is assumed to be orthotropic with equivalent normal thermal conductivity (k eq,n ) and transversal thermal conductivity (k eq,t ), which can be calculated as: The equivalent density (ρ eq ) and specific heat (C p,eq ) of the core are calculated with: C p,eq = C p,iron ρ iron S+C p,oil ρ oil (1−S) where S is the stacking factor.

Coupled and Semi-numerical Framework
The solution begins with an analytical electromagnetic model for the determination of core loss and winding loss.Leveraging the obtained losses, the thermal-fluid simulation is performed using Ansys Fluent for steady-state analysis.The computing domain is discretized and the total number of grids is determined by performing the grid independence test to reduce the computational effort while keeping adequate accuracy [40].The temperature dependence of material properties shown in Tables 1 and 2 are taken into account with the user-defined function (UDF).The SIMPLE algorithm was employed to handle the velocity-pressure coupling in the flow field equations.The pressure discretization is performed with the body force weighted scheme.The second-order upwind scheme was applied for the space derivatives of advection terms in all transport equations.
The analytical model and the CFD model are combined in a sequential and iterative framework, as shown in Figure 2, to address the intrinsic coupling between electromagnetic and thermal behaviors.The overall convergence is achieved when the changes of losses and temperature during two adjacent iterations are no more than their respective thresholds.Based on the proposed coupled, semi-numerical framework, the steady-state losses, oil flow pattern, and temperature distribution can be obtained.
The core is assumed to be a solid compound consisting of flat steel sheets interspersed with oil.The schematic illustration of the core structure is shown in Figure 1, where θoil and θcore represents the layer thickness of oil layer and steel sheet, respectively.The conduction process is assumed to be orthotropic with equivalent normal thermal conductivity (keq,n) and transversal thermal conductivity (keq,t), which can be calculated as: eq n oil core oil oil core iron eq t oil oil iron core oil core The equivalent density (ρeq) and specific heat (Cp,eq) of the core are calculated with: where S is the stacking factor.

Coupled and Semi-numerical Framework
The solution begins with an analytical electromagnetic model for the determination of core loss and winding loss.Leveraging the obtained losses, the thermal-fluid simulation is performed using Ansys Fluent for steady-state analysis.The computing domain is discretized and the total number of grids is determined by performing the grid independence test to reduce the computational effort while keeping adequate accuracy [40].The temperature dependence of material properties shown in Table 1 and Table 2 are taken into account with the user-defined function (UDF).The SIMPLE algorithm was employed to handle the velocity-pressure coupling in the flow field equations.The pressure discretization is performed with the body force weighted scheme.The second-order upwind scheme was applied for the space derivatives of advection terms in all transport equations.
The analytical model and the CFD model are combined in a sequential and iterative framework, as shown in Figure 2, to address the intrinsic coupling between electromagnetic and thermal behaviors.The overall convergence is achieved when the changes of losses and temperature during two adjacent iterations are no more than their respective thresholds.Based on the proposed coupled, semi-numerical framework, the steady-state losses, oil flow pattern, and temperature distribution can be obtained.

Prototypes and Experimental Setup
Two ONAN single-phase transformers with different designs are investigated.The operating power is 5 kVA, while the LV (primary) and HV (high-voltage) (secondary) voltages are 500 V and 5 kV, respectively.The other major system specifications are summarized in Table 3.Both two prototypes use the P-S design and they are placed in the oil tank with the dimension of 200 mm × 200 mm × 200 mm.An experimental setup consisting of ONAN prototypes, programmable power source, load bank, monitoring units, data collection and storage units are built for testing, as shown in Figure 3. Specifically, the LV side is supplied with 1 kHz voltage by using Pacific AC source with programmable controller and output impedance, while the HV side is regulated with programmable load bank to achieve the expected power.The current and voltage are continuously monitored with DSOS054A High-Definition Oscilloscope during the tests and load experiments.The real-time power is measured by aWT3000E Precision Power Analyzer.

Loss Measurement
The core loss is measured by supplying the LV side with the expected voltage and opening the HV side.The open circuit setup allows for full core excitation but extremely low currents in the windings.Therefore, winding losses are negligible and the core loss is approximately equal to the observed total active power, which can be obtained by calculating the mean of the product of current and voltage.The winding loss can be measured by feeding the HV side with sinusoidal voltage while leaving the LV side short circuited.The voltage is stepped up gradually until the current reaches the expected value (1 A in HV side).The flux density in core during the short-circuit test is negligible as the HV side is excited with minimal voltage.In this case, the core loss is ignorable and the observed total loss is the HV winding loss.

Temperature Measurement
Direct temperature measurement, although quite costly, shows a high accuracy thus has been widely used for model validation.A certain number of temperature measurement points have been selected for each prototype, as shown in Figure 4, to collect the reference temperature and compare with the predicted values.J-type thermocouples are inserted to the illustrated positions during the manufacturing process.The accuracy of temperature sensing is ±2.

Validation of Losses Determination
The core losses and winding losses determined by the proposed semi-numerical model and the full numerical model in comparison with the real measurements are shown in Tables 4-5, respectively.It is shown that the modelling results are in reasonable agreement with the measured values, suggesting that the proposed method can simulate the electro-magnetic behaviours of HF transformers well.Meanwhile, the semi-numerical model provides similar accuracy when compared with the full numerical model.Hence, it is validated that the reduction of computing cost does not compromise the modelling accuracy.
The proposed method also scrutinizes the AC effects in each winding layer.Due to the difficulty of measuring layer losses, the proposed semi-numerical model is compared with the FEM-CFD method for validation.The FEM-CFD approach is a full numerical method that uses the same structure as in Figure 2 while the analytical part is replaced by 3D steady-state FEM with Ansys Maxwell.The comparison results of layer losses are shown in Figure 5.It is shown that the

Validation of Losses Determination
The core losses and winding losses determined by the proposed semi-numerical model and the full numerical model in comparison with the real measurements are shown in Tables 4-5, respectively.It is shown that the modelling results are in reasonable agreement with the measured values, suggesting that the proposed method can simulate the electro-magnetic behaviours of HF transformers well.Meanwhile, the semi-numerical model provides similar accuracy when compared with the full numerical model.Hence, it is validated that the reduction of computing cost does not compromise the modelling accuracy.
The proposed method also scrutinizes the AC effects in each winding layer.Due to the difficulty of measuring layer losses, the proposed semi-numerical model is compared with the FEM-CFD method for validation.The FEM-CFD approach is a full numerical method that uses the same structure as in Figure 2 while the analytical part is replaced by 3D steady-state FEM with Ansys Maxwell.The comparison results of layer losses are shown in Figure 5.It is shown that the

Validation of Losses Determination
The core losses and winding losses determined by the proposed semi-numerical model and the full numerical model in comparison with the real measurements are shown in Tables 4 and 5, respectively.It is shown that the modelling results are in reasonable agreement with the measured values, suggesting that the proposed method can simulate the electro-magnetic behaviours of HF transformers well.Meanwhile, the semi-numerical model provides similar accuracy when compared with the full numerical model.Hence, it is validated that the reduction of computing cost does not compromise the modelling accuracy.
The proposed method also scrutinizes the AC effects in each winding layer.Due to the difficulty of measuring layer losses, the proposed semi-numerical model is compared with the FEM-CFD method for validation.The FEM-CFD approach is a full numerical method that uses the same structure as in Figure 2 while the analytical part is replaced by 3D steady-state FEM with Ansys Maxwell.The comparison results of layer losses are shown in Figure 5.It is shown that the layer losses calculated by the semi-numerical method are benchmarked with FEM results.The mean relative deviations of two methods are 8.9% and 4.9% for prototypes 1 and 2, respectively.A scrutiny of the deviation suggests that the proposed method underestimates the losses slightly as compared to the full numerical method.This is due to the existence of edge effect, which is not considered by the analytical model.layer losses calculated by the semi-numerical method are benchmarked with FEM results.The mean relative deviations of two methods are 8.9% and 4.9% for prototypes 1 and 2, respectively.A scrutiny of the deviation suggests that the proposed method underestimates the losses slightly as compared to the full numerical method.This is due to the existence of edge effect, which is not considered by the analytical model.The deviation of two methods can be better explained by analysing the current density and magnetic field distribution on conductors and core window, as shown in Figure 6.It is shown that the proximity effect is quite obvious at the middle part of winding in the vertical direction, and tends to be severe at the boundary between the LV and HV windings.The proposed semi-numerical model shows to well describe the proximity effect.In contrast, the current density distribution is highly non-uniform at the upper and bottom corner of LV winding.This is attributed to the edge effect that is caused by the orthogonal component of magnetic field around the winding edge, as shown in the upper zoom-in figure of Figure 6a,b.While the proposed method well considers the skin and proximity effect, it overlooks the edge effect, so that its underestimation of AC losses is within expectation.
The edge effect observed in prototype 2 is not as strong as in the case of prototype 1.This is because prototype 2 has a larger filling factor (the ratio of winding height to the core window height), which prohibits the tangential component of magnetic field at the edge regions.Therefore, a higher filling factor is always favourable to alleviate the extra losses that are caused by the edge effect.The milder edge effect also coincides with the aforementioned result that the deviation of The deviation of two methods can be better explained by analysing the current density and magnetic field distribution on conductors and core window, as shown in Figure 6.It is shown that the proximity effect is quite obvious at the middle part of winding in the vertical direction, and tends to be severe at the boundary between the LV and HV windings.The proposed semi-numerical model shows to well describe the proximity effect.In contrast, the current density distribution is highly non-uniform at the upper and bottom corner of LV winding.This is attributed to the edge effect that is caused by the orthogonal component of magnetic field around the winding edge, as shown in the upper zoom-in figure of Figure 6a,b.While the proposed method well considers the skin and proximity effect, it overlooks the edge effect, so that its underestimation of AC losses is within expectation.
The edge effect observed in prototype 2 is not as strong as in the case of prototype 1.This is because prototype 2 has a larger filling factor (the ratio of winding height to the core window height), which prohibits the tangential component of magnetic field at the edge regions.Therefore, a higher filling factor is always favourable to alleviate the extra losses that are caused by the edge effect.
The milder edge effect also coincides with the aforementioned result that the deviation of two methods is smaller for the case of prototype 2. It should be mentioned that the skin effect is not obvious for both the two prototypes, as the skin depth under the applied frequency is at the same level with the characteristic length of conductors.

Validation of Temperature Prediction
In seeking to validate the temperature prediction, the proposed method is validated by comparing with both the full numerical model and experimental results.Three simulations have been performed for each prototype to validate the creditability of the simulation results.The simulation result of temperature distribution is shown in Figure 7.It is explicit that a temperature gradient exists along the vertical direction of transformer active components and oil.This is because the hot oil keeps rising to the upper section of tank that is driven by the buoyancy force.The central core and inner winding layers show higher temperature due to significant losses and the poor heat dissipation condition.The predicted temperatures at measurement points are compared with experimental data to evaluate the prediction accuracy.The temperatures obtained from the semi-numerical model and full-numerical model in comparison with measurements are shown in Figure 8.The relative deviations between simulation and measurement calculated by (T simu − T meas )/T meas are summarized in Table 6.It is shown that the predicted temperature by both the two models is benchmarked with the measurement data with reasonable accuracy.The prediction errors for the central core and windings are well confined to approximately 10%.In contrast, the proposed semi-numerical model overestimates the temperature of the external core, as seen from the relatively large prediction errors exceeding 10%.One possible explanation is that the model for thermal conductivity calculation in Figure 1 may be deviated from real prototypes, due to the uncertainties in the manufacturing process.The proposed method is experimental verified with sufficient accuracy in temperature prediction, especially when considering the unavoidable mismatch between designed geometry and real prototypes.Moreover, the proposed semi-numerical model shows to give similar accuracy when compared with the computationally expensive full-numerical model, suggesting that the proposed technique lowers the computing cost without compromising the modelling accuracy.

Discussion
It should be noted the proposed semi-numerical model is a general one that is able to provide accurate prediction of losses and the temperature response of transformer.This section goes further to discuss the extendibility of the proposed model by analysing each sub-model, i.e. core loss model, winding loss model, and thermal-fluid model.The refined core loss model can be easily extended to any magnetic material, structure, and frequency by scrutinizing the special distribution of the flux density.The layer-based winding loss model takes the skin effect, proximity effect, and radial geometric inhomogeneity into account to better describe the resistive feature.From theoretical point of view, the proposed model is applicable at lower frequency, such as 50/60 Hz, since the AC effects are not substantial in this case and they are generally easier for modeling.In contrast, albeit aiming

Discussion
It should be noted the proposed semi-numerical model is a general one that is able to provide accurate prediction of losses and the temperature response of transformer.This section goes further to discuss the extendibility of the proposed model by analysing each sub-model, i.e. core loss model, winding loss model, and thermal-fluid model.The refined core loss model can be easily extended to any magnetic material, structure, and frequency by scrutinizing the special distribution of the flux density.The layer-based winding loss model takes the skin effect, proximity effect, and radial geometric inhomogeneity into account to better describe the resistive feature.From theoretical point of view, the proposed model is applicable at lower frequency, such as 50/60 Hz, since the AC effects are not substantial in this case and they are generally easier for modeling.In contrast, albeit aiming to improve the modeling of AC effects, the proposed model is subject to an accuracy declination at higher frequency, since secondary AC effects, such as the edge effect have been overlooked.This has already been shown in Figures 5 and 6, where explicitly an underestimate of winding losses has been observed.With respect to the thermal-hydraulic model, the numerical model leveraging CFD is performed to obtain the detailed flow characteristic, thermal path, and temperature distribution.This sub-model stands for the state of art for modeling of the thermal behavior of transformer and enjoys a good extendibility in terms of the selection of material, geometry, operating conditions, and cooling medium, provided that the losses can be determined with reasonable credibility.

Conclusions
This paper proposes a coupled, semi-numerical model for the thermal analysis of MF ONAN transformers.An analytical model has been exploited based on spatial distribution of flux density and AC factor to calculate the core and winding losses, while the CFD simulation has been used for thermal-fluid analysis.The electromagnetic and thermal-fluid behaviours are solved iteratively in a close-loop manner by leveraging the analytical and CFD model, respectively.Lab-scale experiments have been carried out on two prototypes with different designs.Results show that the AC effects, losses, and temperature rise can be well predicted by the proposed method.The predicted temperature is benchmarked with experimental data with a reasonable accuracy.The proposed semi-numerical model is technically feasible and experimentally validated, so that it is promising to be used for transformer design optimization and operating performance prediction.It is also preferable when compared to the coupled full numerical approach due to the better trade-off between modelling accuracy and computing cost.

Figure 1 .
Figure 1.Structural diagram for core with flat steel sheets and oil layer.

Figure 2 .
Figure 2. Overall framework of the coupled and semi-numerical method.
2 • C within an operating range from 0 • C to 760 • C. Both prototypes are operated with 5 kVA to obtain the transient temperature response.A laptop that was equipped with BV0006B Data Acquisition Control & Automation software is connected to the 34972A Data Acquisition Switch Unit to acquire the temperature.Load experiments are terminated once the temperature variations are less than 1 • C per hour, which serves as the steady state criterion in this study.Energies 2019, 12 FOR PEER REVIEW 10 software is connected to the 34972A Data Acquisition Switch Unit to acquire the temperature.Load experiments are terminated once the temperature variations are less than 1℃ per hour, which serves as the steady state criterion in this study.

Figure 4 .
Figure 4. Thermocouple locations (Note: the figure does not represent the real geometry of transformer design; it is a diagram to identify the relative locations of temperature measurement points).

Figure 4 .
Figure 4. Thermocouple locations (Note: the figure does not represent the real geometry of transformer design; it is a diagram to identify the relative locations of temperature measurement points).

Figure 4 .
Figure 4. Thermocouple locations (Note: the figure does not represent the real geometry of transformer design; it is a diagram to identify the relative locations of temperature measurement points).

Energies 2019, 12 FOR PEER REVIEW 13 Figure 7 .
Figure 7. Temperature distribution in transformer active components and oil for the two investigated prototypes.

Figure 7 .
Figure 7. Temperature distribution in transformer active components and oil for the two investigated prototypes.

Figure 7 .
Figure 7. Temperature distribution in transformer active components and oil for the two investigated prototypes.

Figure 8 .
Figure 8.Comparison between simulated and measured temperature at different locations: (a) Prototype 1 and (b) Prototype 2.

Figure 8 .
Figure 8.Comparison between simulated and measured temperature at different locations: (a) Prototype 1 and (b) Prototype 2.

Table 1 .
Thermal Parameters of Core and Winding Materials.

Table 3 .
Specifications of transformer prototypes.

Table 4 .
Comparison of modeled core losses and measurements.

Table 5 .
Comparison of modeled winding losses and measurements.

Table 4 .
Comparison of modeled core losses and measurements.

Table 5 .
Comparison of modeled winding losses and measurements.

Table 6 .
Relative deviations between predicted and measured temperatures.

Table 6 .
Relative deviations between predicted and measured temperatures.