Efficient Numerical Modeling of Oil-Immersed Transformers: Simplified Approaches to Conjugate Heat Transfer Simulation

Abstract

The development of digital twins for power transformers has become increasingly important to predict possible operating modes and reduce the likelihood of faults. The accuracy of these predictions relies heavily on the numerical models used, which must be both simple and computationally efficient. This work focuses on creating a simplified numerical model for a template oil-immersed power transformer (100 MVA, 230/69 KV). The study investigates how the number of elements and the strategies used to set up the mesh in the domain of interest influence the results, aiming to identify the key parameters that affect the outcomes. Furthermore, a significant effect of resolving thermal boundary layers on the accurate identification of hot spots is demonstrated. Two approaches to resolving thermal boundary layers are explored in this work. This study presents a comprehensive analysis of three numerical models for conjugate heat transfer simulations, each with distinct features and computational domain compositions. The results show that the addition of extra calculation domains leads to the emergence of new vortex structures, affecting the velocity profile at the channel inlet and altering the location of hot spots. This study provides valuable insights into the configuration and composition of calculated domains in numerical models of oil-immersed power transformers, essential for the accurate prediction of hot spot temperatures and ensuring reliable operation.

1. Introduction

Oil-immersed power transformers are commonly used in distribution networks, industry and domestic applications. Their technical state must be predicted to prevent any faults. A fault in a power transformer might lead to an accident in a plant, network or building and, finally, to high economic costs and dangerous incidents. Controlling the highest temperature in the winding is an important task for maintaining the health index of the transformer [1]. Overheating of the power transformer winding leads to the severest thermal aging of the paper insulation and oil due to cavitation or bubbles appearing in the oil [1,2,3].

One of the possible methods for preventing the failure of a transformer and predicting its current resource is to use digital twins technology based on numerical simulations of physical phenomena in transformers [4,5]. In [6,7], a digital twin model for power transformers was developed to predict hot spot temperature. Modern digital twin modeling methods for transformer temperature prediction encompass both model-driven and data-driven approaches [6]. In model-driven approaches, real-scenario and numerical multi-physics problems are considered. The reliability of predictions made by the digital twin system hinges critically on two key factors: the precision of the underlying numerical model and the swiftness with which it can deliver results to inform timely decision-making.

The physical modeling of power transformers to predict hot spot temperatures necessitates a comprehensive framework that accounts for hydrodynamic, heat transfer and electromagnetic processes within both solid and liquid components. This challenge is particularly pronounced in objects featuring intricate geometries, such as oil-filled transformers. The narrow gaps between windings and magnetic cores give rise to complex hydrodynamic channels, which demand a highly detailed mesh to accurately capture the system’s behavior. Consequently, a compromise must be struck between model complexity and computational time.

Conjugate heat transfer simulation can be achieved through methods based on circuit theory [8,9,10,11] or numerical methods [4]. For instance, a hydraulic network model has been suggested for predicting reverse flow in oil-directed cooling modes for disk-type transformer windings [12]. The application of this method allows for a significant reduction in computation time and exhibits good convergence with computational fluid dynamics (CFD) solutions [12]. However, such approaches become unsuitable for complex transformer geometries or when intricate convection patterns are formed within the oil flows. Furthermore, the resulting model requires reassembly of equivalent circuits for each design iteration, as demonstrated in [4].

Numerical analysis provides an alternative to circuit-based models, enabling more accurate results and a deeper understanding of the study object. However, numerical methods such as finite element and finite volume methods require significantly greater computational resources than those based on circuit theory. For instance, hot spot temperatures in a power transformer can be calculated using hydraulic network models in mere seconds. In contrast, simulations of the natural convection flows in a sector of a power transformer using COMSOL Multiphysics software with approximately 1 million finite elements can take around 90–120 min of computational time on a workstation equipped with two 2.66 GHz processors and 48 GB of RAM [13]. The high computational cost associated with these models becomes particularly apparent when conducting studies that require the evaluation of alternative coolants. The transition from 3D models to 2D numerical models allows the requirement of hundreds of GB of RAM to be reduced to dozens, as was achieved in [14]. A significant reduction in the computational cost for the simulation of the complex convective flows of oil in power transformers was achieved in [15]. The simulation considered only a simplified part of the power transformer design and was conducted in OpenFOAM software. The authors [1,16,17] implemented in COMSOL Multiphysics software a similar approach to define hot spots in high-voltage winding, considering in the numerical model only certain parts of the power transformers. However, the accuracy of the results depends on the settings of the chosen models and physical assumptions, and these questions are still discussed.

It is worth noting that numerical models may be augmented or entirely replaced by reduced-order models, such as those based on proper orthogonal decomposition (POD) [18,19] and dynamic mode decomposition (DMD) [20]. Additionally, semi-analytical approaches [21,22] or physics-informed neural network prediction models [23] can be employed to significantly enhance computational efficiency. For example, the use of physics-informed neural networks in [24] enabled the rapid evaluation of top oil and heat loss distribution in power transformers. This approach facilitates the description of numerical solutions for physical laws formulated on nonlinear partial differential equations.

The computational cost of numerical models can be significantly reduced by using optimized mesh settings or by excluding certain regions from the simulation domain. However, some key mesh parameters are often listed without proper justification in the literature [13,14,15], making it unclear how simplifications impact results. Even though previous studies [14,25,26,27] have investigated reducing the calculation domain, a systematic analysis is still needed. Power transformers are complex systems influenced by various factors that affect fluid dynamics and temperature distribution. Evaluating how different levels of model detail impact results—from simple models that capture only basic physical flow behavior to more detailed ones—can be highly beneficial. This approach helps to illustrate how results change when certain physical phenomena are ignored. The primary focus of this study is to investigate how simplifying the conjugate heat transfer model affects the results from both qualitative and quantitative perspectives.

This paper presents original findings on mesh settings for numerical models of conjugate heat transfer in an oil-immersed transformer. Based on these results, guidelines are provided for configuring the numerical model mesh to minimize computational resources while maintaining accuracy. The specifics of modeling natural convection flows of oil in parts cooled by radiators and the bottom and top domains of a tank are considered by three suggested numerical models. A detailed consideration of these aspects is essential for gaining a deeper understanding of the power transformer’s thermal behavior. The importance of considering the characteristics of the cooling system was highlighted in [28,29]. The suggested three numerical models are described in detail from a mathematical point of view and developed by using a number of new combinations of boundary conditions to reduce simulation domains. This paper is a continuation of the effort to generate new knowledge in the field of numerical simulation of conjugate heat transfer in oil-immersed transformers, which is expected to help future engineers create their own numerical models based on the requirements for the accuracy of results and existing computational resources.

The paper is organized into four main sections: the Introduction (Section 1), which introduces the problem statement and research objectives; the development of simplified conjugate heat transfer models for oil-immersed power transformers (Section 2), which describes the novel numerical approaches proposed in this study for modeling conjugate heat transfer in oil-immersed power transformers; the Results (Section 3), where computational simulations are presented and their features specified; and the Discussion (Section 4) and Conclusion (Section 5), which interpret the findings and provide recommendations for their applications.

2. Development of Simplified Conjugate Heat Transfer Models for Oil-Immersed Power Transformer

2.1. Model Geometry Description

Evaluation of hot spot temperature using numerical models can be achieved by considering only those parts of the transformer where heat generation and dissipation are dominant. According to the research [6], approximately 90% of the total heat energy is dissipated in the windings and core, with heat removal taking place in oil domains. The computational requirements dramatically increase when new domains are included in the numerical model. Therefore, we consider a simplified design of a three-phase oil-immersed power transformer (100 MVA, 230/69 kV), shown in Figure 1.

The transformer consists of a magnetic circuit, depicted as the “Core”, three HV and LV windings, shown as “Coils”, and a tank with nozzles filled with transformer oil proposed for cooling the coils and core. The nozzles support the removal of overheated oil from the tank cavity to an auxiliary cooling system, for example, with radiators or any other installations. This auxiliary cooling system is not drawn in Figure 1. The main geometric parameters of the power transformer can be found in [5,30] and are not listed in this work because that part of the design will not be treated in the numerical model. The design of the power transformer is considered only as a template for preparing a simplified model geometry considering only the main parts of the transformer.

The numerical models assume a symmetrically loaded operation mode for the power transformer. In this case, heat sources in each electrical phase are uniformly distributed, allowing us to consider only one electrical phase for analysis. This simplification enables the examination of three oil flow paths, which can be extruded to study conjugate heat transfer in the power transformer. The first channel is the gap between the core and low-voltage winding. The second channel is the gap between the low- and high-voltage windings. The third channel is the space between the high-voltage winding and the tank wall. The core of the transformer can be modeled as a cylindrical geometry part, resulting in a corresponding cylindrical shape for this channel.

Previous studies, such as those discussed in [31], employed a 3D symmetric formulation to investigate power transformer design. In contrast, this study adopts a 2D symmetric formulation for the numerical model. The transition from a 3D to a 2D model can be justified by the absence of cardboard cylinders separating the high-voltage windings, which would otherwise require treatment as solid obstacles in hydrodynamic simulations. The presence of non-symmetric flows in these channels will primarily result in natural convection or turbulence flows. At this stage of the study, non-symmetric flows are neglected.

In Figure 1, the simulation area, highlighted by the green color, has its edges between the liquid and solid parts marked by the red color. A 2D cut plane covers the core, coils and oil-immersed domain of the tank with nozzles, which are depicted in Figure 2. The geometry parameters listed in

Table 1. The geometry parameters.

Symbol	Value	Description
$R_c$	370 mm	The radius of the core
$R_{w1}$	417 mm	The internal radius of the LV coils
$R_{w2}$	488 mm	The external radius of the LV coils
$R_{w3}$	546 mm	The internal radius of the HV coils
$R_{w4}$	642 mm	The external radius of the HV coils
$D_d$	58 mm	The width of right channel
$D_n$	58 mm	The width of nozzles and pipe branch
$H_w$	1660 mm	The height of windings and core
$H_a$	800 mm	The height of extra domain

correspond to the configuration shown in Figure 2. Notably, the yoke of the core is neglected, and the height of both the core and windings is equal. It should be noted that considering only one electrical phase in a three-phase transformer represents a significant simplification of the system. In reality, each phase operates under different conditions, particularly with respect to heat transfer. The phases are not isolated from each other. For example, the middle phase is influenced by interactions with the other phases, which can affect its thermal behavior. However, this simplification is deliberately adopted to explore how the simplest model, with different levels of detail, can impact the overall results. Consequently, certain physical phenomena, such as the mutual thermal effects between phases and detailed heat transfer mechanisms, are neglected in order to prioritize the computational feasibility and focus on the broader trends of the system’s behavior.

2.2. Different Types of Numerical Models

We examine three types of numerical models, each with distinct features. The primary difference between these models lies in their computational domain composition. The first model includes domains for core, left, middle and right channels, as well as LV and HV windings. It considers only conjugate heat transfer due to the direct interaction between the transformer’s liquid and solid components. This model focuses on the primary heat transfer mechanism within the transformer.

In contrast, the second model is designed to account for vortices and complex flow structures outside the channels. To achieve this, it incorporates the entire domain of the first model plus an additional domain, depicted in Figure 2. The comparison of results calculated by the first and second models enables us to understand the influence of vortex generation on transformer design features.

The third model takes into account flows appearing in pipe branches. For simplification, the closed pipe branch is designed with a simpler shape than its real-life counterpart, solely for evaluating the impact of geometry features on oil flow structures. A summary of these numerical models can be found in

Table 2. The summary of numerical models.

Marker of the Model	Including Domains
Model #1	Core, HV and LV windings, left, right and middle channels
Model #2	Model #1 domains and extra domain
Model #3	Model #2 domains and pipe branch domain

.

2.3. Material Properties

The domains comprising LV and HV windings consist of copper material. The transformer steel serves as the core, while the remaining domains are filled with transformer oil. The physical properties utilized in this study are listed in

Table 3. Physical property parameters.

Materials	Winding	Core	Oil
Mass density, kg/m3	8940	7550	879
Heat capacity, J/(kg·K)	385	446	1711
Thermal conductivity, W/(m·K)	407	72	0.11
Dynamic viscosity, Pa·s	-	-	0.02

. To simplify the thermal and hydrodynamic boundary mesh calculations and analysis of mesh size influence on the results, temperature dependence is ignored. If considered, the physical properties’ dependency on temperature would introduce a complicating factor. As a result, the mesh layer would become non-uniform across different simulation domains, leading to a complex understanding of mesh sensitivity.

It is important to note that all the models described do not account for the insulation wound around the transformer windings. This assumption is commonly made to simplify calculations [32,33]. Insulation influences the resulting physical fields, as discussed in detail in [34]. The authors conclude that, while including insulation reduces the temperature error by 2.5%, it has little impact on the location of the hot spot. For this reason, and because this study focuses on simplifying the primary factors influencing conjugate heat transfer, insulation is not considered in the present analysis.

2.4. Governing Equations

The numerical model is built in COMSOL Multiphysics using the finite element method (FEM) to solve both temperature and hydrodynamic fields. Since the flow under consideration is convective dominated, it is important to implement a stabilization scheme during the calculation using the FEM. In this study, consistent stabilization techniques are used, specifically streamline and crosswind diffusion. The motion of oil domains is treated by solving the continuity mass Equation (1) and the moment Equation (2), which apply to incompressible flow.

$$\begin{array}{c}\hfill \rho \nabla ·\mathbf{u}=0.\end{array}$$

(1)

$$\begin{array}{c}\hfill \hspace{1em}\rho \left[{\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+\left(\mathbf{u}·\nabla \right)\mathbf{u}\right]=\nabla \left[-p\mathbf{I}+\mu \left(\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^T\right)\right]+{\mathbf{f}}_b.\end{array}$$

(2)

Here, $\mathbf{u}$ is the velocity, p is the pressure, ${\mathbf{f}}_b$ is the volumetric buoyancy force, $\rho$ is the mass density, $\mu$ is the dynamic viscosity and $\mathbf{I}$ is the unit matrix. The buoyancy force ${\mathbf{f}}_b$ is simulated using the Oberbeck–Boussinesq approximation (3):

$${\mathbf{f}}_b=\rho \mathbf{g}={\rho }_0\mathbf{g}\left(1-{\beta }_T\right)\left(T-T_0\right),$$

(3)

where $\mathbf{g}$ is the gravity force, ${\rho }_0$ is the mass density at the reference temperature $T_0$, ${\beta }_T$ is the thermal expansion coefficient of the oil (1/K) and T is the temperature field at the current time step. The thermal expansion coefficient can be evaluated by

$${\beta }_T=-{\displaystyle \frac{\partial \rho \left(T\right)}{\partial T}}{\displaystyle \frac{1}{\rho \left(T\right).}}$$

(4)

and set to constant values $7·{10}^{-4}$$K^{-1}$ for most transformer oil in the temperature range of 0–120 °C. According to [35], the Boussinesq approximation is valid when temperature variations in the flow are small. If density variations are at most a few percent, it follows that ${\rho }^{-1}(D\rho /Dt)$ can also be no larger than a few percent of the velocity gradients in $\nabla ·\mathbf{u}$. Consequently, density can be considered constant in both the continuity and momentum equations, except within the gravity term. In this study, the relationship for oil density is expressed as follows (Figure 3):

$$\rho \left(T\right)=1055.05-0.58·T-6.4·{10}^{-5}·T^2$$

In the considered temperature range, the density change does not exceed 7%. Such an approach is also implemented in a number of works [33,36,37,38].

The energy Equation (5), as expressed for the temperature field, is calculated to evaluate both the hot spot temperature and the buoyancy force (3) within the moment Equation (2). In particular, it should be noted that the term $\mathbf{u}·\nabla T$ is not considered for the solid parts of the numerical model of the power transformer.

$$\rho c_p\left[{\displaystyle \frac{\partial T}{\partial t}}+\mathbf{u}·\nabla T\right]-\lambda {\nabla }^{2}T=q_{heat}$$

(5)

In (5), $c_p$ represents the specific heat capacity, $\lambda$ denotes the thermal conductivity and $q_{heat}$ signifies a volumetric heat source produced by Joule losses in the windings and core. Reference values for the Joule losses are preliminary calculated using a 3D formulation of the numerical model for the harmonic magnetic field expressed from the magnetic vector potential. The reference values, denoted as $q_{ref}$, are kW/m³ for the LV and HV windings and 1.1 kW/m³ for the core domain. Simulation of the operating conditions of the power transformer in the numerical model is implemented by varying the heat load coefficient $k_{heat}$ according to the expression

$$q_{heat}=k_{heat}{q}_{ref}.$$

(6)

The implicit discretization scheme is employed, specifically the backward difference formula (BDF) of order 1–2. Time stepping is restricted by the Courant–Friedrichs–Lewy (CFL) condition to ensure numerical stability and accuracy.

2.5. Boundary Conditions

General boundary conditions applicable to all models (described in Section 2.2) are considered further. First, non-slip boundary conditions ($\mathbf{u}=0$) are applied to all edges separating solid and liquid domains, which are colored black. Secondly, a symmetric boundary condition is imposed on the red dashed edge labeled “BC #1”, and, finally, the tank walls are assumed to be perfectly thermally insulated, as described by Equation (7).

$$-\mathbf{n}·\nabla T=0$$

(7)

In models #1 and #2, the simulation domain of liquid is restricted by the edges colored turquoise (BC #2) and blue (BC #3) in Figure 2, respectively. On these edges, the open boundary condition is set as for pressure–velocity fields and temperature fields. This boundary condition allows natural consideration of flows caused by convection, particularly appearing as backflows in a number of modes. The hydrodynamic Equations (1) and (2) are then treated together with Equation (8) on the “BC #2” edges (depicted in Figure 2) for model #1, and the “BC # 3” edges (in Figure 2) for model #2.

$$\left(-p\mathbf{I}-\mu \left[\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^T\right]\right)·\mathbf{n}=-{\rho }_0\mathbf{g}·\left(\mathbf{r}-{\mathbf{r}}_{\mathbf{ref}}\right)·\mathbf{n},$$

(8)

where $\mathbf{r}$ and ${\mathbf{r}}_{\mathbf{ref}}$ are the positional vector and reference position vector, respectively.

The open boundary condition of the temperature field depends on the sign of velocity at the edges labeled “BC #3” and “BC #2”. This condition is described by (9) when the oil flow is outgoing from the simulation domain, i.e., $\mathbf{u}·\mathbf{n}<0$. Conversely, for incoming flows where $\mathbf{u}·\mathbf{n}\ge 0$, the boundary condition is set according to (7).

$$-\mathbf{n}·\lambda \nabla T=\rho {\int }_{T_0}^T{c}_{p}dT\mathbf{u}·\mathbf{n}$$

(9)

The upper and bottom horizontal edges of the core and LV and HV windings are described by a thermal insulation boundary condition (7) for model #1. The black dashed edges of the pipe branch used in model #3 are simulated as a non-slip boundary condition for the hydrodynamic equations and prescribing a heat flux (10) for the temperature field equation.

$$-\mathbf{n}·\mathbf{q}=h\left(T_0-T\right)$$

(10)

Here, $h={10}^4$ W/(m·K) is expressed as the heat transfer coefficient. An extreme value of the heat transfer coefficient is specifically chosen to simulate an ideal cooling system and intensify temperature gradients for the purposes of this study. This approach helps us to better understand the relative importance of different thermal mechanisms in the model. The boundary condition allows the oil temperature in pipe branches to be reduced instead of considering the real design of the cooling system, including the complex geometry of pipes and radiators.

2.6. Meshing and Computational Information

The simulation domain, shown in Figure 4, uses a consistent mesh configuration across all models. The regions above and below the channels to the nozzle are divided into 20 elements ($N_z^{ext}=20$), with each subsequent element increasing by a factor of 5. The remaining exterior domain is also divided into 20 elements along the z-axis ($N_z^{ext2}=20$). The upper and lower nozzles have 20 elements in radial direction ($N_r^{nozzle}=20$), increasing by a factor of 5 towards the pipe branch, which is resolved with 5 uniform radial elements ($N_r^{pipe}=5$). The elements in the azimuthal direction for the channels and windings range from 20 to 400 ($N_z$), while those in the radial one in the windings ($N_{solid}$) range from 5 to 40. The core is divided into 30 elements ($N_{core}=30$), with each increasing by a factor of 5 towards the axis of symmetry.

Additionally, a specific thermal boundary layer (TBL) forms at the area of the largest temperature gradient. This TBL should be resolved by a sufficient number of mesh elements to accurately capture its effect. The thickness of the thermal boundary layer, denoted as ${\Delta }_{TBL}$, can be expressed using the values of the hydrodynamic boundary layer, denoted as ${\Delta }_{HBL}$, and the Prandtl number, defined as $Pr={\displaystyle \frac{\mu c_p}{\lambda }}$:

$$\begin{array}{c}\hfill {\Delta }_{TBL}={\Delta }_{HBL}·P{r}^{-1/3}\end{array}$$

(11)

$$\begin{array}{c}\hfill {\Delta }_{HBL}=4.91·L·R{e}^{-1/2}\end{array}$$

(12)

where the Reynolds number $Re=\rho UL/\mu$, and U and L are the characteristic velocity and length, respectively. The characteristic length may be set to half of the channel width, while the characteristic velocity is typically taken as the average velocity of oil in the channel. The thermophysical properties of transformer oil used in industry result in thermal boundary layers that comprise 2.5–38% of the hydrodynamic layer, depending on the type of oil and its temperature. Resolving the thermal boundary layer also resolves the hydrodynamic layer.

Building meshes in channels involves two approaches to resolve thermal boundary layers. The first approach does not explicitly resolve the thermal boundary layers in the channel. Instead, a non-uniform distributed number of elements $N_r$ is used to describe the mesh, as shown in Figure 4b. This distribution follows a geometric progression with a growth coefficient of $k_g=20$ where the element size grows from the walls to the center of the channel. The resulting mesh is referred to as a “non-uniform mesh”. Results obtained using this mesh will be presented in Section 3.1.

In contrast, the second approach explicitly resolves the thermal boundary layers in the channels. To achieve this, the domain of the channels is divided into a flow core area and a boundary layer region. The number of elements in the boundary layer $N_{TBL}$ and flow core $N_r$ are varied within the range of 5–40 and distributed as shown in Figure 4c. This mesh is referred to as a “non-uniform mesh with resolved TBL”.

The reference characteristics of the available computational resources are as follows: 2 × Intel(R) Xeon(R) E5-2630 v4 processors, each running at 2.20 GHz, providing a total of twenty cores across two sockets. Additionally, the system has 240 GB of RAM available.

3. Results

3.1. Study of Settings of Uniform Structured Mesh Without Explicit Resolution of Thermal Boundary Layer for Model #1

The simulation results for the non-uniform mesh are presented in this section, as described in Section 2.6. The sensitivity of these results to radial ($N_r$) and azimuthal ($N_z$) discretization of the mesh can be evaluated by examining the average velocities at the channel outlets and the maximum temperatures in the high- and low-voltage coils. Figure 5 shows the calculated average velocities at the channel outlets for a heat load coefficient $k_{\mathrm{heat}}=1$. The color of the lines represents the number of radial mesh elements, while the line style indicates the number of azimuthal mesh elements. This representation will be used throughout the paper. A coarse resolution in the azimuthal direction leads to the appearance of artificial wavy behavior in the average velocity curves for all channels considered. This wavy behavior signifies a fluctuating velocity curve over time. Increasing the number of mesh elements in the z-direction reduces this wavy behavior, and it is most pronounced for cases where $N_z=100$. In contrast, radial discretization only affects the values of velocity and does not change their shape. All presented results can be evaluated by comparing them to the finest mesh settings ($N_r=30$, $N_z=300$). The relative deviation of velocities calculated using the finest mesh and mesh cases with increasing radial elements is reduced by about 1.5–2.6 times. Conversely, an increase in the azimuthal elements leads to a decrease in relative velocity values from the smallest grid case by an order of magnitude.

The accurate prediction of hot spot temperatures plays a crucial role in the development and implementation of digital twin technologies, which rely on numerical simulations of physical phenomena within transformers. Figure 6 illustrates the maximum temperatures of low- and high-voltage windings for various mesh ratios, as a function of time, with a heat coefficient $k_{\mathrm{heat}}=1$. Notably, the radial number of elements in the mesh does not significantly impact the dynamic behavior of the maximum temperatures within both the high- and low-voltage windings. Conversely, it is observed that the influence of the mesh on the calculated temperature values weakens when compared to the maximum velocity values. Specifically, as the size of the elements in the azimuthal direction decreases, the relative deviation of the maximum temperature drops from 0.01 to 0.003, corresponding to absolute deviations of 0.36 °C and 0.074 °C, respectively.

The considered thermal mode of the power transformer leads to non-significant heating-up windings. It is important to consider mesh sensitivity under stronger modes, for example, when the heat load coefficient $k_{heat}$ equals 10. The curves over time of the maximum velocity at the inlets of the left, middle and right channels are shown in Figure 7. Notably, a significant deviation can be observed between results calculated for $N_z=100$ and those obtained from other mesh settings. Specifically, the deviation amounts to approximately $7.5\mathrm{mm}/\mathrm{s}$ for $N_z=100$, regardless of the radial number of elements used in the finest mesh case. Conversely, calculations conducted with other mesh settings yield velocities that fall within a narrower range of deviations, namely, 2–0.5 mm/s. It is also observed that the wavy behavior characteristic of $k_{\mathrm{heat}}=1$ becomes weakly expressed for $k_{\mathrm{heat}}=10$. Interestingly, in the right channel, the wavy behavior of the maximum velocity curves appears to be slightly more pronounced than in other channels when the azimuthal number of mesh elements is set to $N_z=100$.

The dynamic curves of the maximum temperature of the high- and low-voltage windings calculated for heat load coefficient $k_{heat}=10$ are presented in Figure 8. Increasing the number of elements in the azimuthal direction results in a reduction in maximum velocity, as depicted in Figure 7, and an increase in maximum temperature within the windings, as shown in Figure 8. Notably, the absolute deviation between values calculated by the finest mesh and other mesh cases for any radial number of mesh elements and azimuthal number of elements $N_z=100$ amounts to 4–6 °C. Furthermore, this deviation grows with the simulation time. Conversely, when the azimuthal number of elements is increased to $N_z=200$, the absolute deviation reduces to less than 1 °C. When the azimuthal number of elements is further increased to $N_z=300$, the deviation becomes even smaller, falling into the range of hundredths of a degree Celsius.

3.2. Study of Settings of Non-Uniform Structured Mesh Taken into Simulation of Explicit Resolution of Thermal Boundary Layer for Model #1

The results obtained in Section 3.1 demonstrate that the velocity and temperature fields exhibit a weak dependence on the number of elements along the channel radius. This phenomenon is primarily attributed to insufficient discretization in this direction near the channel walls, which can lead to inaccurate temperature predictions. Furthermore, when utilizing a geometric progression to control mesh size and capture the thermal boundary layer, an excessively large number of elements is required to resolve this region accurately. However, by adopting the mesh settings described in Section 2.6 for the “non-uniform mesh with resolved TBL”, it becomes possible to explicitly capture the thermal boundary layer with a reasonable number of elements and optimal mesh size. This approach eliminates the need for excessive refinement in regions where it is not necessary, thereby reducing computational costs and improving overall simulation efficiency.

The optimal number of elements can be achieved by considering mesh convergence for distinct parts of the simulation domain. In particular, we investigate the influence of discretization degrees on thermal boundary layers, flow cores and solid regions in both the radial and azimuthal directions. The mesh is structured, and any change to the azimuthal direction affects other parts of the simulation. Figure 9 illustrates the relative error in percentage for velocity and temperature, indicating deviations from the results obtained using the finest mesh. As depicted in this figure, the resolution of the thermal boundary layers plays a significant role in achieving accurate results. Notably, the relative error in velocity drops significantly from 5% to hundredths, as shown in Figure 9a, while temperature decreases dependent on the number of elements within the thermal boundary layer, ranging from 2.5% to hundredths (Figure 9a). Furthermore, our analysis reveals that the radial number of mesh elements has a minimal impact on temperature and velocity, as shown in Figure 9b and Figure 9c, respectively. Most importantly, we observe that the number of mesh elements in the azimuthal direction plays a crucial role in determining the magnitude of velocity, with less significance for temperature, as illustrated in Figure 9d.

Understanding the velocity and temperature profiles is crucial for identifying regions within the simulation where results diverge due to different mesh setups. Figure 10 presents the velocity and temperature profiles of the outlet of the middle channel for 4, 10 and 40 mesh elements in the thermal boundary layer, with the distributions plotted as a function of the relative width of the channel based on the width of the thermal boundary layer ${\Delta }_{TBL}$. Notably, a contradiction in calculated velocity is observed between different discretization elements from $r/{\Delta }_{HBL}=0.2$ and beyond. Specifically, the maximum deviations of velocity calculated on meshes corresponding to the number of mesh elements in the thermal boundary layer $N_{TBL}=4$ and $N_{TBL}=10$ are 3 mm/s (approximately 0.1%) and 1 mm/s (about 0.075%), respectively. It is worth noting that increasing the number of elements in the thermal boundary layer does not lead to linear growth in the accuracy of the results. Rather, temperature distributions along the outlet of the middle channel exhibit similar behavior for all considered numbers of mesh elements in the thermal boundary layer. However, it is observed that using a small number of elements (e.g., fewer than 4–6) in the thermal boundary layer leads to non-smooth temperature curves. This issue can be resolved by increasing the number of elements or using different types and orders of finite elements for the mesh.

In our previous discussion (Section 3.1), we considered the influence of the number of elements in the azimuthal direction on the simulation results. In this section, we focus on evaluating the impact of the number of elements in the azimuthal direction when the mesh setup in the radial direction has a weak effect on the results. Specifically, we examine the scenarios where the number of elements in the thermal boundary layer ($N_{TBL}=10$), flow core ($N_r=5$) and solid parts ($N_{solid}=3$) is fixed. The velocity distribution on the outlet of the middle channel is plotted as a function of the relative distance between the outlet channel width and the hydrodynamic boundary layer and is shown in Figure 11. It should be noted that the maximum difference in the calculated results between the finest mesh ($N_z=400$) and the coarser case ($N_z=20$) amounts to 4 mm/s, which corresponds to approximately 0.1% in relative units. Furthermore, when compared with other mesh setup cases, the differences are found to be less than 1.5 mm/s, equivalent to about 0.05% in relative units from the finest mesh result.

3.3. Domain Extension Study: Influence on Numerical Modeling

Excluding a domain from a numerical model allows the computational efforts and complexity of building the model to be reduced, but it might lead to significant distortion of results. In this paragraph, the influence of taking into account the simulations of the upper and lower domains of the power transformer’s tank and pipe branches in the area of the radiators is considered. The legends used in the figures of the models with different simulation domains are explained in Section 2.2 and summarized in Table 2. The thermal boundary layer (${\Delta }_{TBL}$) is resolved using 10 elements ($N_{TBL}=10$), while the remaining flow domain and solid parts are constructed with 5 elements in the radial direction ($N_r=5$). The number of elements in the azimuthal direction $N_z$ is set to 200. By examining the effects of domain exclusion on the simulation results, we aim to provide insight into the potential consequences of simplifying complex systems and to inform strategies for improving model accuracy. This investigation will help to identify the most critical domains to include in simulations, thereby optimizing computational resources and minimizing distortions in results.

The time-dependent average velocity profiles at the outlets of the left, middle and right channels for heat load coefficients $k_{heat}=1$, 5 and 10 are compared in Figure 12. The inclusion of the exterior domain and upper and bottom parts of the tank in the numerical model does not lead to significant velocity changes, as can be observed by comparing the velocity curves calculated by models #1 and #2 in Figure 12. The maximum difference between the results calculated by the previously mentioned models for $k_{heat}=10$ amounts to 2.5 mm/s, as shown in Figure 12c. However, when including the pipe branch in the area of the radiators, significant differences in velocity values are observed. The velocity curve calculated by model #3 exhibits an exponential dependence but also contains oscillatory components in time. The best convergence between the velocities calculated by models #1 and #3 is observed in the left channel, as shown in Figure 12a. In contrast, the maximum divergence of the results calculated by model #3 from those of models #1 and #2 in the left and right channels lies within a range of 20 mm/s to 35 mm/s, as depicted in Figure 12a,c. Notably, the trend of increasing velocity with growing heat load coefficients is maintained in all channels, albeit weakly in the middle channel.

The maximum temperature curves for the low- and high-voltage windings, dependent on simulation time, are presented in Figure 13. The deviations between the maximum temperatures calculated by model #3, which takes into account the most comprehensive domain of calculation, and those obtained by models #1 and #2 are provided in

Table 4. The comparison of maximum temperature values in low- and high-voltage windings calculated by model #3 with results obtained by models #1 and #2.

Heat Coefficient	Comparing Model	LV		HV
		Absolute	Relative	Absolute	Relative
$k_{heat}=1$	Model #1	0.2	0.6%	1	3%
	Model #2	0.05	0.2%	0.1	0.2%
$k_{heat}=5$	Model #1	4.25	7%	8	11%
	Model #2	1	1.6%	5	7%
$k_{heat}=10$	Model #1	11	12.5%	21	20%
	Model #2	3	3.5%	10	10%

. The results demonstrate that the influence of considering simulation domains in the numerical model of a power transformer plays a crucial role in determining temperature values, with an increasing effect observed as the heat load coefficient is raised. Comparing the maximum temperatures obtained for low- and high-voltage windings, it can be seen that they lie within different ranges: 88–100 °C and 100–120 °C, respectively. In contrast, the results for maximum temperature calculated on a uniform mesh yield values of 80–85 °C for low-voltage winding and 95–100 °C for high-voltage winding. These findings underscore the importance of accurately modeling simulation domains in power transformer design as they have a significant impact on the accuracy of temperature predictions.

The difference in the maximum temperature of low- and high-voltage windings can be explained by the velocity and temperature profiles shown in Figure 14 at the latest time step of the simulation: 7200 s. In Figure 14, temperature profiles are plotted within the solid parts of the power transformer, while velocity distributions are shown for the oil. It is worth noting that the temperature distribution exhibits slight variations depending on the chosen simulation model. Upon closer examination, the temperature distribution in the lower and upper parts of the windings reveals distinct differences between the models. The temperature profile obtained using model #1 increases monotonically along the height of the winding. In contrast, the temperature distributions calculated using model #2 and model #3 also increase along the height of the winding but eventually decrease in the upper parts of the winding. This phenomenon can be associated with the emergence of specific vortex structures within the oil flow when it is simulated by model #2 or #3. The addition of extra calculation domains leads to the appearance of new vortex structures, as observed in Figure 14c,f. A serpentine-shaped behavior of oil motion is evident in the upper parts of the channels. This phenomenon occurs at $z=0.4$ m for $k_{heat}=1$ and at $z=0.2$ m for $k_{heat}=10$, indicating that the location of turbulent flow generation is strongly dependent on the intensity of the power transformer thermal mode and may appear anywhere within the domain of the channel. Notably, model #3 exhibits a greater number of vortex structures compared to other models. Furthermore, these vortices appear in the upper nozzle and propagate towards the pipe branch, leading to backflow and reducing the effective cross-sectional area for liquid pumping.

To facilitate a detailed examination of the velocity profiles depicted in Figure 14 for models #2 and #3, the upper parts of the tank are zoomed in on and replotted in Figure 15. The velocity distribution calculated by model #2 reveals that the transformer oil flows into the branch only through half the cross-sectional area of the upper nozzle, with the remaining portion exhibiting an opposite direction of oil motion. A similar behavior is observed for model #3 when simulated at a heat load coefficient of $k_{heat}=10$, as shown in Figure 15c. In this case, fully developed vortices form within the nozzle, indicating a complex and turbulent flow regime. These findings suggest that the velocity profiles and vortex formation within the power transformer are highly dependent on the heat load coefficient and simulation model employed. The detailed examination of these phenomena provides valuable insights into the complex fluid dynamics at play within the power transformer and highlights the importance of accurate modeling and simulation techniques for predicting flow behavior under various operating conditions.

A quantitative assessment of velocity profiles can be made by comparing the results calculated by the three models for a range of heat load coefficients ($k_{heat}$ = 1–10) at the inlets and outlets of the channels, as shown in Figure 16. The curves of velocity magnitude obtained by model #1 exhibit a symmetrical distribution at all inlets of the channels for all considered heat load coefficients. These curves display a parabolic shape for low heat load coefficients and transform into an M-shape for higher coefficients. The magnitude of velocity calculated by model #1 becomes non-symmetrically distributed along the width of the outlet channels, with different peak velocities near the walls. For example, in Figure 16d, the curve calculated for $k_{heat}=10$ at the left channel outlet shows approximately twice the velocity values near the right wall compared to the left one. Moreover, the velocity curves exhibit an even more pronounced M-shape characteristic at the middle channel outlet compared to the inlet. This phenomenon is attributed to the growing influence of natural convection due to increasing temperature along the height of the windings. Asymmetric profiles of velocity are observed when adding an exterior domain to the tank, as seen in models #2 and #3, for any heat load coefficients. It is noteworthy that these asymmetric profiles are more pronounced in model #3. Special attention should be paid to the curves of velocity magnitude shown in Figure 16f. In this case, the curves undergo a V-shape transition near the center of the channel, indicating that the magnitude must only be positive and suggesting the presence of backflows.

4. Discussion

The results obtained by calculation of model #1 using a non-uniform mesh without explicitly resolving the thermal boundary layers show underestimated values of temperature and overestimated values of velocity compared to results calculated by a non-uniform mesh resolving both thermal and hydrodynamic boundary layers (results from Section 3.2). A comparison of temperature and velocity at the last time step of the simulation using different meshes is presented in Table 5. The temperature and velocity values presented in this table were calculated using the finest sizes of the “non-uniform mesh” and “non-uniform mesh with resolved TBL” settings, which correspond to conditions where the relative error is less than 1%. Based on the results shown in Figure 9, the “non-uniform mesh with resolved TBL” settings for the simulations are characterized by the following parameters and dependencies:

• The number of elements in the thermal boundary layer is $N_{TBL}=10$. This parameter has a significant influence on both thermal and hydrodynamic designs, as reflected in the accurate temperature values obtained.
• The number of elements needed to resolve core flow is $N_c=5$. The parameter reflects a weak influence on thermal design and does not introduce an impact on the hydrodynamic one.
• The number of elements in the solid parts is $N_{solid}=3$ in the radial direction. The parameter reflects a weak influence on both thermal and hydrodynamic calculations.
• The number of elements in the solid parts is $N_z=200$ in the azimuthal direction. This parameter has a significant influence on hydrodynamic calculation.

Table 5. The comparison of maximum temperature in low- and high-voltage windings and velocity at the inlets of the left, middle and right channels. The data are calculated by model #1 using the finest uniform and non-uniform meshes.

Type of Mesh	Heat Coefficient	Temperature, °C		Velocity, mm/s
		HV	LV	Left	Middle	Right
Mesh # 1 from Section 3.1	$k_{heat}=1$	31	29	12	20	9
	$k_{heat}=10$	100	80	38	50	20
Mesh # 2 from Section 3.2	$k_{heat}=1$	33.6	30	10	15	5
	$k_{heat}=10$	120	100	35	40	17

Table 5. The comparison of maximum temperature in low- and high-voltage windings and velocity at the inlets of the left, middle and right channels. The data are calculated by model #1 using the finest uniform and non-uniform meshes.

Type of Mesh	Heat Coefficient	Temperature, °C		Velocity, mm/s
		HV	LV	Left	Middle	Right
Mesh # 1 from Section 3.1	$k_{heat}=1$	31	29	12	20	9
	$k_{heat}=10$	100	80	38	50	20
Mesh # 2 from Section 3.2	$k_{heat}=1$	33.6	30	10	15	5
	$k_{heat}=10$	120	100	35	40	17

The relative differences for the maximum values of the temperature and velocity calculations in the case of optimal mesh parameter settings do not exceed 0.4% and 5%, respectively, for the non-uniform mesh. These results were obtained relative to the finest non-uniform mesh described in Section 3.2. They are presented in Figure 17. The heat load coefficients considered cover underloaded and overloaded modes, indicating that the recommended mesh parameters can be applied to a wide range of problem scenarios involving the simulation of power transformers.

The resolution of thermal boundary layers is crucial for providing accurate simulations and predicting desired results. Even with a mesh containing more elements, as described in Section 3.1, compared to some of the mesh parameter cases outlined in Section 3.2, significant distortions in results can occur. Therefore, the evaluation of thermal and hydrodynamic boundary layers is an essential part of preliminary simulation. To conduct this evaluation, the average velocity in the channel is required, which can be obtained through simulations under coarse meshes or by using the following well-established expression, which is a viable option:

$$U_0=\sqrt{g\beta (T_w-T_{amb})L}$$

(13)

where $T_w$ is the maximum temperature of the winding and $T_{amb}$ is the temperature of the inlet oil. Alternatively, a preliminary simulation can be run on a coarse mesh to estimate the velocity. It should be noted that using (13) may lead to an overestimated velocity. Therefore, running a preliminary simulation on coarse mesh is generally preferable. If this method is used, it is recommended that a coarser mesh based on the initial results obtained from expression (13) is prepared.

The use of a non-uniform mesh with an explicitly resolved TBL allows for a significant reduction in the number of elements required for model #1 (the reduction is approximately 20%). This approach also enables the layer to be resolved, whereas a non-uniform mesh with geometrically growing elements necessitates more elements and fails to resolve the layer in the case of an insufficient mesh. As previously noted, sufficient resolution of the layer has a substantial impact on the identification of hot spots.

Implementing a mesh technique to resolve layers near walls leads to a sharp increase in element size. In such cases, it is essential to employ conservative numerical methods, such as the finite volume method, and special polynomials to describe the desired unknown functions, or to select an optimal growth coefficient for the elements within the layer of interest.

The configuration and composition of the calculated domains in a numerical model of power transformers play an important role. Excluding areas where oil does not directly interact with windings or a steel core from calculations can lead to results being distorted by up to 20%. The positions of hot spots change when additional domains are introduced into the models. In these extra domains, oil flow paths form and turbulent flows may occur, affecting the velocity profile at the channel inlet and, consequently, the flow behavior in the channel itself.

The nature of liquid flows can significantly affect the intensity of heat removal from cooled surfaces and lead to increased calculation complexity for temperature distribution. Calculated velocity profiles at the inlets (Figure 16) become more non-symmetrically distributed as additional domains are added to the calculation model, primarily due to the formation of vortices at the top and bottom parts of the tank (Figure 14 and Figure 15). Additionally, velocity profiles at the outlets exhibit more non-symmetrical curves compared with inlet profiles (Figure 16).

The increasing non-symmetrical behavior of velocity can be illustrated by a detailed explanation of the serpentine-shaped distribution of velocity along the length of the channels, as shown in Figure 14c,d. To accurately assess this phenomenon, it is essential to separate natural convection into two distinct types: vertical and horizontal. Vertical convection occurs due to temperature differences along the height of windings, whereas horizontal convection arises from temperature differences between solid parts. When the calculation area is expanded, the hot spot temperature for both high- and low-voltage windings changes their coordinates. In model #1, the hot spot temperature is located at the highest point of the winding height. In model #2, this point is shifted slightly lower, but still remains at the top. However, for model #3, the maximum temperature point moves to a slightly upper-middle part of the windings. The shift of the hot spot temperature leads to an increase in the temperature gradient in the vertical direction, which enhances the effect of vertical natural convection. Furthermore, the difference in temperatures in the horizontal direction between solid parts grows as well when shifting the hot spot temperature to the model winding part. This introduces an increasing impact of horizontal convection. These results indicate that natural convection plays a significant role as the height of windings increases. Therefore, when assessing natural convection intensity, it is crucial to consider both vertical and horizontal convection, using similarity numbers and characteristic dimensions that account for these factors.

The time spread can be attributed to varying intensities of physical phenomena over time for different load factors. The simulation times of the models under consideration are as follows:

• For model #1, 1–2 min.
• For model #2, 2–4 h.
• For model #3, 10–13 h.

Specifically, calculations require more time when the heat load coefficient equals 5.

It should be noted that oil flows forming at the top and bottom parts of the transformer tank can cause turbulent flows in three-dimensional space that propagate through the channels. This phenomenon was not explored in this study.

The presented results and recommendations for building the mesh can be used for the simulation of similar systems characterized by a Grashof number in the range of $5·{10}^8$ to $5·{10}^{11}$, as expressed by Equation (14). Additionally, these results and recommendations apply to systems with a Reynolds number ($Re$) of 100–1000, where $Re=\rho UL/\mu$.

$$Gr={\displaystyle \frac{g\beta \left(T_w-T_{amb}\right)L^3}{{\nu }^2}}.$$

(14)

5. Conclusions

In this paper, a mesh sensitivity analysis for the simulation of natural convection in an oil-immersed transformer with an external cooling system has been provided. The study emphasizes the importance of resolving the thermal boundary layer using appropriate mesh settings. It is shown that, in order to resolve the thermal boundary layer, it is necessary to conduct a preliminary assessment before carrying out the modeling. Several procedures are described for conducting the preliminary assessment. Recommendations for building a mesh based on the obtained results are provided. A brief comparison of the computation time of the numerical model with different mesh settings is presented. In further research, it would be of great interest to investigate the application of adaptive meshing techniques to tailor the mesh configuration for various thermal modes, thereby minimizing computational expenses and maximizing simulation efficiency. However, when employing adaptive meshing techniques, it is essential to ensure that the additional computational effort required for mesh assembly does not exceed the overall simulation runtime, thereby offsetting any potential benefits.

Three suggested numerical models are described in detail from a mathematical point of view and developed using new combinations of boundary conditions to reduce the simulation domains. The primary difference between these models lies in their computational domain composition, which affects the accuracy of temperature distribution and velocity profiles. The first model focuses on conjugate heat transfer due to the direct interaction between liquid and solid components, while the second and third models account for vortices and complex flow structures outside the channels and flows appearing in pipe branches, respectively. The analysis reveals that the emergence of specific vortex structures within the oil flow when simulated by models #2 or #3 leads to alteration of hot spot positions. The addition of extra calculation domains leads to the formation of vortices, which affect the velocity profile at the channel inlet and, consequently, the flow behavior in the channel itself. The work also provides the required calculation time for the models under consideration. These studies and results are very important when selecting the necessary models for integration into a digital twin system.

The findings presented in this study can be applied not only to digital twin technology but also to the development of design conclusions and remarks, as demonstrated in [1]. Furthermore, numerous scientific and engineering problems related to natural convection in power transformers require parametric studies to be conducted, which can be time consuming. In a similar vein, Ref. [39] presents a parametric study based on a simplified numerical model of a power transformer to examine the correlations between the Prandtl number, Reynolds number and hot spot factor. The results of this study provide valuable insights into the underlying relationships governing these parameters.

This research article presents a comprehensive analysis of mesh settings for numerical models of conjugate heat transfer in oil-immersed transformers. The study offers original insights into selecting optimal mesh configurations that efficiently balance computational resource demands while ensuring reliable model performance. The findings contribute to enhancing simulation practices by highlighting key considerations for mesh design that support effective thermal management analysis in transformer applications.