Review of Modeling Approaches for Conjugate Heat Transfer Processes in Oil-Immersed Transformers

Abstract

This review addresses the modeling approaches for heat transfer processes in oil-immersed transformer. Electromagnetic, thermal, and hydrodynamic thermal fields are identified as the most critical aspects in describing the state of the transformer. The paper compares the implementation complexity, calculation time, and details of the results for different approaches to creating a mathematical model, such as circuit-based models and finite element and finite volume methods. Examples of successful model implementation are provided, along with the features of oil-immersed transformer modeling. In addition, the review considers the strengths and limitations of the considered models in relation to creating a digital twin of a transformer. The review concludes that it is not feasible to create a universal model that accounts for all the features of physical processes in an oil-immersed transformer, operates in real time for a digital twin, and provides the required accuracy at the same time. The conducted research shows that joint modeling of electromagnetic and thermal processes, reducing the dimensionality of models, provides the most comprehensive solution to the problem.

1. Introduction

Power transformer is the key element of distribution networks, where the most common type is an oil-immersed transformer. Its technical state determines the reliability and safety of the entire power system. Transformer defects or deviations in their operating parameters can be caused due to various reasons, such as humidification, oil and insulation contamination, and the transformer’s mismatch in terms of its electrodynamic resistance to short-circuit currents. Another group of defects, mechanical, is introduced by excessive influences: short-circuit currents, thermal overloads, electrical overvoltages, and overexcitation of magnetic systems. The third group is design and manufacturing defects, etc. [1,2]. It should be noted that modern operating experience shows that the main transformer defects are associated with a lack of information about their current technical state. For this reason, monitoring the main parameters of transformer functioning plays an important role [1,2,3,4,5]. Additionally, the problem of transition to maintenance based on the actual state of the transformers remains relevant. This requires technical diagnostics by obtaining the most complete and reliable diagnostic data synchronously or within a minimum amount of time. A solution to this problem might be the use of digital twin (DT) technology for modeling equipment. It is described in detail by the authors of the article [6]. This direction combines systems of equipment operating parameters monitoring and systems of their processing and analysis. Digital twin technology can be implemented using either artificial intelligence (AI) methods or models based on the physical description of the object. AI-based models require a large volume and completeness of the data, as the observed parameters of deep retrospective allow for correct generalization and interpretation of the results. In work [7], the authors describe the problem of power transformers’ technical state assessment based on data preprocessing and machine learning. The drawbacks of DT based on physical modeling include the complexity of mathematically describing the physical processes of various nature occurring in transformers. However, a DT implemented using mathematical models of physical processes is more robust to random events and requires fewer operating parameters of equipment for its implementation compared to AI-based models. Therefore, it is advisable to consider an approach to implement DT technology of an oil-immersed transformer based on a mathematical description of its physical processes.

Nowadays, quite a lot of works are devoted to the implementation of DT based on physical models [8,9,10,11,12,13,14]. Work [8] describes the effectiveness of using virtual sensors that expand the capabilities of monitoring the state of transformers. The authors note that installing physical sensors is quite a difficult task due to the problems related to the accessibility of installation sites caused by complex equipment design. They also draw attention to the fact that local sensor readings do not provide detailed information about the current distribution of the controlled volume quantity. However, this is necessary to predict the state of the object more accurately. For this reason, the authors propose to use virtual sensors based on the multiphysics modeling of physical fields in a transformer, namely, temperature and hydrodynamics. Similar multiphysics modeling is also used in the works [9,10,11], where the authors use multiphysics modeling via numerical methods. Finite element and finite volume methods are used to predict the temperature change of the transformer more accurately. Another concept is to use simpler equivalent circuits as a physical model [12,13,14]. The authors of the work [12] create a DT of a transformer based on the mechanical equivalent circuit. It makes it possible to detect various emergency situations associated with the mechanical deformation of parts of the transformer. The authors of work [13] represent a DT, based on an electrical equivalent circuit of a transformer. Such a design of the DT makes it possible to obtain the parameters of the high-voltage winding without actual measurements, considering information only about the parameters of the low-voltage winding. The use of an electrical equivalent circuit for the purpose of monitoring unmeasured (additional) parameters, such as the magnetizing current of a transformer, is also considered in [14]. Summarizing the research data, it can be argued that physical models of transformers are a key element in creating a full-fledged DT of a transformer.

It is evident that the state of the transformer is influenced by a complex interplay of physical processes. One such process is the conversion of electromagnetic energy from the primary circuit to the secondary one, which can be described using electromagnetic theory. The electromagnetic field’s action in a physical medium results in the release of thermal energy and the presence of electromagnetic forces. It is also crucial to consider the impact of the transformer’s thermal state on its operation, as this alters its physical properties, such as electrical and thermal conductivity. To accurately describe the temperature state of the transformer, convection in the cooling fluid (transformer oil) must be taken into account. The electromagnetic forces can cause deformation of the transformer parts and vibration, which is accompanied by the generation of acoustic waves. Additionally, chemical processes like insulation aging and copper corrosion affect the transformer’s condition, and the interactions between these processes are depicted in a flow chart in Figure 1.

As it can be seen, for a complete physical description of all processes in a transformer, a rather complex multiphysics model is required. It can be implemented in various ways, including the use of analytical dependencies together with a circuit theory-based model. Another way is to use numerical solution of distribution problem and interaction of physical fields. The final mathematical model depends on the required physical processes and the detalization level of results. It is also worth emphasizing that the description of approaches to modeling all physical processes present in a transformer goes beyond the scope of one scientific work. For this reason, the presented review is limited to considering the modeling of the electromagnetic problem and the problem of coupled heat transfer as the main phenomena describing the operation and state of an oil-immersed transformer. This choice can be justified by a significant number of works [6,7,15,16,17,18,19,20,21,22,23], dedicated to assessing the state of the transformer and aging of insulation, where the temperature distribution and the values of the hottest points are of key importance. As illustrative example, the work [17] is devoted to the study of the degradation of insulating paper under the action of a non-uniform temperature field. The authors use a second-order kinematic model to simulate the aging of insulating paper. Its difference from the traditional kinematic model is that it takes into account the uneven heating of the transformer windings, so this model can be considered more accurate from a physical point of view than the traditional one. As conclusions, the research team notes the following: the moisture content of the paper is different along the height of the winding and opposite to the direction of the temperature gradient, and the model verification results show the superiority of the proposed model over the traditional one. Similar results are obtained in work [18]. The authors argue that most models do not take into account some of the phenomena that are present in a real transformer. For instance, different parts of the transformer have different degrees of aging due to the high temperatures of the windings and its non-uniform distribution. As the result, the concentration of furfurals in windings is also non-uniform. However, due to oil circulation as schematically shown in Figure 2, the concentration of furfurals becomes uniform over the volume of the tank, so the content of furfurals reflects the average state of the insulation. But the furfural generation rate and temperature distribution are non-uniform along the length of the winding (L) as shown in Figure 2. For that reason, it is crucial to take into account in the mathematical model in order to precisely predict insulation aging.

This review focuses on the most commonly used methods of mathematical modeling for electromagnetic, thermal, and hydrodynamic processes in core-type oil-immersed power transformers with two-winding configurations based on natural convection cooling. The accurate modeling of heat transfer processes in oil-immersed transformers is crucial for ensuring their reliable operation, especially when considering state monitoring through robust real-time digital twin technology. The aim of this review is to address the challenge of developing a comprehensive digital twin of a transformer by analyzing various modeling approaches. The complexity arises from the coupled interaction of physical processes occurring within an oil-immersed transformer. To create an effective digital twin, it is necessary to develop simplified mathematical models that sufficiently capture the behavior of these processes while maintaining accuracy. The review compares different methods, such as circuit-based models, finite element analysis, and finite volume methods, considering factors like implementation complexity, calculation time, accuracy, and practical considerations. The physical fields of electromagnetics, thermal dynamics, and hydrodynamics are considered critical aspects of transformer modeling. By examining the strengths and limitations of these modeling methods, the review assesses their ability to create an accurate and robust digital twin for core-type oil-immersed power transformers. The remaining sections of this paper are structured as follows. In Section 2, we provide a brief mathematical description of each considered mathematical method and discuss their advantages and disadvantages for solving specific problems, based on their mathematical features and formulations. Moving on to Section 3, we offer a comprehensive review of contemporary methods for electromagnetic simulation of transformers, considering the unique aspects of coil and steel part modeling. In Section 4, we assess the impact of auxiliary structural components on the results of transformer simulation. Section 5 addresses the simulation of transformer temperature conditions by solving the conjugate heat transfer task. Section 6 contains remarks drawn from the reviewed simulation methods and their distinct characteristics. Finally, in Section 7, the paper is concluded.

2. Review of Mathematical Modeling

Mathematical modeling represents the object of study in mathematical abstractions based on physical laws, using a number of assumptions. Various methods and approaches used in mathematical modeling have been developed. Power transformer modeling can be implemented by these methods. For instance, conjugate heat transfer models are built using numerical, analytical, or circuit theory methods. It is important to note that modern numerical models of thermal processes in a transformer often combine the advantages of several methods, called hybrid models. A classification is proposed as illustrated in Figure 3, based on the approach to describing partial differential equations for various types of physical fields associated with the processes of modeling phenomena in oil-immersed transformers.

The effective application of mathematical methods and tools in practical settings poses a complex challenge due to the multitude of available variants. Therefore, studies that focus on assessing the complexity and usability of these tools are crucial. Several works have been dedicated to exploring this topic [24,25,26]. In [24], a review of mathematical methods for using real-time simulation is given, but the paper does not consider issues related to the creation of complex numerical models. The authors consider the possibilities of implementing mathematical models in real-time simulation to estimate the electrical parameters of a transformer. In the study [25], the author conducts a comprehensive analysis of power transformer cost optimization methods. The focus is on exploring the main trends and directions in power transformer design within the context of optimization problems. The study concludes that novel optimization methods for transformers take into account more complex transformer models and incorporate recent advancements in optimization theory. There is a clear trend towards the evolution of not only transformer models but also optimization aspects in transformer design.

The selection of numerical modeling methods for research on power transformers is influenced by various factors, including the specific research question, the availability of computational resources, the expertise of the researchers involved, the flexibility and availability of software tools, and the traditional methods used in the field. This flexibility results in various classifications of mathematical models from different perspectives, such as spatial regions or model formalization approaches. In this review section, we will discuss the characteristics of different mathematical approaches to simulating physical phenomena in power transformers.

2.1. Analytical Approaches

Analytical approaches involve obtaining solutions to differential equations that describe specific physical phenomena. The use of analytical methods has a significant advantage in terms of ease of application, as the obtained dependencies can be readily integrated into control systems and software applications for engineering calculations. This is particularly important for the diagnostic of power transformers, where high calculation performance is essential. However, a major drawback of analytical models is the difficulty in developing a universal procedure for creating an analytical solution that can accurately describe complex physical phenomena dependent on numerous variables. However, despite these limitations, several scientific studies have successfully employed analytical approaches to simulate power transformers, yielding results that demonstrate sufficient accuracy and rapid calculations.

In the study [27], the authors address the problem of formulating the short circuit impedance (SCI) analytically. They focus on optimizing core-type power transformers and autotransformers by considering the value of the SCI as a key parameter. The SCI plays a critical role in power transformer design, as it influences the proper functioning of the electrical grid and helps mitigate short-circuit stresses in the winding. The authors extensively examine the challenges associated with implementing geometric programming for core-type power transformers and autotransformers. To tackle these challenges, the authors propose a solution method that combines geometric programming with the Branch-and-Bound method. This approach effectively addresses the problem of SCI for large core-type power autotransformers. In [28], the results of the analytical method, proposed by the authors for calculating the mutual inductance of electric coils, are compared with 2D and 3D models based on finite elements. The difference between the results does not exceed 1%, and the calculation time for the analytical model takes 0.36–0.75 s, for two-dimensional models, 9–9.6 s, and for three-dimensional models, 8–10 min. Analytical approaches are well suited for preliminary calculations in the design of transformers. It is shown in [29] that the reluctance calculated by approximate analytical expressions differ from those measured by 30–32%, while the calculations based on the FEM (finite element method) by 4.2%. In the same work, the analytical methods agree better with the measurements in terms of the values of mutual induction. Damaged conditions can be determined using the analytical model from [30] by instantaneously modifying the currents drawn in the sound condition. The papers [31,32] present a methodology for designing high-frequency transformers based solely on analytical methods. This technique enables high-precision analysis of magnetic and electrical losses in the coils and steel of the transformer. The accuracy of the model is evaluated by comparing experimental results. While analytical methods can significantly expedite calculation times, they may not provide the desired level of accuracy for certain applications. As such, analytical approaches are often used to streamline calculations in other mathematical models by replacing complex components with simpler, more efficient alternatives.

2.2. Models Based on Circuit Theory

The method of representing a transformer with the use of equivalent circuits is widely used in practice to define heat sources and other important parameters for real-time simulations. The essence of the method lies in the abstraction of the representation of an object in the form of a connection of chain elements to describe a physical phenomenon. With the help of this method, it is possible to describe the processes of electromagnetic energy conversion, analyze the temperature modes of operation, and predict the mode of operation of the cooling system. The application of this approach can significantly reduce the computation time, the requirements for the computing station, and the necessary qualifications for implementing the model based on such methods. The disadvantages of this method are that it has low accuracy, especially for complex non-traditional transformer designs, the non-linearity of the system is poorly taken into account, and there is no possibility of analyzing the local parameters of the system.

Let us consider the mathematical principle of replacing partial differential equations with the help of equivalent circuits, using the example of differential Equation (1) for a two-dimensional magnetic field written relative to the vector magnetic potential $\mathbf{A}$. In Equation (1), $\sigma$ and $\mu$ are scalar quantities describing the physical properties of space, while $\mathbf{J}$ is the current density vector, given in space (Figure 4a):

$$\sigma \mu \frac{\partial \mathbf{A}}{\partial t}+{\nabla }^2\mathbf{A}=\mu \mathbf{J},$$

(1)

Assume that the sampling region is finite and discrete. Discrete point indices in the x direction correspond to i and in the y direction are expressed using j. Using the Stokes theorem and the finite-difference approximation of derivatives, we can rewrite Equation (1) for the node $(i,j)$ in the following form:

$$\begin{array}{c}-\frac{\sigma \Delta x\Delta y}{L}·\frac{\partial {\Phi }_{i,j}}{\partial t}+{\Phi }_{i+1,j}\left[\frac{\Delta y}{\Delta xL\mu }+v_x\frac{\sigma \Delta y}{2L}\right]+{\Phi }_{i-1,j}\left[\frac{\Delta y}{\Delta xL\mu }{-v}_x\frac{\sigma \Delta y}{2L}\right]+{\Phi }_{i,j+1}\left[\frac{\Delta x}{\Delta yL\mu }{+v}_y\frac{\sigma \Delta x}{2L}\right]+\hfill \\ \hfill {\Phi }_{i,j-1}\left[\frac{\Delta x}{\Delta yL\mu }{-v}_y\frac{\sigma \Delta x}{2L}\right]=F_{i,j}\end{array}$$

(2)

and represent the constant coefficients as resistances (Figure 4b). In this case, these will be reluctance to magnetic flux. The work [33] considers the process of deriving detailed equivalent circuits in more detail for electromagnetic, temperature and hydrodynamic fields.

Using the example of [34], the principle of drawing up a detailed equivalent circuit for modeling an electromagnetic field can be demonstrated. The main sources of generation and trajectories of the magnetic field closure in the parts of the transformer are singled out. In this case, the winding flux of three-phase coils, leakage fluxes and zero-sequence fluxes are considered (Figure 5). For each flow trajectory, depending on its geometrical parameters and the physical properties of the region, through which the magnetic field passes, the values of equivalent reluctances are calculated. In conclusion, a magnetic equivalent circuit is drawn up from magnetic field generators and reluctances. The model in [34] is based on a field programmable gate array (FPGA) with a high-resolution non-linear magnetic equivalent circuit (HR-MEC). The article [34] highlights the following advantages of this approach:

• Accounting for the main flow paths in the transformer core using detailed magnetic field equivalent circuit;
• Modeling the hysteresis of a transformer core using Preisach theory;
• Modeling eddy currents with the use of a frequency-dependent equivalent circuit;
• Reduced model latency through the use of full hardware parallelism and pipelining on the FPGA, and thus, smaller temporary steps can be used to capture high-frequency transients.

Equivalent circuits can be drawn up with details up to the structural parts of the transformer. In [35], T-shaped equivalent circuits are considered for calculating the electrical and magnetic parameters of a transformer. In T-shaped equivalent circuits, magnetic and electrical resistances are calculated for the primary and secondary coils, and for their magnetic flux linkage. The solution of such a circuit is represented by the integral parameters of the fields. In [36], a comparison of two equivalent circuits of a transformer (a “T” type circuit and a “$\Pi$” type circuit) is considered from the point of view of calculating the inrush current of the transformer. The authors note that both circuits have sufficient convergence in the steady state, but the inrush current obtained using the “$\Pi$” circuit is closer to the experimental values. This discrepancy is caused by the fact that in the “T” circuit, the inductance of the primary coil is before the inductance of the magnetizing branch. In this case, the inductance of the primary coil incorrectly limits the flow of current through the magnetization branch when the magnetic circuit is saturated. Work [37] is an example of using detailed equivalent circuits, where the level of detail is determined by each individual turn of the transformer coils. This level of detail is necessary to account for the phenomenon of partial discharge, for which the electrical parameters of each turn are determined, and the occurrence of a partial discharge is introduced by an additional voltage source. The study [38] proposes a technique for compiling equivalent circuits based on traditional transformer equivalent circuits to determine the turn-to-turn short circuit. Equivalent circuits are also used in the calculation of the temperature conditions of the transformer. In studies [39,40,41,42], the thermoelectric analogy method based on thermal equivalent circuits is used. The advantage of this approach is a shorter calculation time due to a simplified mathematical representation compared to partial differential models.

2.3. Numerical Models

The physical phenomena of a transformer are in most cases described by partial differential equations. The search for an analytical solution of these equations is a complex mathematical problem, and for a number of cases, there are no approaches to their solution analytically. Numerical methods make it possible to obtain approximate solutions of such equations with sufficiently high accuracy. The main idea of numerical models is to move from a partial differential equation to a system of linear equations written for a discrete domain. The accuracy of the results should theoretically increase with increasing discretization, but this increases the number of equations that need to be solved. Numerical methods are quite easy to formalize using software algorithms. Therefore, numerical methods are very common in modern engineering programs for calculation (CAE—Computer Aided Engineering).

The development of computer technology has significantly influenced the development of numerical methods. The most commonly used methods in practice are finite volumes, finite elements, finite differences, boundary elements, and meshless methods. Each of the methods has its own advantages and disadvantages. For example, the discretization process requires less time when using the finite difference method compared to other methods. However, the accuracy of the results is significantly reduced when calculating a space with complex geometry. In [43], the author considers the finite difference, finite volume, and finite element methods for solving physical problems. As a conclusion, it is reported that the methods of finite differences and volumes have a high stability of solutions to dynamic problems with homogeneous physical properties in space. The finite element method, on the contrary, shows high efficiency in solving problems with heterogeneous physical properties and a complex computational domain. Therefore, the author notes that it is very important to combine these methods to obtain the best results. As part of the review, the finite element and finite volume methods are considered, being the most common in modern research related to the numerical simulation of a transformer.

Let us consider the mathematical interpretation of the methods using the example of a partial differential equation for an abstract stationary physical field (3):

$$\nabla ·\left(\nabla u\right)=F(u,\mathbf{x}),$$

(3)

where u is the scalar value dependent on the spatial distribution, F is the source of the scalar field, and $\mathbf{x}$ is the vector of spatial coordinates belonging to space $\Omega$.

2.3.1. Finite Element Method

In the finite element method (FEM), Equation (3) is transformed into an integral form with multiplication by the test function

$${\int }_{\Omega }\nabla ·(\nabla u)\varphi dV={\int }_{\Omega }F\varphi dV,$$

(4)

where $\varphi$ is the test function.

Applying the Gauss theorem and making some transformations, the following formula can be obtained:

$${\int }_{\partial \Omega }\left(\nabla u\right)·\mathbf{n}\varphi dS-{\int }_{\Omega }\left(\nabla \mathrm{u}\right)·\nabla \varphi dV={\int }_{\Omega }F\varphi dV,$$

(5)

where $\partial \Omega$ is the boundary of the computational domain and S is the surface area of the computational domain.

Expression (5) is called a weak formulation, which must be satisfied for the entire variety of test functions. Test functions can be represented by various functions, but most often they are polynomials. In this case, discretization is due to the choice of a finite number of test functions.

Since in the end we must obtain an approximate solution for the desired value u, we can state that

$$u\approx u_h\left(\mathbf{x}\right)=\sum _i{u}_i{\psi }_i\left(\mathbf{x}\right),$$

(6)

where $u_h$ is the approximate solution, and ${\psi }_i$ is the basis function of the finite element.

Using expression (6), the Equation (5) could be transformed as

$$\sum _i{u}_i{\int }_{\Omega }\nabla {\psi }_i·\nabla {\psi }_{j}dV+\sum _i{\int }_{\partial \Omega }{u}_i\nabla {\psi }_i·\mathbf{n}{\psi }_{j}dS={\int }_{\Omega }F\left(\sum _i{u}_i{\psi }_i\right){\psi }_{j}dV,$$

(7)

The resulting equation is the basis for compiling a system of linear algebraic equations, with the help of which the required physical quantity is found.

Setting up the computational mesh is a key step in creating a model. The authors of [44] pay special attention to this, having considered the mesh tuning process using the example of the problem of thermal convection in the cooling tank of a transformer. The mesh consists of about 1 million finite elements of various types for the 7.5 degree sector of the transformer. The calculation takes from 90 to 120 min when using a workstation with two 2.66 GHz processors and 48 GB of RAM (Random Access Memory). In [45], a similar study is carried out using three-dimensional and two-dimensional models, based on the FEM with the number of elements 17304800 and 988906, respectively. A three-dimensional model requires hundreds of GB of RAM for calculations, and a two-dimensional model requires dozens. In [46], the possibility of reducing the computing power from 30 GB to 4 GB of RAM for calculating finite element models for calculating the three-dimensional magnetic field of a transformer is considered, which leads to a decrease in the calculation time from 43 h 44 min to 3 h 32 min.

Among the advantages of finite elements, one can distinguish the accuracy of modeling complex object geometries, the significant development of the method, and the ability to take into account local effects.

2.3.2. Finite Volume Method

The formulation for solving Equation (3) by the finite volume method (FVM) is a special case of Equation (5), where the test functions are represented as a constant equal to 1, then

$${\int }_{\partial \Omega }\left(\nabla u\right)·\mathbf{n}dS={\int }_{\Omega }FdV.$$

(8)

It can be seen here that discretization is performed by choosing a finite number of control volumes, for which expression (8) must be satisfied. It can also be argued that expression (8) is a conservation law for the volume under consideration, which is the main difference from formulation (5), where the conservation law is valid for the entire computational domain, and not for a discrete element.

When using this method, special attention must be paid to the type of used mesh. In [47], a similar problem of thermal convection, considered in the finite element section in [45], is solved using the finite volume method on a mesh, consisting of about ${10}^6$ nodes and $2.8\times {10}^6$ elements. In [48], the authors consider, how the use of different meshes affects the calculated values. In the study, different options are compared respectively to the resulting discretization error and the accuracy of the gradient reconstruction. The authors report that:

• The accuracy of the gradient reconstruction deteriorates on meshes with high ratio parties;
• The lack of mesh regularity significantly affects the approximation error;
• The discretization error is almost insensitive to the regularity of the mesh on triangular elements, but the sensitivity is greater when using quadrilateral elements;
• The dependence of the solution accuracy on the characteristics of the mesh has a complex dependence.

2.3.3. Comparison of Considered Methods

In work [47], the multi-physics modeling of an oil-immersed power transformer was conducted using the finite element method by Ansys Maxwell and the finite volume method realized in Ansys Fluent for the three-dimensional case. The authors utilized finite element analysis to calculate the transformer winding losses and employed the finite volume method to build a thermal model for forecasting the hot-spot temperature. In [49], a comparison is made between the finite volume method and the finite element method using the example of hydrodynamic problems. The study finds that software packages using the finite element method are more sensitive to the type of mesh and its quality, which is a result of the discretization formulation used. In contrast, the finite volume method shows a faster solution speed in problems of hydrodynamics, with a difference of about five times compared to the finite element method. This suggests that the finite element method may struggle with tasks that depend on time. However, the advantage of the finite element method lies in its flexibility in using various mesh types and the continuity of the resulting solution, as the solution is sought at each point of the computational mesh, rather than in discrete volumes. This difference can be crucial in certain problem classes, such as electrodynamics. For the sake of brevity and clarity, the main advantages, disadvantages, and areas of application of the mathematical methods discussed here are summarized in

Table 1. Comparison of mathematical methods.

	FEM	FVM	Equivalent Circuit	Analytical
Application	Calculation of the magnetic field, mechanics, heat transfer tasks	Hydrodynamics, convection	Digital twins, control systems, diagnostics, preliminary and auxiliary calculations	Transformer design, preliminary calculations, auxiliary calculations
Programs	Comsol Multiphysics, ANSYS Maxwell, Elmer, Quick Field, Agros2D, CalculiX, FEniCS, Hermes, GetFEM++, GetDP, JMAG, ADINA, EMWorks, SIMULIA, NX Nastran, Altair Flux, INTEGRATED Engineering Software, TRAFALO, Code_Saturne	OpenFOAM, ANSYS Fluent, OpenFVM, Clawpack, Gambit	Simulink, MatLab, SciLab, EMTP, PSCAD/EMTDC, Python libraries, almost any environment that allows to implement mathematical calculations	Proprietary code, MathCad, in commercial programs to accelerate a number of procedures
Calculation time	Long (strongly depends on the PC performance)	Long (strongly depends on the PC performance)	Fast	Extremely fast
Embeddability in software	Medium	Medium	High	High
Accuracy	High	High	Medium, requires clarification	Depends on the object
PC Requirement	High	High	Below medium	Low
Works in which this method is used	[45,46]	[47,50,51]	[33,34,35,36,37,38,39,40,41,42]	[52]

.

3. Modern Approaches to Simulation of Electromagnetic Processes in a Transformer

The modeling of electromagnetic processes serves to define heat losses in different parts of the power transformer, which can then be used in heat transfer simulation and forecasting thermal conditions. Additionally, the simulation of electromagnetic fields in transformers is employed to calculate refined parameters for simpler and more rapid models that take into account complex geometry and physical features. The accuracy of most simulating processes in the transformer depends on the numerical model of the magnetic field. Therefore, it is crucial to consider the main approaches used in modeling electromagnetic processes for various types of problems.

In [53], electromagnetic forces in coils and eddy current losses in the steel of a transformer are studied. For these purposes, the electromagnetic field is calculated using the $\mathbf{A}-\phi$ formulation of Maxwell’s equations (9) with Gaussian calibration (10). Gaussian calibration improves the convergence of the solution:

$${\nabla }^2\mathbf{A}+\sigma \mu \nabla \phi +\sigma \mu \frac{\partial \mathbf{A}}{\partial t}=\mu {\mathbf{J}}_{\mathbf{ext}}$$

(9)

$$\nabla ·\left.\left(\sigma \frac{\partial \mathbf{A}}{\partial t}+\sigma \nabla \phi \right.\right)=0$$

(10)

In these equations, $\mathbf{A}$ is the vector magnetic potential, defined by the expression $\mathbf{B}=\nabla \times \mathbf{A}$, $\phi$ is the scalar electric potential, expressed in terms of the electric field strength $\mathbf{E}=\nabla \phi$, $\sigma$ is the electrical conductivity of the medium, and ${\mathbf{J}}_{\mathbf{ext}}$ is the external current density. The use of potentials instead of magnetic induction for calculation makes it possible to increase the stability of solutions due to a smoother desired function. The number of model unknowns can be reduced by using the $\mathbf{A}$-formulation of the magnetic field

$$\sigma \mu \frac{\partial \mathbf{A}}{\partial t}+{\nabla }^2\mathbf{A}=\mu {\mathbf{J}}_{\mathbf{ext}}$$

(11)

The formulation of Equation (11) for the stationary case of a magnetic field, when $\sigma \mu \frac{\partial \mathbf{A}}{\partial t}=0$, is often used to analyze the mechanical characteristics in transformer windings during a short circuit in [54,55,56,57,58,59,60]. Accounting for the influence of forces is built in by introducing current in the transformer windings as an analytical expression, for example, $i_d\left(t\right)={\sqrt{2}I}_D[cos\alpha e^{-\frac{r_d}{L_d}t}-cos(\omega t+\alpha )]$ in [54], where $i_d\left(t\right)$ is a transient short-circuit current, $I_d$ is the RMS value of the steady current, $\alpha$ is the initial phase, and $r_d$ and $L_d$ are the resistance and inductance of the transformer coils, respectively. In [61,62], the dynamic formulation of the $\mathbf{A}$-formulation equation is used when $\sigma \mu \frac{\partial \mathbf{A}}{\partial t}\ne 0$ for similar studies. In problems where the magnetic field can be considered harmonic in time, it is advantageous from the point of view of computational resources to use the harmonic formulation of the magnetic field $\sigma \mu \frac{\partial A}{\partial t}=-j\omega \underline{\mathbf{A}}$, as in [63]. In this formulation, $\omega$ describes the angular frequency of the magnetic field, $\underline{\mathbf{A}}$ is the vector magnetic potential in a complex form, and $j=\sqrt{-1}$ is a complex number. The $\mathbf{T}-{\varphi }_m$ formulation of the magnetic field equations is also used, for example, to analyze the short circuit of separated transformer windings in [64]. This formulation is described by two differential equations:

$$\begin{array}{cc}\hfill \nabla \times \rho \nabla \times {\mathbf{T}}_{\Omega }-\nabla \rho \nabla ·{\mathbf{T}}_{\Omega }+\mu \frac{\partial {\mathbf{T}}_{\Omega }}{\partial t}-\mu \nabla \frac{\partial {\varphi }_m}{\partial t}& +\mu \frac{\partial {\mathbf{H}}_{\mathbf{s}}}{\partial t}=0,\hfill \\ \hfill \nabla ·\left({\mathbf{T}}_{\Omega }{-\varphi }_m\right)& +\nabla ·\mu {\mathbf{H}}_{\mathbf{s}}=0,\hfill \end{array}$$

(12)

where ${\mathbf{T}}_{\Omega }$ is the electric vector potential, ${\varphi }_m$ is the scalar magnetic potential, ${\mathbf{H}}_{\mathbf{s}}$ is the magnetic field intensity generated by the source current in an infinite space. The numerical model using the $\mathbf{T}-{\varphi }_m$ formulation matches well with models of the external circuit of a transformer based on equivalent circuits. In a number of problems, to reduce computational resources, models based on the magnetic scalar potential ${\varphi }_m$ (13) are used. In [65], this formulation is used to analyze the electromagnetic forces of a dry-type power transformer in three-dimensional space:

$$\mathbf{H}={\mathbf{H}}_j-\nabla {\varphi }_m,$$

(13)

where $\mathbf{H}$ is the magnetic field strength vector, and ${\mathbf{H}}_j$ is the source magnetic field strength vector.

3.1. Features of Calculation of the Transformer Coils

Maximum heat losses occur in transformer windings. The accuracy of the simulation will impact the forecasting of hot-temperature spots and equivalent parameters of transformers, which are crucial for calculating the convection of tank oil or building digital twins of power transformers. Therefore, the simulation of these parts receives special attention. The possibility of accounting for physical phenomena, such as local expansion, the interaction between coils with current, and others, depends on the method of mathematical representation of coils in models. However, many effects in coils are impractical to simulate directly using numerical simulation due to high computational resources. In such cases, combined models with numerical and analytical methods are employed. In this part of the work, the authors examine the main approaches to describing transformer coils and explore possibilities for model simplification.

The study of the electromagnetic behavior of the transformer is complicated by the need to model each part of the electromagnetic converter. A high voltage coil in some cases consists of about 3000 elementary thin conductors. In [66], it is proposed to replace the procedure for setting the current density in these turns by equivalent values of the coil load currents. This proposal is based on the assumption that the reaction field of high-voltage thin elementary conductors does not affect the distribution of the entire magnetic field. This assumption is confirmed by comparison with the experimental results of two different high-voltage coils. When using a current load, the procedure for creating a model mesh is greatly simplified.

Analytical methods are employed in transformer models to calculate inductance in [28,67]. The proposed models are verified by comparing their results with those of FE (finite element) models. The authors of [28] argue that their analytical model can be used to calculate transformers of complex shape, making it useful for design purposes. In [67], the Rogowski method is used to calculate inductance. In [31,32], analytical models for Litz-wire winding geometry are presented, which also suggest a means of evaluating losses. The main advantage of the analytical model is its ability to find an optimized geometry of transformer winding in a short amount of time. The results show that the use of the analytical method significantly expedites calculations but at the cost of accuracy in the results.

The accuracy of determining the losses in the coils of a transformer is important for evaluating its performance. Works [66,68,69] are devoted to this issue. In [68], the authors use a three-dimensional statement of the transformer model based on the finite element method to calculate eddy currents and losses in the transformer winding. The obtained results are compared with the test data and the analytical solution using simulation. The study [66] considers the use of an analytical approach for modeling additional losses in the transformer coil, which makes it possible to increase the accuracy of calculating the parameters of the transformer. In [70], a FE model is created to calculate losses in transformer coils due to Joule heating. From the results of this work, it can be concluded that the current density in the turns of the coils is unevenly distributed; therefore, the heat losses will also be distributed unevenly. The effect of uneven distribution of electrical losses can be taken into account only in spatially distributed models. In [69], the authors propose to introduce a correction factor $k_c$ to calculate the electrical losses in the transformer coils $P=R{I}^2{k}_c$. This coefficient makes it possible to take into account the influence of proximity effects and the skin effect, depending on the main power frequency and the size of the winding turns.

Another intricate aspect of simulating transformer windings is accounting for electrical insulation. The primary reason for considering electrical insulation in numerical models is to analyze insulation breakdown. Incorporating insulation elements into the transformer windings substantially complicates the mesh construction, increases the number of unknowns in the system, and requires the consideration of additional physical equations. However, this allows for a more comprehensive simulation of the physical process. For instance, in [71], the analysis of insulation breakdown is analyzed by estimating the field lines of force. In [72], the dielectric characteristics of oil–paper insulation under copper sulfide conditions are discussed based on changes in capacitance, insulation resistance, and partial discharges. An electric field analysis using FEM and partial discharge results is used to develop an empirical model for predicting the occurrence of partial discharge in paper–oil insulation in the presence of copper sulfide. Accounting for capacitive and inductive effects is possible using only the $\mathbf{A}-\phi$ formulation of the electromagnetic field or segregated models from different formulations. The effect of winding damage on the change in electrical parameters is considered in [73] using the finite-volume method. Also, accounting for electrical insulation is important in the analysis of the temperature conditions and humidity of the insulation of the transformer winding. In this type of problem, insulation is considered an obstacle with a certain temperature. This space region does not require significant mesh adjustments. Therefore, it is easier to take into account insulation in these problems than in the breakdown analysis. In [74,75,76,77,78], the influence of the skin effect, the effect of proximity and isolation on the calculations of losses in the transformer windings for the frequency range from 60 Hz to 2 MHz is analyzed. In these works, the features of mesh and model tuning are considered. Accounting for capacitive effects in numerical models is discussed in [77]. In [79], the influence of turn-to-turn insulation on magnetic and fluid-temperature fields has been presented. A numerical model implemented in ANSYS, which took into account the influence of each turn of the transformer windings, showed a difference in hot spot temperature ranging from 16.37 to 13.30 K compared to the same model without considering the influence. The authors argue that the influence introduces a ’wavy-shape’ distribution of magnetic flux, which does not affect the winding loss, average temperature of the winding, or oil flow distribution.

Calculating coil parameters such as capacitance and inductance is a critical stage in transformer simulation, as highlighted by [80]. The authors point out that a wide range of failures are related to transient oscillatory overvoltages and winding internal resonances, making it essential to model transformer windings in detail. In addition, the authors of [81] emphasize that modeling transient voltage distributions in the windings of electrical equipment is crucial for electrical insulation sizing and design optimization. They note that in some cases, such as lightning shock, the impact of the turn-to-turn capacitance, disk-to-disk capacitance, and capacitance of winding to ground cannot be ignored. A study by [82] discusses various approaches to calculating winding parameters. They divided the methods into three categories: Multi-conductor Transmission Line (MTL) [80,83], Lumped Circuit, and FEM. According to [82], the FEM approach is more universal since the same model can be used for both Frequency Response Analysis (FRA) and overvoltage analysis, and this approach has the lowest application complexity. However, it should be noted that the FEM method requires a 3D calculation model in most cases, and the simulation time is longer than in other approaches. In accordance with [80,81,82], the capacitance and inductance matrix can be derived from electrostatic and magnetostatic calculations, respectively. The capacitance matrix is calculated for cases where one conductor is under test voltage (e.g., 1 V), and the others are grounded. The calculation can be performed using the formulas provided:

$$C_{ij}=\frac{q_i}{V_j},$$

(14)

$$q_i=-\epsilon {\int }_{S_i}\nabla V·\mathbf{s},$$

(15)

where $q_i$ is the charge of conductor i, $S_i$ is the surface of the conductor, and $\epsilon$ is the permittivity of the medium. Similarly, the inductance matrix can be derived when one turn is under the test magnetic vector potential, and the others are set to zero. The inductance matrix can be calculated using the formulas provided:

$$\mathbf{\Lambda }={\mathbf{L}}^{-1},$$

(16)

$${\Lambda }_{ij}=\frac{I_i}{{\psi }_j},$$

(17)

$$I_i=-\frac{1}{\mu }\oint \nabla \times \mathbf{A}·d\mathbf{l},$$

(18)

where $I_i$ is the current of conductor i, ${\psi }_j$ is the magnetic flux between conductor j and the ground, and $\mu$ is the permeability of the medium.

3.2. Approaches to the Calculation of Transformer Steel

Losses in steel make up a significant part of all losses in a transformer. The accuracy of calculating this type of loss affects the accuracy of the transformer efficiency estimate. The modeling of losses in transformer steel consists in a complex mathematical description of the relationship between the magnetic flux density and magnetic field strength, which is called the material magnetic field (Equation (19)). The complexity of the mathematical description depends on the type of steel and how the hysteresis curve is described. In [84], it is proposed to divide the models representing the non-linear behavior of the magnetic field in steel into models that ignore the physics of the material and take into account the features of the material. In [85,86], reviews of the main methods for calculating losses in steel with their conditional classification are carried out. Based on the classifications in [85,86] and the own experience of the authors of this work in numerical modeling, a classification of models for calculating losses in steel is presented in Figure 6. In this section, we consider the main approaches to modeling the magnetic properties of steel parts of a transformer based on the presented classification.

The magnetic permeability of a medium is often used in describing the material magnetic field, shown in Equation (19)

$$\mathbf{B}={\mu }_0\mu \mathbf{H}.$$

(19)

Here, ${\mu }_0=4\pi ·{10}^{-7}$ is the absolute magnetic permeability of the vacuum, and $\mu$ is the relative magnetic permeability. The relative magnetic permeability $\mu$ is taken as a constant value in cases where the considered range of deviations of the magnetic flux density does not significantly affect the value of the magnetic field strength. In problems in which it is necessary to take into account the saturation of the magnetic field, the relative magnetic permeability is used as a function of flux density or strength. When using this approach, hysteresis losses are not taken into account. These losses can be taken into account in models using the semi-empirical formulas of Steinmetz (20) or Bertotti (21) [87] or using complex values of magnetic permeability [88]:

$$p_{fe}=C_{SE}{f}^{\alpha }{B}^{\beta }$$

(20)

$$p_{total}=C_{SE}{f}^{\alpha }{B}^{\beta }+k_c\sigma f^2{B}^2+{8\sqrt{\sigma GS{V}_0}f}^{1.5}{B}^{1.5}$$

(21)

The Steinmetz’s equation only calculates the hysteresis loss for an ideal sinusoidal field. The Bertotti equation makes it possible to take into account eddy current losses, and to take into account the non-linearity of the magnetic field, the Steinmetz component must be modernized, for example, as in (22):

$$p_{fe}{{\frac{1}{T}{\int }_0^T{C}_{SE}\left|\frac{dB}{dt}\right|}^{\alpha }\left|B\left(t\right)\right|}^{\beta -\alpha }dt$$

(22)

A detailed description of the development of the Steinmetz equation is displayed in many modern articles, for example, in [85,86]. In [89], three variants of models for calculating losses in steel based on the Bertotti theory are considered. The studies are carried out for sinusoidal and non-sinusoidal magnetic field forms using FE models. The paper provides fitted coefficients for Bertotti loss models based on measured data. The transformer core loss models designed for harmonic magnetic fields use an effective non-linear magnetic saturation curve $B_{eff}\left(H\right)={\mu }_{reff}H$. The values of ${\mu }_{reff}$ in this curve are overestimated compared to the real values of $\mu$ to increase the convergence of the problem solution. The effective values of ${\mu }_{reff}$ are determined from the simple energy method [90]. The obtained results of the magnetic induction $B_{eff}$ using this method should be recalculated to real values of B. The calculation of magnetic losses in a transformer using this method is implemented in [91].

The work [60] considers modeling the behavior of the magnetic field in three-phase transformers, the cores of which are made of isotropic and anisotropic steel. The work takes into account the saturation of steel with a magnetic field using the dependence $B\left(H\right)$. The paper also takes into account hysteresis losses depending on the value of the magnetic flux density modulus. Anisotropy is taken into account using a family of curves that are built for different rolling angles. Data on iron losses in the saturation region at high flux densities are not available in the catalog, so they are obtained using linear interpolation. The Piecewise Linear Saturation Model [24,92,93] is commonly employed to calculate steel losses and account for the saturation phenomenon. This model only requires values for the saturation curve, which are specified for each steel grade. By using this model, transformer might be simulated in real time [24]. The distribution of average losses over time in the active zones of a conventional magnetic circuit, a magnetic circuit with an oblique connection, and a magnetic circuit with a butt connection using anisotropic steel plates is studied in [44].

References [94,95,96] present an approach to modeling the magnetization of a transformer core under constant bias conditions using the inverse Jiles–Atherton (J-A) model. The implementation of this model consists in determining five parameters that describe the behavior of the magnetic hysteresis in steel. In the J-A model, the magnetization M is represented as the sum of the reversible $M_{rev}$ and the irreversible component $M_{irrev}$

$$M=M_{rev}+M_{irr}.$$

(23)

The reversible component (24) can be expressed in terms of the hysteresis-free magnetization $M_{an}$ and the constant coefficient c.

$$M_{rev}={c(M}_{an}+M_{irrev}).$$

(24)

The hysteresis-free magnetization $M_{an}$ is described by function (25) using the effective component of the magnetic field strength $H_e=H+\alpha M$:

$$M_{an}=M_s\left[coth\frac{H_e}{a}-\frac{a}{H_e}\right].$$

(25)

In (25), $a=k_{b}T/{\mu }_{0}m$ is the parameter of the Jiles–Atherton model, which can be interpreted as a measure of the connection between the neighboring magnetic domain, and m is the magnetic moment. Then, the magnetization is represented as a function of two derivatives based on the effective magnetic induction $B_e={\mu }_0{H}_e$:

$$\frac{dM}{dB}=\frac{\left(1-c\right)\left(\frac{d{M}_{irr}}{d{B}_e}\right)+\frac{c}{{\mu }_0}\frac{d{M}_{an}}{d{H}_e}}{1+{\mu }_0\left(1-c\right)\left(1-\alpha \right)\left(\frac{d{M}_{irr}}{d{B}_e}\right)+c\left(1-\alpha \right)\left(\frac{d{M}_{an}}{d{H}_e}\right)},$$

(26)

$$\frac{d{M}_{irr}}{d{B}_e}=\frac{M_{an}{M}_{irr}}{{\mu }_{0}k\delta }.$$

(27)

$\alpha$, a, k, c, $M_s$ are the desired five parameters that describe the hysteresis curve. The coefficient $\delta$ describes the sign of the derivative of the magnetic induction, for $\frac{dB}{dt}>0,\delta =1$, and for $\frac{dB}{dt}<0,\delta =-1$. The equations mentioned above can be reduced to the final equation

$$\begin{array}{c}\hfill \frac{dM}{dH}=c\frac{M_s}{a}[1{-(coth\frac{H_e}{a})}^2+{\left(\frac{a}{H_e}\right)}^2](1-c)\frac{(M_{an}-M)}{k\delta (1-c)-\alpha (M_{an}-M)}/(1-\alpha c).\end{array}$$

(28)

The Jiles–Atherton model and its variations is one of the most accurate numerical models for reproducing the behavior of a magnetic field in ferromagnetic parts. The disadvantage of this model is the complexity of its settings, as it is necessary to select parameters that are based on experiments and significantly increase the computational requirements for solving the problem. Many works are devoted to the creation of machine learning algorithms for finding the coefficients of this model, for example, in [84,94], and in [97], genetic optimization algorithms are used; in [98] the method of stochastic optimization “imitation of annealing” is used; and in [99], hybrid optimization methods are used.

Another group of models commonly used for capturing hysteresis behavior is the history-dependent hysteresis models, which include the Preisach, Play, and Stop models [100,101,102]. These models are effective in characterizing the behavior of laminated steels in iron cores as well as permanent magnets. Among them, the Preisach model is widely adopted, while the play model provides a simpler alternative. In a comprehensive analysis conducted by the authors in [100], it was concluded that the stop model might be more suitable, particularly when finite element method (FEM) computations solve Ampere’s law using the inverse form. This preference stems from the fact that the stop model represents hysteresis in an opposite manner to H(B), in contrast to the common form B(H). Furthermore, the stop model offers simplicity in directly linking to iron losses, as the energy is defined as the integral of H(B) instead of the coenergy B(H).

The developed dynamic model for the accounting of hysteresis is based on the fuzzy Gustafson–Kessel (GK) approach. The use of normalized collected data from measured dynamic cycles of a C-core transformer, made of 0.33 mm thick cold-rolled SiFe plates, is presented in [103]. This paper discloses an algorithm for finding and compiling a dynamic hysteresis curve based on measurements.

In [104], the dependencies of hysteresis losses and eddy currents are given as a function of the magnetic flux density values. Characteristics are given for isotropic, anisotropic and bidirectional steels. These curves are derived from expression (29)

$$P=\frac{f}{\rho }{\int }_0^T\left(\frac{d{B}_x\left(t\right)}{dt}{H}_x\left(t\right)+\frac{d{B}_y\left(t\right)}{dt}{H}_y\left(t\right)\right)dt.$$

(29)

where $\rho$ is the material density; $B_x\left(t\right)$, $B_y\left(t\right)$ are the time dependencies of the magnetic flux density in the x and y directions; $H_x\left(t\right)$, $H_y\left(t\right)$ are the time dependencies of the magnetic field strength in the x and y directions; f is the field frequency; and T is the network period. Loss versus induction results for isotropic steel are shown in Figure 7.

Machining of steel products affects the magnetic properties of the transformer. Modern models allow taking into account the effect of machining. In [105,106], approaches to modeling these effects are considered.

It is expedient to finish the discussion of the paragraph on the calculation of losses in steel by comparing the considered mathematical models. Comparison of model options for calculating losses in steel, such as the ability to take into account the complex shape of the magnetic field (non-sinusoidal), the need for preliminary calculations, and the accuracy of describing the hysteresis phenomenon can be summarized in

Table 2. Comparison of loss models.

Loss Model	Complex Field Shape	Prior Knowledge of Material Properties	Accuracy
The Steinmetz equation	−	Low	Low
Modified Steinmetz equation	+	Low	Low-Medium
Improved generalized Steinmetz equation	+	Low	Low-Medium
Bertotti losses	+	Medium	Medium
Dynamic hysteresis model	+	High	Good
Surface loss model	+	High	Good
A model based on magnetodynamic viscosity	+	High	Good
Friction as a hysteresis model	+	High	Good
Hysteresis model based on energy	+	High	Good
Separation of losses after magnetization processes	+	Depends on the model	Depends on the model

. This table was developed on the basis of works [85,86].

3.3. Selection of Equivalent Parameters Based on Numerical Simulation

Transformer equivalent circuits are widely used in control and diagnostic systems, in many design stages and in tasks where a simple mathematical apparatus and fast calculation are required. The accuracy of this method in some problems is not acceptable. In [73], the parameters of the equivalent circuit are refined based on the results of the finite element model for diagnosing axial and radial deformation. In the study [37], the frequency-dependent parameters for partial discharge diagnostics are improved with the FE model. The FEM is used to assess the influence of the transformer core and transformer windings equivalent circuits parameters on the frequency analysis, for example, in [77,78,107,108]. In [109], a genetic algorithm is used to determine the parameters of a laminated core of an equivalent tiling circuit based on frequency analysis. A similar study is carried out in [110] to improve the process of simulation of the winding of a power transformer with the bacteria accumulation algorithm and frequency analysis.

4. Auxiliary Part Simulation

In numerical modeling, by using the finite element method and the finite volume method, there is always the problem of reducing the detail of the object of study since each structural element complicates the meshing process due to the increased number of calculated nodes. Also, it increases the number of degrees of freedom, which ultimately increases the calculation time. Therefore, it is crucial to leave only those structural elements that can affect the final result. For that reason, frequently, only the main elements are considered, such as the magnetic circuit and windings. However, it is essential to understand whether additional elements can change the distribution of values under study. Under such consideration, additional metal elements might include the yoke clamps, radiator plates, and an oil tank, which are schematically presented in Figure 8.

In papers [50,111,112,113], issues of modeling auxiliary parts are considered. The authors of the work [50] argue that the additional internal metal parts of the transformer, in which heating occurs due to induced currents, can affect the accuracy of calculating the thermal field. To obtain the distribution of heat release power due to electromagnetic influence, the authors carry out an additional calculation of the electromagnetic problem. In accordance with the results from [50], the distribution of the field is uneven, and there are areas of local increase in magnetic flux density, in which there is increased heat release due to electromagnetic losses. For convenience, the relative values of heat release are summarized in

Table 3. Thermal output power distribution.

Location	Percentage of Total Losses, %
Oil tank	1.22
Internal metal structural parts	2.27
High voltage winding	54.61
Low voltage winding	42.39

. The integral values clearly show that the additional parts are not the main source of heat generation, but even such values can affect the result. The verification of this approach shows that the difference in the steady-state temperature of the installation with the experimental data is no more than 3 °C, which indicates a fairly high convergence of the results.

The paper [112] also considers the auxiliary parts but for the auto-transformer case. The main purpose of this work is to test the approach with a simplified calculation of the hydrodynamic problem, which consists in using the equivalent heat transfer coefficient, which is calculated for each of the key surfaces. But in this case, the results for the auxiliary parts are of greater interest. As in [50], the authors obtain an uneven distribution of the heat release density in the parts of the transformer, those areas that are closer to the windings are under greater electromagnetic influence. In the result, the main losses are concentrated in the transformer tank. As a consequence, the uneven distribution of temperature in the structural parts of the transformer is obtained. In the conclusion, the authors state that the comparison of measured values and calculated results show a relative error not exceeding 9%.

The study [111] considers the use of magnetic shunts and shields from induced current in order to reduce electromagnetic losses and the highest temperature. The results of the impact of various loss reduction measures are summarized in

Table 4. Comparison table of various loss reduction measures [111].

Measure Applied	Losses, kW	Tank Temperature, °C	End-Frame Temperature, °C	Flitch Plate Temperature, °C
Without screens	13.18	85.18	87.17	160.98
With wall shunts	6.07	71.91	82.21	156.23
With copper screens	11.51	74.88	89.43	160.26
Combined use of wall shunt and copper screen	5.71	72.01	75.88	155.67
Combined use of wall shunt and shunt in magnetic circuit	6.09	71.11	75.81	155.33

.

As a result, the authors claim that magnetic shunts reduce losses and the highest temperature to a greater extent than aluminum or copper shields. When using shunts in the magnetic core, the highest temperature in the end frames is significantly reduced. This paper shows how modeling additional parts of the transformer can significantly affect the temperature field in the transformer. In [63], with a combination of the FEM and analytical method, the same efficiency for vertical and horizontal placements of the transformer shunts is shown, but the former weigh 25% less.

5. Conjugate Heat Transfer Simulation Review

Understanding the thermal dynamics of a transformer is a key factor in ensuring the reliability of its operation [15,114,115]. Considering that almost all properties of materials (thermal conductivity, heat capacity, density, electrical conductivity, etc.) depend on temperature, the dynamics of changes will be non-linear. In addition, the analysis of condition is complicated by the fact that the transformer load is not constant during the day and throughout the year. Mathematical models for transformer thermal simulation can be broadly classified into two categories: those based on circuit theory [39,40,41,42], and those that solve field problems using numerical methods [116,117,118,119,120,121,122,123]. The first group has the advantage of allowing for real-time calculations, as they do not require significant computational resources. However, these models rely on empirical or analytical dependencies with multiple assumptions to approximate ongoing physical processes, which can limit their accuracy [124]. They also fail to provide a comprehensive understanding of the internal workings of the transformer, such as fluid flows and temperature distribution. The second group, on the other hand, offers a more complete picture of the ongoing processes in the transformer, as the output data can include temperature field distributions, velocity fields, and the like. However, these models require significantly more computational resources and time to set up and run, making them less practical for real-time applications.

Modeling using equivalent circuits (dynamic models) involves forming an ordinary differential equation (ODE) that includes terms responsible for heat sources, such as heat generation in windings and magnetic circuits, as well as passive elements like heat capacity and thermal resistance. This approach enables the calculation of temperature changes over time in selected parts. The simplest example of such an equation is relationship (30), which describes the heat exchange between oil and air for a core-type two-winding oil-immersed transformer (Figure 9) [40]. It should be acknowledged that the described dependence may undergo significant changes in the case of special design configurations or under extreme conditions, such as in the event of a fault.

$$q_{fe}+q_{cu}=C_{coil}·\frac{d{\theta }_{oil}}{dt}+\frac{1}{R_{oilR}}·{\left[{\theta }_{oil}-{\theta }_{amb}\right]}^{\frac{1}{n}},$$

(30)

where $q_{fe}$ is the heat generated by losses in steel, $q_{cu}$ is the heat generated by losses in copper, $C_{coil}$ is the heat capacity of oil, ${\theta }_{oil}$ is the oil temperature, $R_{oilR}$ is the thermal resistance of oil at nominal conditions, and ${\theta }_{amb}$ is the ambient temperature. The technique described in the IEC 60076-7 [125] standard has a similar implementation, but according to [126], dynamic models require less input data. It is worth noting that such models can be more complex, including more parameters [126]. For example, in work [127], the authors present a gray-box model of an oil-filled transformer based on a thermoelectric analogy. The resulting model can be used for the online monitoring of hot-spot temperature (HST) and oil temperature. From a future perspective, this information can help the operator load the equipment more optimally. Work [128] shows a comparison of the modified dynamic model with the methods from the standards (IEC 60076-7, IEEE C57.91-2011 [129]). The authors claim that the described approach shows better results than those described in the standards. A further development of the dynamic model can be found in work [130], where the thermal resistance of paper insulation is additionally considered, taking into account changes due to humidity. Through comparison with models without moisture, the authors prove that their model shows more accurate results. Another work [131] presents an approach to take into account weather influence. A top-oil thermal model based on a thermal–electrical scheme analogy is developed to define the hot-spot temperatures of two different power transformers under various conditions.

The next step in constructing thermal equivalent circuits is the detailed thermal equivalent circuits or Thermal Network Model. This approach combines the advantages of simple equivalent circuits and models based on the finite element or volume method. The resulting model does not require significant computational resources, and its parameters are spatially dependent. Spatial discretization allows for the consideration of the vertical temperature gradient, which is a weakness of simple equivalent circuits [132]. However, this approach differs from FVM and FEM in that it uses analytical and empirical dependencies for regions with liquid, and therefore, it does not account for the uneven distribution of fluid velocity along the hydrodynamic channel (flow profile). As a result, some phenomena associated with heat transfer are not taken into account [133]. Examples of the implementation of such an approach can be found in various works [132,134,135], such as the modeling of hydraulic circuits for oil-immersed power transformers [134], which demonstrate a deviation of less than 1 °C in hot spot temperature and 1.4 °C in average temperature between computational fluid dynamics (CFD) and hydraulic circuit modeling results. In work [135], the authors propose a dynamic Thermal Network Model of an oil-immersed on-board traction transformer. As a result, the described model helps to save more than 99% in time compared to the CFD model. Additionally, the average relative error is lower than 2.57%, verifying the accuracy of this approach.

Numerical calculation of the temperature field using the partial differential equation for energy conservation enables a more accurate prediction of the thermal conditions of the transformer. In most engineering problems, this equation is written as:

$$\begin{array}{c}\hfill \rho C_p\left[\frac{\partial T}{\partial t}-\mathbf{u}·\nabla T\right]=-\nabla (\lambda \nabla T)+Q,\end{array}$$

(31)

where $\rho$ is the fluid density, $\lambda$ is the thermal conductivity, Q is the volumetric heat source, $C_p$ is the specific heat capacity, T is the temperature, and $\mathbf{u}$ is velocity of a flow. The convective term $\mathbf{u}·\nabla T$ is accounted for only in fluids parts, for instance oil, and is equal zero for all solid parts of the power transformer. In order to calculate the velocity $\mathbf{u}$, it might be performed by the Navier–Stokes equation, being consisted of conservation of momentum (32) and mass Equation (33):

$$\begin{array}{c}\hfill \frac{\partial \rho \mathbf{u}}{\partial t}-\left(\rho \mathbf{u}·\nabla \right)\mathbf{u}=-\nabla p+\mu {\nabla }^2\mathbf{u}+\rho \mathbf{g},\end{array}$$

(32)

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}+\nabla ·\mathbf{u}& =0.\hfill \end{array}$$

(33)

Here, p is the pressure field, $\mu$ is the dynamic viscosity, and $\mathbf{g}$ is the gravity acceleration. The motion of liquid flows is caused by changes in mass density $\rho$ and the non-uniform distribution of buoyancy force $\rho \mathbf{g}$. The above-mentioned Equations (32) to (33) are written for compressible flows, but a solution of this task requires significant computational resources. The mass density of transformer oil varies slightly over the entire operating temperature range. Therefore, the equations can be simplified using Boussinesq’s approximation ${\rho }_0\beta \left(T-T_0\right)\mathbf{g}$ under the condition $\Delta \rho =\rho \left(T_0\right)-\rho \left(T\right)\ll \rho \left(T_0\right)$. Here, the thermal expansion coefficient equals $\beta =\frac{1}{{\rho }_0}\frac{\partial \rho }{\partial T}$ and can often be set as a constant value for each fluid. It should be noted that ${\rho }_0$ is the reference value of mass density at the corresponding temperature $T_0$. As a result, Equations (32) to (33) can be replaced with:

$$\begin{array}{c}\hfill \frac{\partial \mathbf{u}}{\partial t}-\left(\mathbf{u}·\nabla \right)\mathbf{u}=-\nabla p_k+\nu {\nabla }^2\mathbf{u}+\beta \left(T-T_0\right)\mathbf{g},\end{array}$$

(34)

$$\begin{array}{cc}\hfill \nabla ·\mathbf{u}& =0,\hfill \end{array}$$

(35)

where $\nu =\mu /\rho$ is the kinematic viscosity, and $p_k=p/\rho$ depicts the kinematic pressure. Solving this type of equation requires significantly fewer computational resources because it does not require taking into account the non-linearity of the density. In [116], a study is conducted on a 30 MVA power transformer radiator operating in oil using Boussinesq’s approximation to calculate convection. A novel approach that takes into account secondary flows in “Code_Saturne” is used to decouple thermal and fluid-dynamics issues. In [117,118], a coupled thermal and hydrodynamic calculation based on Boussinesq’s approximation is used to determine the thermal conditions of the transformer, aiming to identify weaknesses in the structure and optimize it further.

Numerical simulation of the temperature field of general parts of an oil-immersed power transformer, including conjugate heat transfer, requires significant computational resources and high-quality engineering. One way to reduce the complexity of the numerical model is to consider heat transfer only in the solid parts of the power transformer, without calculating hydrodynamic processes using convective correlation based on empirical, semi-empirical, or calculated relationships of heat fluxes on the boundaries of solid parts. In this case, the heat removal by transformer oil from solids is considered by the boundary condition of convective heat flux $q_{conv}$:

$$\begin{array}{c}\hfill q_{conv}=h·(T_{oil}-T_{wall}),\end{array}$$

(36)

where $T_{oil}$ and $T_{wall}$ are the temperatures of the oil and wall of the solid part, and h is the heat transfer coefficient which is dependent on the geometry parameters, such as shape and its measurements, and Nusselt ($Nu$) and Prandtl ($Pr$) dimensionless numbers. The model in [131] can be applied to calculate the heat transfer coefficient of a radiator in a large power transformer using the Nusselt number. It describes the relation between the total convective flux and the conductive flux, and is a complex dependency of Rayleigh ($Ra$) and Prandtl number for natural convection. For instance, in [51], the Nusselt number for vertical (37) and horizontal (38) natural convection can be defined as

$$\begin{array}{c}\hfill N{u}_y=0.68+\frac{0.67R{a}^{1/4}}{{\left[1+{\left(0.492/Pr\right)}^{9/16}\right]}^{4/9}},\end{array}$$

(37)

$$\begin{array}{c}\hfill N{u}_x=0.59R{a}^{1/4},\end{array}$$

(38)

The above-mentioned expressions are valid for laminar fluid flow around a long plate. When calculating a transformer, the Nusselt number is often expressed in another way, for example, in [136,137] as (39), and in [51,138,139] as (40):

$$\begin{array}{c}\hfill Nu=0.85\left(0.74R{e}^{0.2}{\left(Pr·Gr\right)}^{0.1}P{r}^{0.2}\right),\end{array}$$

(39)

$$\begin{array}{c}\hfill Nu=0.15{\left(Gr·Pr\right)}^{1/3},\end{array}$$

(40)

The Rayleigh number $Ra=Gr·Pr$ allows for evaluating the relationship between buoyancy and viscosity within a fluid, while the Prandtl number $Pr=C_p\mu /\lambda$ describes the relationship between convective and diffusion terms of heat transfer. In (39), the Reynolds number $Re=UL/\nu$ is used, which defines the ratio of inertial force to viscous forces, where U and L are the characteristic velocity of flow and length. The Grashof number $Gr$ is introduced to calculate Rayleigh number $Ra$, which represents the ratio of buoyancy body force to viscosity force. Unlike the Prandtl number, which is uniquely determined, the Grashof number (41) is described in various works by different characteristic parameters, which can significantly affect the value of the heat transfer coefficient h in (36):

$$\begin{array}{cc}\hfill Gr& =\frac{g\beta \Delta T{L}^3}{{\nu }^2}.\hfill \end{array}$$

(41)

In [136,137], the Grashof number is calculated using $L=D_h=4A_c/P_c$, which represents the hydraulic diameter, namely, four cross-sectional areas $A_c$ divided by the wetting perimeter $P_c$. The characteristic temperature difference $\Delta T$ is defined as the surface temperature drop from the oil to the cooler wall. In [140], $\Delta T=T_{avg}-0.5(T_b+T_t)$ is expressed through the average value of the winding temperature $T_{avg}$ and the temperature in the upper $T_t$ and lower parts $T_b$ of the tank. In [51,138,139], the characteristic temperature $\Delta T$ is defined as the temperature difference between the winding and the cooled oil, and the characteristic length L is expressed by the height of the transformer coils. In [138], the characteristic size L is taken to be the wire length. In [138], the heat transfer coefficient $h=Nu·\lambda /L$ depends on the Nusselt number and the characteristic length L, so it is important in future works to consider the influence of the value of the calculated coefficient h on the selected characteristic parameters when calculating it.

It should be noted that the approach using the heat transfer coefficient may yield slightly distorted results, as heat transfer is non-uniformly distributed on the walls of the transformer winding as shown in [141]. Additionally, the heat transfer coefficient varies greatly depending on temperature and other system parameters. Therefore, studies are being conducted to improve and correlate heat transfer coefficients for different geometries and operating conditions. In [42,142,143], the authors develop simple models to predict high-temperature spots, but the heat transfer coefficient’s value significantly depends on the system’s geometry, hydrodynamics, and heat parameters, which are not considered in those studies. In [144], the authors propose a method to predict high-temperature spots using a self-written MATLAB code based on the finite element method for the heat conduction equation. They evaluate the heat transfer coefficient using semi-analytical expressions, similar numbers, and experimental data. The suggested approach allows for the prediction of high-temperature spots that correspond reasonably well with simulation and experimental data. The presented results are demonstrated on a 35 MVA transformer with layer windings. In [145], the author conducts a study on the dependence of heat transfer coefficient for five different types of transformer windings. The study finds that each winding type needs to be considered separately for both natural and forced convection. The author concludes that hydraulic resistance is the most significant factor in describing universal heat transfer in transformer windings. This conclusion highlights the importance of developing numerical conjugate heat transfer models or conducting experimental tests to define the hydrodynamic parameters of power transformer systems. In [119], a numerical model is developed based on the finite element method, considering only heat transfer equations. The model is correlated using measurements from 108 thermocouples and is asserted to have an accuracy of $\pm 1$ °C by the authors. However, the model only weakly describes the non-linear physical properties of transformer oil, which is a limitation, as different oils are used in power transformers and exhibit diverse heat transfer behavior as demonstrated in [146].

The problem of coupled heat transfer simulation is quite difficult to solve because it requires solving the thermal and hydrodynamic problems together. This becomes especially difficult in objects with complex geometry, such as an oil-filled transformer. The gaps between the windings and the magnetic core form narrow hydrodynamic channels that require a high-resolution mesh. As a result, researchers often use a two-dimensional axisymmetric approximation of the problem to simplify the calculation process. In a three-phase oil-filled transformer, only one of the electrical phases is considered, which significantly reduces the calculation time and the complexity of setting up the finite element mesh. This approximation is based on the assumption that the load in all phases is symmetrical and that all hydrodynamic and thermal processes are also symmetrical in nature. An example of such an approach is work [120]. In it, the authors explore both two-dimensional and three-dimensional formulation of the problem of coupled heat transfer in a disk-type winding. The authors compare the results of the 2D and 3D models with benchmarks and find that both models show good convergence. However, they also note that three-dimensional effects can significantly influence the results, such as the location and value of the hottest point. These conclusions are obtained in [121], where the authors study the heating of a disk-type low-voltage winding in 2D and 3D. In the three-dimensional formulation, the details of the winding design are taken into account, which have a significant impact on the result. The authors propose using corrective dependencies for two-dimensional results to simplify the calculation process. This approach shows good convergence with the three-dimensional results. However, it is important to understand under what conditions it is permissible to use simplified models without adjustments. The work [45] demonstrates the calculation of the hydrodynamic and temperature fields through numerical simulation, with both two-dimensional and three-dimensional models. In the course of the study, the authors vary the fluid velocity, thereby changing the Reynolds number. One result is that at certain Reynolds numbers (greater than 1000), backflows occur in certain areas of the transformer. This phenomenon is not taken into account by the two-dimensional formulation of the problem, so the temperature of the transformer windings will be calculated incorrectly. In the considered case, the winding temperature is 20 degrees higher in the two-dimensional formulation than in the experiment and three-dimensional calculation. Also, with the entrainment of the flow velocity, the overlap of some channels is observed due to the non-uniformity of the liquid flow and the occurrence of reverse flows, shown in Figure 10. In this case, the hydrodynamic similarity criterion serves as a convenient indicator of the possibility of using simpler models. In subsequent works, the authors continue their research focused on describing processes in a transformer using similarity criteria. In the second part of the study [140], the authors concentrate on examining free cooling modes. A notable result is that the authors propose using the Richardson number to describe the phenomena occurring in an oil-filled transformer. The work [140] demonstrates that this similarity criterion correlates strongly with the degree of flow uniformity, as well as the highest temperature of the windings.

It is worth noting that although three-dimensional models are more challenging to implement, they are often employed in transformer design [122,123]. These models can be utilized to evaluate the design or optimal placement of measuring sensors. For instance, in study [122], the authors propose the use of spatially distributed temperature sensors, and the three-dimensional modeling approach helps to analyze physical processes in greater detail. In another study [123], numerical modeling is used to determine the temperature of the hottest point, which can aid in transformer design. However, if more accurate results are required, a more comprehensive model of the physical processes responsible for heat transfer may be necessary. In addition to associated heat exchange, such a model would also include the modeling of electromagnetic processes, enabling a more accurate description of the heat sources distribution within the transformer. In [51], the authors propose an approach to transformer simulation by compiling three numerical models: an electromagnetic field distribution model, a fluid dynamics model between transformer windings, and an oil tank hydrodynamics model. The authors use an interesting approach; in addition to direct numerical calculations, they use a reduced-order model. This makes it possible to simplify the task in those areas where high accuracy of results is not required. The authors note that as a result of the calculation, the discrepancy with the experimental data does not exceed 0.1%. The models described possess high accuracy, but their computationally intensive nature poses a challenge when creating digital twins that must operate in real time. A feasible solution in this scenario is to employ reduced-order models, which can be derived from the results of high-precision field calculations through regression analysis or projection onto a lower-dimensional space. The work [147] can be considered one of the examples of such an approach. With numerical models, the authors calculate lots of oil-cooling alternatives in radiators under various initial and boundary conditions, after which, based on machine learning algorithms, they obtain dependencies for calculating the heat transfer coefficient on transformer radiators, which is given in (42):

$$h_{air}=1.499\times {\left(T_{avo}-T_{amb}\right)}^{0.4158},$$

(42)

where $h_{air}$ is the equivalent heat transfer coefficient, $T_{avo}$ is the average temperature, and $T_{amb}$ is the ambient temperature. This study allows further replacement of complex numerical calculations of radiators with a reduced-order model, which is demonstrated in [51]. This approach can be applied to other areas in transformer simulation, for example, in the calculation of electromagnetic losses. Another research team [9] utilizes the results of coupled heat transfer in a transformer to develop a prediction model for hot spot temperature (HST). By combining calculation results, data from oil temperature sensors, environmental factors, and load values, the authors train a support vector machine (SVM) prediction model, resulting in an average error rate of 0.587%. While these findings are valuable for modeling transformers and predicting specific lumped parameters, creating a digital twin capable of simulating a transformer’s state requires accounting for the distribution of physical quantities within the volume as demonstrated in [8]. This paper presents the development of virtual sensors based on the results of coupled heat transfer calculations using numerical methods. Similar distributed results can be attained through dimensionality reduction techniques, which significantly expedite the calculation process. For instance, in [11], the authors mention that employing a model based on the finite element method can be challenging to meet calculation requirements. Therefore, they suggest using the proper orthogonal decomposition (POD) method to construct a reduced-order model of the transient fluid–solid coupling temperature field. This approach results in a model that requires 192 times less calculation time while maintaining sufficient accuracy of the results. A similar approach is presented in [148], where the authors utilize a numerical model of a transformer that accounts for electromagnetic, thermal, and hydrodynamic processes as the foundation for training their algorithm. The algorithm is based on the dynamic mode decomposition (DMD) procedure, which can reduce the number of degrees of freedom in numerical simulation results. This allows for efficient training of the model, enabling subsequent forecasting of any parameters depending on the initial conditions. A similar mathematical model for state forecasting is presented in (43):

$$x_f={\Phi }_d{\Lambda }_d^{f-1}{B}_d,$$

(43)

where $x_f$ is the future state of the system, ${\Phi }_d$ is the new state matrix including dominant modes, ${\Lambda }_d$ is the new eigenvalue matrix advising dominant modes, and $B_d$ is the amplitudes corresponding to dominant modes. The authors compare their low-order model with the numerical model results, and the forecasted values of the highest winding temperature have good convergence (the difference in absolute values is no more than 2 K).

The temperature distribution simulation is an important step in monitoring the condition of a transformer. A literature review shows that there are different approaches to simulation. When using the thermoelectric analogy, it is possible to obtain the highest temperature of the windings without significantly studying the model, and the output will be a model that allows real-time monitoring of the state of the installation. However, this approach necessitates the refinement of empirical dependencies, and it does not yield high accuracy. Therefore, we can say with confidence that the simulation of an oil-filled transformer is inextricably linked with the calculation of the hydrodynamic processes in it. This fact greatly complicates the stage of creating a model but allows a more detailed analysis of thermal processes. It is also worth noting that the use of machine learning algorithms opens up new prospects in this area since the combined use with numerical models makes it possible to obtain simple dependencies at the output as in the thermoelectric analogy.

6. Discussion

Analytical approaches are highly valuable due to their ease of integration into software applications and control systems, making them extremely useful for diagnostic purposes and the preliminary design stage. However, the development of a universal procedure for creating analytical solutions that accurately describe complex phenomena poses significant challenges. Analytical expressions often rely on simplifying assumptions and idealized models to derive equations. Since power transformers are intricate devices that encompass various non-linearities, such as magnetic saturation, hysteresis, and eddy current losses, accurately capturing these effects within analytical expressions can be arduous and may necessitate additional approximations. Another limitation arises from the fact that analytical expressions are typically derived for specific idealized scenarios or simplified models. They may not fully account for the complexities and variations encountered in real-world power transformers. Consequently, their applicability can be restricted to certain design configurations or operating conditions. These expressions often assume specific design constraints or limitations to simplify calculations. However, these constraints may not align with the unique requirements or restrictions of a particular transformer design. As a result, the accuracy and validity of the analytical results can be compromised in such cases. Power transformers exhibit non-linear behavior, particularly under load variations or during fault conditions. Analytical expressions may struggle to accurately capture these non-linear effects, leading to less precise predictions of transformer performance under dynamic or abnormal operating conditions. While analytical methods can expedite calculations, they may not provide the desired level of accuracy for certain applications. Therefore, analytical approaches are often used to simplify calculations in other mathematical models by replacing complex components with simpler alternatives. The method of representing a transformer using equivalent circuits is widely used in practice for real-time simulations. This approach can significantly reduce computation time and requirements for calculation. However, it has disadvantages, including low accuracy, especially for complex non-traditional transformer designs, poor accounting for non-linearity, and no possibility of analyzing local system parameters. Another approach is the implementation of numerical methods, which can provide approximate solutions with high accuracy. However, the number of equations to be solved increases with the discretization level. Finite element and finite volume methods are two well-known numerical methods used to solve partial differential equations describing physical processes in oil-immersed transformers. The finite element method is more sensitive to mesh quality and type, while the finite volume method is faster in solving problems involving hydrodynamics. However, the finite element method offers more flexibility in using various mesh types and provides a continuous solution. This difference can be crucial in certain problem classes, such as electrodynamics. A comparison of the considered mathematical methods from different aspects of their application is provided in Table 1 of this review.

The accurate modeling of electromagnetic processes is crucial in predicting heat losses in power transformers, which can then be used for heat transfer simulation and forecasting thermal conditions. Different approaches are used in modeling electromagnetic processes for various types of problems. Using potential formulations instead of magnetic induction improves solution stability and optimizes the calculation process. For instance, the harmonic $\mathbf{A}$-formulation of the magnetic field reduces the number of model unknowns and reduces requirements for computational resources in problems where the magnetic field is harmonic in time. Also, in cases where it is necessary to take into account short-circuit conditions with combinations of external circuits, it might be useful to use the $\mathbf{T}-{\phi }_m$ formulation of the magnetic field equations. Simulating transformer windings is a difficult task. In this case, it is helpful to use combined models with numerical and analytical methods to streamline calculations. Additionally, it is useful to use some assumptions, such as using equivalent values of the coil load currents instead of setting the current density in thin conductors, or not considering each turn of the transformer winding. According to considered studies, the finite element method is the most universal and has the lowest application complexity in the case of calculating windings parameters. The modeling of losses in transformer steel requires a sophisticated mathematical description of the relationship between magnetic flux density and magnetic field strength, which is commonly referred to as the material magnetic field equation. In certain cases, the relative magnetic permeability is assumed to be a constant value, as the variations in magnetic flux density do not significantly impact the magnetic field strength. However, when the saturation of the magnetic field needs to be considered, the relative magnetic permeability is treated as a function of flux density or strength. Various approaches can be employed for this purpose, including the Piecewise Linear Saturation Model, the Jiles–Atherton model and its variations, the Preisach model, and the Gustafson–Kessel (GK) approach. The choice of model depends on the desired level of accuracy and the availability of computational resources. A comparison of these models can be found in Table 2.

In the case of auxiliary parts simulation, it could be concluded that simplifying the object of study is crucial when using methods such as the finite element and finite volume methods. It is important to consider only the structural elements that can significantly impact the final result, such as the magnetic circuit and windings. However, this simplification can sometimes lead to inaccuracies. While additional elements may not be the primary source of heat generation, modeling these parts of the transformer can significantly affect the temperature field, making them important for design purposes and others.

Mathematical models for transformer thermal analysis can be broadly classified into two categories: circuit theory-based models and numerical field problem-solving models. The first group allows for real-time calculations but relies on empirical or analytical dependencies with multiple assumptions, limiting their accuracy. A simple example from this group is dynamic models, which have better accuracy but are comparable to those described in standards (IEC 60076-7, IEEE C57.91-2011). A more accurate method is the Thermal Network Model, which combines the advantages of simple equivalent circuits and models based on the finite element or volume method. This approach does not require significant computational resources, and its parameters are spatially dependent. However, it uses analytical and empirical dependencies for regions with liquid, and therefore, some phenomena associated with heat transfer are not taken into account. Another approach is to use numerical methods to solve partial differential equations for conjugate heat transfer, which requires solving the heat transfer equation for solid and fluid domains coupled with the Navier–Stokes equation. However, the full formulation for compressible flow requires significant computational resources, and in the case of oil-immersed transformers, a simplified Boussinesq’s approximation can be used. An alternative way to reduce the complexity of the numerical model is to consider heat transfer only in the solid parts of the power transformer, without calculating hydrodynamic processes using convective correlation based on empirical, semi-empirical, or calculated relationships of heat fluxes on the boundaries of solid parts. However, it should be noted that the approach using heat transfer coefficients may yield slightly distorted results, as heat transfer is non-uniformly distributed on the walls of the transformer winding, and the heat transfer coefficient varies greatly depending on temperature and other system parameters. Additionally, researchers often use a two-dimensional axisymmetric approximation of the problem to simplify the calculation process, but three-dimensional effects can significantly influence the results, such as the location and value of the hottest point. Three-dimensional models can be utilized to evaluate the design or optimal placement of measuring sensors, possessing high accuracy but a computationally intensive nature, which poses a challenge when creating digital twins that must operate in real time. A feasible solution in this scenario is to employ reduced-order models, which can be derived from the results of high-precision field calculations through regression analysis or projection onto a lower-dimensional space, significantly expediting the calculation process.

Based on the results of the studies, the recommendations regarding the selection of mathematical approaches can be derived. The choice of approach depends on the desired output results. If the goal is to predict integral parameters for traditional oil-immersed power transformer design configurations, analytical and circuit-based methods are suitable. These methods offer fast calculations and have minimal computational requirements. However, for more complex and distributed results, such as predicting moisture content and degree of polymerization, physical modeling using finite element and finite volume methods becomes necessary.

However, it should be noted that the list of the considered research is limited to the scope of this review. Creating a physics-based digital twin for an oil-immersed power transformer should cover many other aspects. Future research should delve into detailed approaches for constructing a digital twin and also take into account other physical interactions, including mechanical effects.

7. Conclusions

Based on the studies considered, it can be concluded that the primary defects in modern power transformers are often a result of insufficient information about their current technical state. To address this problem, digital twin (DT) technology proves to be useful. While data-driven approaches are valuable, the integration of physical models is crucial for creating a comprehensive digital twin. In this context, multiphysics models are necessary to describe all the processes occurring in a transformer. These models can be implemented using analytical dependencies, circuit theory-based models, or numerical solutions. Therefore, this review focuses on the most commonly used methods of mathematical modeling for electromagnetic, thermal, and hydrodynamic processes in oil-immersed power transformers. By adopting this approach, stakeholders in the power industry can gain valuable insights into the technical state monitoring of transformers. This enables proactive maintenance and minimizes the occurrence of primary defects. The analysis of various modeling methods presented in the review highlights the complexity of accounting for the interaction of different physical field. The analysis of various approaches to modeling physical processes in a transformer reveals that it is not feasible to create a universal model that accounts for all the features of physical processes in an oil-immersed transformer and operates in real time when developing a digital twin. This is due to the need for different levels of detail in the results according to different purposes. For instance, lumped models based on circuit theory can significantly reduce the calculation time, but they do not allow for the uneven distribution of heat dissipation in the transformer structural parts to be taken into account. Moreover, mathematical modeling is complicated by the interaction of different physical fields, which cannot be sufficiently described by analytical expressions due to their nature and non-linearity. Therefore, to implement digital twin technology, it is necessary to select the main physical processes that guarantee the adequacy of the reproduction of the real functioning of the oil-immersed power transformer.

The conducted research shows that the joint modeling of electromagnetic and thermal processes provides the most comprehensive solution to the problem. However, as demonstrated in this review, modeling individual electromagnetic, thermal, or hydrodynamic processes while accounting for the spatial distribution of the results can be computationally intensive. Therefore, a comparative analysis of methods for modeling physical processes in a transformer is conducted, taking into account their computational complexity and potential ways to reduce it. Based on the analysis results, we can conclude that reducing the dimensionality of each model separately, when used in combination, is a promising approach. At the same time, the decrease in accuracy of individual models with a decrease in their dimension and/or simplification is compensated for by an increase in the completeness and adequacy of reproducing the processes of the actual functioning of a power transformer. Indeed, combining the results and approaches from power transformer modeling with digital twin technology holds great promise. The active research in power transformer modeling, as is evident from the numerous publications in international journals, provides a solid foundation for integrating digital twin technology into the field. By leveraging the insights and advancements in power transformer modeling, a comprehensive digital twin can be developed to monitor and optimize the performance of transformers in real-time.

However, it is important to note that modeling the physical processes in a complex object such as a power transformer is a multifaceted task that requires further study. Moreover, there are numerous variations in power transformer design configurations, each of which necessitates an individualized approach to modeling. To specify the scope of this review, this study specifically focuses on the electromagnetic, thermal, and hydrodynamic processes of core-type oil-immersed power transformers with two-winding configurations, utilizing natural convection cooling.