# Optimization Design of Oil-Immersed Iron Core Reactor Based on the Particle Swarm Algorithm and Thermal Network Model

Authors: Fating Yuan; Kai Lv; Bo Tang; Yue Wang; Wentao Yang; Shihong Qin; Can Ding
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, according to the design parameters of oil-immersed iron core reactor, the thermal network model of windings is established by the thermo-electric analogy method, and the temperature distribution of the windings can be obtained. Meanwhile, a fluid-thermal coupled finite element model is established, the temperature and fluid velocity distribution are extracted, and the simulation results show that the error coefficient of temperature is less than 3% compared with the thermal network model, so the correctness of thermal network model has been verified. Taking the metal conductor usage and loss of windings as the optimization objects, the optimization method based on the particle swarm algorithm and thermal network model is proposed, and the Pareto optimal solutions between the metal conductor usage and loss of windings are given. The optimization results show that the metal conductor usage is reduced by 23.05%, and the loss is reduced by 20.25% compared with the initial design parameters, and the maximum temperature of winding does not exceed the expected value. Thus, the objects of low metal conductor usage and loss of windings are conflicted and cannot be optimized simultaneously; the optimization method has an important guiding significance for the design of oil-immersed iron core.


---

## Body

## 1. Introduction

The oil-immersed iron core reactor is an indispensable equipment in the power system, which plays the role of reactive power compensation, limiting the fault current and improving the power factor of the system [1]. However, the reactor encounters overheating and even burns out in the operation process, and the relevant studies show that the excessive temperature of the windings is the main reason [2]. The oil flow of the oil-immersed reactor enters the oil channels of the disc-type windings; the hot spot is formed as the convection heat dissipation capacity of each disc which is different; when the hot spot exceeds the limit, it will directly affect the safe and stable operation of the reactor. In order to reduce the temperature of reactor, the methods, such as increasing the cross-sectional area of conductor and oil channel width, are adopted; they can reduce the loss and temperature of the windings, but the metal conductor usage is increased. The loss and the metal conductor usage of windings are two key parameters in the design of reactor; in order to realize the low loss and the metal conductor usage simultaneously, the temperature calculation is needed and the optimization design method about the reactor is essential.Currently, the relevant studies mainly include (1) the temperature calculation of the reactor. The finite-difference method is used to calculate the temperature of the reactor [3, 4], and it is merely suitable for the situations that the temperature and heat flux distribution of the coil surface are the same, which limits its practical application. The finite element model is adopted [5–7], and the detailed temperature and fluid velocity distribution of the reactor can be obtained, but the calculation process is complex and calculation time is long. The thermal network model is put forward, and it was applied to the transformer windings [8]. On the basis, the optimization model of transformer windings is proposed [9–12]. In [13], the complex model of oil channels in windings is simplified to an approximate matrix of hydraulic channel, and the correctness has been verified. (2) The optimization design of the reactor. In [14], the finite element method is adopted to improve the heat dissipation efficiency of the windings by adjusting the oil channel width and the thickness of the coils, but the calculation time is large and the principle is not clear. In [15, 16], the thermal-electromagnetic combined optimization method is proposed to reduce the metal conductor, but the influence of air duct’s width on the heat dissipation efficiency is not considered. In [17–21], the genetic algorithm and particle swarm optimization are applied to minimize the metal conductor usage or the loss of the transformer, but it can only realize the optimization of a single target, the overall optimization of the reactor cannot be achieved.In this paper, taking the metal conductor usage and loss of windings as the optimization objects, the thermal network model of windings is established by the thermo-electric analogy method, the temperature distribution of the windings can be obtained, and the correctness has been verified by the fluid-thermal coupled finite element method. Meanwhile, the optimization method based on particle swarm algorithm and thermal network model is proposed, and the Pareto optimal solutions between the metal conductor usage and loss of windings are given. The optimization results show that the metal conductor usage is reduced by 23.05%, and the loss is reduced by 20.25% compared with the initial design parameters, and the maximum temperature of windings does not exceed the expected value, and the correctness has been verified.


## 2. The Structure and Equivalent Model of Windings

### 2.1. The Structure and Parameters of Core Reactor

The oil-immersed core reactor with a rated voltage of 26.4 kV and rated capacity of 66 MVA is selected as a research object. The disc-type windings are made of wires with insulating layer, which are surrounded on the iron core column. In the windings, the horizontal oil channels are formed between disc-type windings in the axial direction; the inside and outside of the windings are the vertical insulation paper tubes. Meanwhile, a baffle for oil is set up at intervals of a certain number of discs, so as to lead the oil flow in the vertical direction to the horizontal direction, which can enhance the cooling effect of windings. According to the structure characteristic of iron reactor, the influence of iron column and iron yoke on the heat dissipation of windings can be ignored; thus, the windings of iron reactor are only considered in this paper, and the basic structure of windings is given in Figure1.Figure 1
The basic structure of windings.


### 2.2. The Equivalent Model of Windings

The total windings consist of 60 layers of disc, and they were divided into 4 channels by the oil guiding washer. Each channel consists of 15 layers of disc, each layer of disc consists of 23 turns of conductor with insulating layers, and the thickness of the insulating layer is 0.4 mm. Considering the temperature and fluid velocity distribution in each channel of the windings are basically the same, thus, one of the channels is selected to research the temperature distribution of windings, and the equivalent model and parameters are given in Figure2.Figure 2
The equivalent model of windings.The characteristics of metal conductor and insulating materials are given in Table1.Table 1
The characteristics of metal conductor and insulating materials.
MaterialsAttributeValueCopperHeat conductivity (W·M−1·K−1)400Specific heat (W·kg−1·K−1)385Density (Kg M−3)8940Insulating materialsHeat conductivity (W·M−1·K−1)0.21Specific heat (W·kg−1·K−1)1250Density (Kg M−3)870


## 2.1. The Structure and Parameters of Core Reactor

The oil-immersed core reactor with a rated voltage of 26.4 kV and rated capacity of 66 MVA is selected as a research object. The disc-type windings are made of wires with insulating layer, which are surrounded on the iron core column. In the windings, the horizontal oil channels are formed between disc-type windings in the axial direction; the inside and outside of the windings are the vertical insulation paper tubes. Meanwhile, a baffle for oil is set up at intervals of a certain number of discs, so as to lead the oil flow in the vertical direction to the horizontal direction, which can enhance the cooling effect of windings. According to the structure characteristic of iron reactor, the influence of iron column and iron yoke on the heat dissipation of windings can be ignored; thus, the windings of iron reactor are only considered in this paper, and the basic structure of windings is given in Figure1.Figure 1
The basic structure of windings.


## 2.2. The Equivalent Model of Windings

The total windings consist of 60 layers of disc, and they were divided into 4 channels by the oil guiding washer. Each channel consists of 15 layers of disc, each layer of disc consists of 23 turns of conductor with insulating layers, and the thickness of the insulating layer is 0.4 mm. Considering the temperature and fluid velocity distribution in each channel of the windings are basically the same, thus, one of the channels is selected to research the temperature distribution of windings, and the equivalent model and parameters are given in Figure2.Figure 2
The equivalent model of windings.The characteristics of metal conductor and insulating materials are given in Table1.Table 1
The characteristics of metal conductor and insulating materials.
MaterialsAttributeValueCopperHeat conductivity (W·M−1·K−1)400Specific heat (W·kg−1·K−1)385Density (Kg M−3)8940Insulating materialsHeat conductivity (W·M−1·K−1)0.21Specific heat (W·kg−1·K−1)1250Density (Kg M−3)870


## 3. The Temperature Calculation of Windings Based on the Thermal Network Model

### 3.1. The Establishment of Thermal Network Model

The heat transfer model of windings can be equivalent to the circuit node network, which adopts the analogy method between the circuit and the heat transfer system. The basic idea of thermo-electric analogy is to compare the transfer of heat between media to the transfer of electric current in a conductor. The related physical essence of the temperature field is analogized to the heat transfer parameters, such as temperature difference, heat flux, heat resistance, and heat capacity. So, the topology of the thermal network can be established, and it is similar to the circuit [22]. In order to obtain the temperature distribution of each conductor and oil flow node, a thermal network model with each conductor as the basic unit is established based on the equivalent model of windings, as shown in Figure 3.Figure 3
The thermal network model.In Figure3, Tdisc and Toil are the node temperature of conductor and oil flow, respectively, Rcir is the conduction thermal resistance between conductors and insulating layers, Roor is the conduction thermal resistance of oil flow, Rior is the convection thermal resistance between the insulation layers and oil, and q is the heat source of the conductor.


### 3.2. The Calculation of Thermal Resistance

In the actual operation process of oil-immersed iron core reactor, the influence of the leakage magnetic on the eddy current loss of each turn is basically the same, so it is considered that each turn conductor has the same heat source. Considering the temperature difference between the adjacent conductors is small, the radiation heat transfer can be ignored, so the heat transfer modes of heat flux in the four directions of the conductor are merely the heat conduction and heat convection. In order to establish the thermal network model, the calculation of thermal resistance is essential.(1)
The conduction thermal resistance: in the inner of conductors and insulating material, the conduction thermal resistance depends on the geometric structure and material properties. When the heat flux is transferred to the oil channel, the heat transfer will occur in the oil flow, and the conduction thermal resistance of the oil flow is related to the geometric shape of the oil channel and physical properties of the oil. The conduction thermal resistance in the conductors, insulation layers and oil can be calculated, as shown in the following equation:(1)Rcir/oor=LqkmS,whereLq is the thickness of the conductor in the direction in which the heat flows, km is the thermal conductivity of material, S is the surface area, which is perpendicular to heat flow direction, Rcir is the conduction thermal resistance of metal conductor and insulating layers, and Roor is the conduction thermal resistance of oil.(2)
The convection thermal resistance: the convection heat transfer is formed between the discs and the horizontal and vertical oil channel, and the thermal resistance is related to the fluid velocity and temperature of oil. According to the Newton’s law of cooling, the convection thermal resistance between the surfaces of the insulating layers and the oil channels can be expressed, as shown in the following equation:(2)Rior=1αA,whereRior is the convection thermal resistance between the insulation layers and oil, A is the area of the insulating layer, and α is the convection coefficient, which is related to the geometric structure, material, and fluid characteristics of oil, and it can be written as follows:(3)α=NukoilLio,whereLio is the geometric characteristic length of the convection heat transfer surface, koil is the thermal conductivity of transformer oil, and Nu is the Nusselt number, and the expression is [23](4)Nu=1.86RePrDL1/3μμs0.14,whereL is the horizontal length, D is the diameter of oil channel, μ is the kinetic viscosity of oil with the ambient temperatures, μs is the kinetic viscosity of oil, Re is the Reynolds number, and Pr is the Prandtl number. The parameters, Re and Pr, can be calculated by the following equation:(5)Re=ρoilviDμs,Pr=cpμskoil,whereρoil and cp are the density and specific heat capacity of oil, respectively, and vi is the fluid velocity of horizontal oil channel in the ith layer.


### 3.3. The Temperature Calculation of Windings

According to the thermal network model in Figure3, it contains N nodes, so the n-1 node equations can be listed [21], which is similar to the node voltage equation in the circuit. The temperature equations of each node can be written [24], as shown in the following equation:(6)∑1Rself−impedanceTdisc,i,j−∑1Rmutual−impedanceTdisc,mi,j−∑1Rmutual−impedanceToil,i,nj=q,∑1Rself−impedanceToil,i,j−∑1Rmutual−impedanceTdisc,mi,j−∑1Rmutual−impedanceToil,i,nj=0,where Rself−impedance is the node self-thermal resistance, Rmutual−impedance is the node mutual thermal resistance, Tdisc,i,j and Toil,i,j are the conductor and oil flow nodes of the ith row and jth column in the model, respectively, mi and nj are the numbers of adjacent nodes, mi=i−1,i+1, nj=j−1,j+1. The above equation can be solved by the MATLAB software, and the node temperature can be expressed as(7)T=G−1Q,where T is the temperature of node, Q is the heat source matrix of the conductor, and G is the heat conductivity matrix, which is composed by the self-thermal resistance and mutual thermal resistance, as shown in the following equation:(8)G=1Rsa_1−1Rma1_10⋅−1Rma2_1⋅00−1Rma1_11Rsa_2−1Rma1_20⋅⋅⋅00−1Rma1_21Rsa_2−1Rma1_3⋅⋅⋅⋅⋅0−1Rma1_3⋅⋅⋅⋅−1Rma2_t−1Rma2_1⋅⋅⋅⋅⋅0⋅⋅⋅⋅⋅⋅1Rsa_n−2−1Rma1_n−100⋅⋅⋅0−1Rma1_n−11Rsa_n−1−1Rma1_n00⋅−1Rma2_t⋅0−1Rma1_n1Rsa_n,where, 1/Rma1_i and 1/Rma2_j are the sum of mutual admittance of the ith and jth nodes, respectively, and the parameter, 1/Rsa_i, is the sum of self-admittance of the node.According to the structure and electrical parameters of windings, combined with the above calculation method, the temperature of each node can be obtained by solving equation (7).


### 3.4. The Calculation Results Based on Thermal Network Model

The settings of thermal network model are as follows: the loss of each turn conductor is 34.5 W, the temperature in the inlet of the oil channel is 303.15 K, and the velocity is 0.3 m/s; the boundary conditions of the upper and lower surfaces of each channel are heat insulation as the low thermal conductivity of the oil guide baffle. In the thermal network model, there exist 437 conductor nodes and 460 oil channel nodes; the number of the discs is sorted from top to bottom; combined with the thermal resistance calculation method and node temperature equation, the temperature of each node can be obtained, and the maximum temperature of each disc is given, as shown in Table2.Table 2
The maximum temperature of each disc.
Disc number12345678Temperature (℃)40.7945.5850.4755.4959.3361.3562.3262.65Disc number9101112131415Temperature (℃)62.5662.1761.5460.7359.7758.6857.47In Table2, it can be seen that the maximum temperature of discs is gradually increased from the bottom to the top, but there is a downward trend when nearing the top, the maximum temperature is 62.65°C, and the maximum temperature area is mainly concentrated in the upper middle position of winding; the main reason is that the fluid velocity of oil is large in the inlet and outlet, and the convection heat transfer degree is high. When the oil flows in the long channel, the temperature is increased and the fluid velocity of oil is decreased, which leads to the accumulation of heat in the middle of windings, so the temperature of the middle discs is higher than the upper and lower discs.


## 3.1. The Establishment of Thermal Network Model

The heat transfer model of windings can be equivalent to the circuit node network, which adopts the analogy method between the circuit and the heat transfer system. The basic idea of thermo-electric analogy is to compare the transfer of heat between media to the transfer of electric current in a conductor. The related physical essence of the temperature field is analogized to the heat transfer parameters, such as temperature difference, heat flux, heat resistance, and heat capacity. So, the topology of the thermal network can be established, and it is similar to the circuit [22]. In order to obtain the temperature distribution of each conductor and oil flow node, a thermal network model with each conductor as the basic unit is established based on the equivalent model of windings, as shown in Figure 3.Figure 3
The thermal network model.In Figure3, Tdisc and Toil are the node temperature of conductor and oil flow, respectively, Rcir is the conduction thermal resistance between conductors and insulating layers, Roor is the conduction thermal resistance of oil flow, Rior is the convection thermal resistance between the insulation layers and oil, and q is the heat source of the conductor.


## 3.2. The Calculation of Thermal Resistance

In the actual operation process of oil-immersed iron core reactor, the influence of the leakage magnetic on the eddy current loss of each turn is basically the same, so it is considered that each turn conductor has the same heat source. Considering the temperature difference between the adjacent conductors is small, the radiation heat transfer can be ignored, so the heat transfer modes of heat flux in the four directions of the conductor are merely the heat conduction and heat convection. In order to establish the thermal network model, the calculation of thermal resistance is essential.(1)
The conduction thermal resistance: in the inner of conductors and insulating material, the conduction thermal resistance depends on the geometric structure and material properties. When the heat flux is transferred to the oil channel, the heat transfer will occur in the oil flow, and the conduction thermal resistance of the oil flow is related to the geometric shape of the oil channel and physical properties of the oil. The conduction thermal resistance in the conductors, insulation layers and oil can be calculated, as shown in the following equation:(1)Rcir/oor=LqkmS,whereLq is the thickness of the conductor in the direction in which the heat flows, km is the thermal conductivity of material, S is the surface area, which is perpendicular to heat flow direction, Rcir is the conduction thermal resistance of metal conductor and insulating layers, and Roor is the conduction thermal resistance of oil.(2)
The convection thermal resistance: the convection heat transfer is formed between the discs and the horizontal and vertical oil channel, and the thermal resistance is related to the fluid velocity and temperature of oil. According to the Newton’s law of cooling, the convection thermal resistance between the surfaces of the insulating layers and the oil channels can be expressed, as shown in the following equation:(2)Rior=1αA,whereRior is the convection thermal resistance between the insulation layers and oil, A is the area of the insulating layer, and α is the convection coefficient, which is related to the geometric structure, material, and fluid characteristics of oil, and it can be written as follows:(3)α=NukoilLio,whereLio is the geometric characteristic length of the convection heat transfer surface, koil is the thermal conductivity of transformer oil, and Nu is the Nusselt number, and the expression is [23](4)Nu=1.86RePrDL1/3μμs0.14,whereL is the horizontal length, D is the diameter of oil channel, μ is the kinetic viscosity of oil with the ambient temperatures, μs is the kinetic viscosity of oil, Re is the Reynolds number, and Pr is the Prandtl number. The parameters, Re and Pr, can be calculated by the following equation:(5)Re=ρoilviDμs,Pr=cpμskoil,whereρoil and cp are the density and specific heat capacity of oil, respectively, and vi is the fluid velocity of horizontal oil channel in the ith layer.


## 3.3. The Temperature Calculation of Windings

According to the thermal network model in Figure3, it contains N nodes, so the n-1 node equations can be listed [21], which is similar to the node voltage equation in the circuit. The temperature equations of each node can be written [24], as shown in the following equation:(6)∑1Rself−impedanceTdisc,i,j−∑1Rmutual−impedanceTdisc,mi,j−∑1Rmutual−impedanceToil,i,nj=q,∑1Rself−impedanceToil,i,j−∑1Rmutual−impedanceTdisc,mi,j−∑1Rmutual−impedanceToil,i,nj=0,where Rself−impedance is the node self-thermal resistance, Rmutual−impedance is the node mutual thermal resistance, Tdisc,i,j and Toil,i,j are the conductor and oil flow nodes of the ith row and jth column in the model, respectively, mi and nj are the numbers of adjacent nodes, mi=i−1,i+1, nj=j−1,j+1. The above equation can be solved by the MATLAB software, and the node temperature can be expressed as(7)T=G−1Q,where T is the temperature of node, Q is the heat source matrix of the conductor, and G is the heat conductivity matrix, which is composed by the self-thermal resistance and mutual thermal resistance, as shown in the following equation:(8)G=1Rsa_1−1Rma1_10⋅−1Rma2_1⋅00−1Rma1_11Rsa_2−1Rma1_20⋅⋅⋅00−1Rma1_21Rsa_2−1Rma1_3⋅⋅⋅⋅⋅0−1Rma1_3⋅⋅⋅⋅−1Rma2_t−1Rma2_1⋅⋅⋅⋅⋅0⋅⋅⋅⋅⋅⋅1Rsa_n−2−1Rma1_n−100⋅⋅⋅0−1Rma1_n−11Rsa_n−1−1Rma1_n00⋅−1Rma2_t⋅0−1Rma1_n1Rsa_n,where, 1/Rma1_i and 1/Rma2_j are the sum of mutual admittance of the ith and jth nodes, respectively, and the parameter, 1/Rsa_i, is the sum of self-admittance of the node.According to the structure and electrical parameters of windings, combined with the above calculation method, the temperature of each node can be obtained by solving equation (7).


## 3.4. The Calculation Results Based on Thermal Network Model

The settings of thermal network model are as follows: the loss of each turn conductor is 34.5 W, the temperature in the inlet of the oil channel is 303.15 K, and the velocity is 0.3 m/s; the boundary conditions of the upper and lower surfaces of each channel are heat insulation as the low thermal conductivity of the oil guide baffle. In the thermal network model, there exist 437 conductor nodes and 460 oil channel nodes; the number of the discs is sorted from top to bottom; combined with the thermal resistance calculation method and node temperature equation, the temperature of each node can be obtained, and the maximum temperature of each disc is given, as shown in Table2.Table 2
The maximum temperature of each disc.
Disc number12345678Temperature (℃)40.7945.5850.4755.4959.3361.3562.3262.65Disc number9101112131415Temperature (℃)62.5662.1761.5460.7359.7758.6857.47In Table2, it can be seen that the maximum temperature of discs is gradually increased from the bottom to the top, but there is a downward trend when nearing the top, the maximum temperature is 62.65°C, and the maximum temperature area is mainly concentrated in the upper middle position of winding; the main reason is that the fluid velocity of oil is large in the inlet and outlet, and the convection heat transfer degree is high. When the oil flows in the long channel, the temperature is increased and the fluid velocity of oil is decreased, which leads to the accumulation of heat in the middle of windings, so the temperature of the middle discs is higher than the upper and lower discs.


## 4. Model Validation

### 4.1. The Finite Element Simulation Calculation

According to the structure parameters and material characteristics of the windings, the two-dimensional model is established based on the COMSOL, as shown in Figure4.Figure 4
The simulation model of windings.The governing equations and boundary conditions settings are the key steps of simulation calculation.


#### 4.1.1. Governing Equation

The heat conduction and heat convection are the main heat transfer processes in the windings region, and the governing equation of heat conduction can be expressed by the following equation [25]:(9)∂2T∂r2+∂2T∂z2+qkm=0.The heat convection between windings and oil channels satisfies mass, momentum, and energy conservation, as shown in the following equation [26]:(10)1r∂ur∂r+1z∂vr∂z=0,ρoilu∂u∂r+v∂u∂z=Fr−∂p∂r+μs∂2u∂r2+∂2u∂z2,ρoilu∂v∂r+v∂v∂z=Fz−∂p∂z+μs∂2v∂r2+∂2v∂z2,u∂T∂r+v∂T∂z=koilρoilcp∂2v∂r2+∂2v∂z2,where u and v are the radial and axial fluid velocities of oil, respectively, Fr and Fz are the radial and axial components of fluid gravity, and p is the pressure of fluid.


#### 4.1.2. The Boundary Conditions Setting

The boundary conditions are set as follows: (1) the heat transfer between the insulating paper and the outer transformer oil can be ignored, and it sets as an adiabatic boundary. (2) The discs surface is regarded as the stationary, the axial and radial velocity are zero, the inlet velocity of oil is 0.3 m/s, and the temperature is 303.15 K. (3) The outlet boundary condition of oil flow is natural pressure outlet boundary and the axis of rotation is axisymmetric boundary condition. The heat source of each turn conductor is setting the same; the value is 34.5 W. Meanwhile, the characteristic parameters of oil are given in Table3.Table 3
The characteristics of oil.
MaterialsAttributeValueOilHeat conductivity (W·M−1·K−1)0.134 −8.05 × 10 − 5TSpecific heat (W·kg−1·K−1)−13408.15 + 123.04T − 0.33T ^ 2Density (Kg M−3)1055.05 − 0.58T − 6.4 × 10 − 5T ^ 2Dynamic viscosity (kg·M−1·S−1)91.45 − 1.33T + 0.01T ^ 2


### 4.2. The Simulation Results

The transient field is applied in the simulation, the total calculation time is 10 h, and the time step is 0.1 h, which can be considered that it has reached the steady state. The temperature of windings and fluid velocity distribution of oil are given in Figure5.Figure 5
Simulation results: (a) temperature distribution; (b) fluid velocity distribution.
(a)(b)In Figure5, the temperature of discs in the inlet is lower than in the middle, and the maximum temperature reaches 66.20°C in the 8th disc; the fluid velocity of oil is relatively large in the inlet and outlet. In order to analyze the heat transfer process of discs, the radial and axial direction paths are extracted, as shown in Figure 6.Figure 6
The extract path of the model.According to the paths selected and simulation results, the temperature distribution with the different paths is plotted, as shown in Figures7–9.Figure 7
The temperature rise distribution along the radial paths.Figure 8
The temperature rise distribution along the axial paths.Figure 9
The fluid velocity distribution along the horizontal paths.From Figures7–9, it can be seen that the temperature of discs is increased first and then decreased in axial and radial directions. On the contrary, the fluid velocity along the axial direction is decreased first and then increased, which is opposite to the trends of temperature. Meanwhile, the maximum temperature locates the eighth disc, and the fluid velocity is the lowest on the both sides of the discs with the maximum temperature, the main reasons are that the fluid velocity in lower channels is relatively large when the oil passes the inlet, and the temperature is low. The heat dissipation conditions are weak because of the high temperature and low velocity of the oil in the middle discs; the temperature of the middle discs is the highest. In the upper discs, the oil converges in the outlet; it leads to the increase of the fluid velocity, so the temperature of the upper discs is lower compared with the middle discs.


### 4.3. Model Validation

In order to verify the accuracy of the thermal network model, the maximum temperature of each disc is obtained in both the thermal network model and finite element method, as shown in Figure10.Figure 10
The maximum temperature of each disc.In Figure10, it can be seen that the maximum temperature is basically the same based on the thermal network model and finite element method, and the maximum error coefficient of the disc is only 5.30%, which locates the 8th disc, so the correctness of the thermal network model is verified.The error coefficient is defined as follows:(11)ε=Thot−TCOMTCOM,where ε is the error coefficient, Thot is the temperature based on the thermal network model, and TCOM is the temperature based on the finite element method. Meanwhile, the thermal network model is a multivariate linear equation related to thermal resistance, and the calculation speed is fast compared with the finite element method; thus, the thermal network model is adopted to optimize the windings in this paper.


## 4.1. The Finite Element Simulation Calculation

According to the structure parameters and material characteristics of the windings, the two-dimensional model is established based on the COMSOL, as shown in Figure4.Figure 4
The simulation model of windings.The governing equations and boundary conditions settings are the key steps of simulation calculation.


### 4.1.1. Governing Equation

The heat conduction and heat convection are the main heat transfer processes in the windings region, and the governing equation of heat conduction can be expressed by the following equation [25]:(9)∂2T∂r2+∂2T∂z2+qkm=0.The heat convection between windings and oil channels satisfies mass, momentum, and energy conservation, as shown in the following equation [26]:(10)1r∂ur∂r+1z∂vr∂z=0,ρoilu∂u∂r+v∂u∂z=Fr−∂p∂r+μs∂2u∂r2+∂2u∂z2,ρoilu∂v∂r+v∂v∂z=Fz−∂p∂z+μs∂2v∂r2+∂2v∂z2,u∂T∂r+v∂T∂z=koilρoilcp∂2v∂r2+∂2v∂z2,where u and v are the radial and axial fluid velocities of oil, respectively, Fr and Fz are the radial and axial components of fluid gravity, and p is the pressure of fluid.


### 4.1.2. The Boundary Conditions Setting

The boundary conditions are set as follows: (1) the heat transfer between the insulating paper and the outer transformer oil can be ignored, and it sets as an adiabatic boundary. (2) The discs surface is regarded as the stationary, the axial and radial velocity are zero, the inlet velocity of oil is 0.3 m/s, and the temperature is 303.15 K. (3) The outlet boundary condition of oil flow is natural pressure outlet boundary and the axis of rotation is axisymmetric boundary condition. The heat source of each turn conductor is setting the same; the value is 34.5 W. Meanwhile, the characteristic parameters of oil are given in Table3.Table 3
The characteristics of oil.
MaterialsAttributeValueOilHeat conductivity (W·M−1·K−1)0.134 −8.05 × 10 − 5TSpecific heat (W·kg−1·K−1)−13408.15 + 123.04T − 0.33T ^ 2Density (Kg M−3)1055.05 − 0.58T − 6.4 × 10 − 5T ^ 2Dynamic viscosity (kg·M−1·S−1)91.45 − 1.33T + 0.01T ^ 2


## 4.1.1. Governing Equation

The heat conduction and heat convection are the main heat transfer processes in the windings region, and the governing equation of heat conduction can be expressed by the following equation [25]:(9)∂2T∂r2+∂2T∂z2+qkm=0.The heat convection between windings and oil channels satisfies mass, momentum, and energy conservation, as shown in the following equation [26]:(10)1r∂ur∂r+1z∂vr∂z=0,ρoilu∂u∂r+v∂u∂z=Fr−∂p∂r+μs∂2u∂r2+∂2u∂z2,ρoilu∂v∂r+v∂v∂z=Fz−∂p∂z+μs∂2v∂r2+∂2v∂z2,u∂T∂r+v∂T∂z=koilρoilcp∂2v∂r2+∂2v∂z2,where u and v are the radial and axial fluid velocities of oil, respectively, Fr and Fz are the radial and axial components of fluid gravity, and p is the pressure of fluid.


## 4.1.2. The Boundary Conditions Setting

The boundary conditions are set as follows: (1) the heat transfer between the insulating paper and the outer transformer oil can be ignored, and it sets as an adiabatic boundary. (2) The discs surface is regarded as the stationary, the axial and radial velocity are zero, the inlet velocity of oil is 0.3 m/s, and the temperature is 303.15 K. (3) The outlet boundary condition of oil flow is natural pressure outlet boundary and the axis of rotation is axisymmetric boundary condition. The heat source of each turn conductor is setting the same; the value is 34.5 W. Meanwhile, the characteristic parameters of oil are given in Table3.Table 3
The characteristics of oil.
MaterialsAttributeValueOilHeat conductivity (W·M−1·K−1)0.134 −8.05 × 10 − 5TSpecific heat (W·kg−1·K−1)−13408.15 + 123.04T − 0.33T ^ 2Density (Kg M−3)1055.05 − 0.58T − 6.4 × 10 − 5T ^ 2Dynamic viscosity (kg·M−1·S−1)91.45 − 1.33T + 0.01T ^ 2


## 4.2. The Simulation Results

The transient field is applied in the simulation, the total calculation time is 10 h, and the time step is 0.1 h, which can be considered that it has reached the steady state. The temperature of windings and fluid velocity distribution of oil are given in Figure5.Figure 5
Simulation results: (a) temperature distribution; (b) fluid velocity distribution.
(a)(b)In Figure5, the temperature of discs in the inlet is lower than in the middle, and the maximum temperature reaches 66.20°C in the 8th disc; the fluid velocity of oil is relatively large in the inlet and outlet. In order to analyze the heat transfer process of discs, the radial and axial direction paths are extracted, as shown in Figure 6.Figure 6
The extract path of the model.According to the paths selected and simulation results, the temperature distribution with the different paths is plotted, as shown in Figures7–9.Figure 7
The temperature rise distribution along the radial paths.Figure 8
The temperature rise distribution along the axial paths.Figure 9
The fluid velocity distribution along the horizontal paths.From Figures7–9, it can be seen that the temperature of discs is increased first and then decreased in axial and radial directions. On the contrary, the fluid velocity along the axial direction is decreased first and then increased, which is opposite to the trends of temperature. Meanwhile, the maximum temperature locates the eighth disc, and the fluid velocity is the lowest on the both sides of the discs with the maximum temperature, the main reasons are that the fluid velocity in lower channels is relatively large when the oil passes the inlet, and the temperature is low. The heat dissipation conditions are weak because of the high temperature and low velocity of the oil in the middle discs; the temperature of the middle discs is the highest. In the upper discs, the oil converges in the outlet; it leads to the increase of the fluid velocity, so the temperature of the upper discs is lower compared with the middle discs.


## 4.3. Model Validation

In order to verify the accuracy of the thermal network model, the maximum temperature of each disc is obtained in both the thermal network model and finite element method, as shown in Figure10.Figure 10
The maximum temperature of each disc.In Figure10, it can be seen that the maximum temperature is basically the same based on the thermal network model and finite element method, and the maximum error coefficient of the disc is only 5.30%, which locates the 8th disc, so the correctness of the thermal network model is verified.The error coefficient is defined as follows:(11)ε=Thot−TCOMTCOM,where ε is the error coefficient, Thot is the temperature based on the thermal network model, and TCOM is the temperature based on the finite element method. Meanwhile, the thermal network model is a multivariate linear equation related to thermal resistance, and the calculation speed is fast compared with the finite element method; thus, the thermal network model is adopted to optimize the windings in this paper.


## 5. The Optimization Design of Windings Based on Particle Swarm Algorithm and Thermal Network Model

According to actual engineering requirements of the oil-immersed iron core reactor, the metal conductor usage and loss of windings are the two major optimization objects; however, these goals are in conflict with each other, the performance improvement of one object results in the reduction of the other object. In order to achieve the optimal design parameters of windings, the optimization design of iron reactor is needed.


### 5.1. The Particle Swarm Optimization Algorithm

The particle swarm optimization is one of the most influential multiobjective algorithms, including the prominent characteristics: the optimal solution of multiobjective can be easily obtained for its efficient search ability, and the algorithm has good multipeak search ability, which means it has a good ability to search for a global optimal value and it converges quickly. The flow chart of particle swarm optimization method is given in Figure11.Figure 11
The flow chart of particle swarm optimization algorithm.The formulas of particle update in particle swarm optimization are as follows [27]:(12)Vk+1=wVk+c1r1Pik−Xk+c2r2Pgk−Xk,Xk+1=Xk+Vk+1,where w is the inertia weight, r1,r2 are the random numbers distributed in [0, 1], k is the current iteration number, Pi is the individual optimal particle position, Pg is the global particle position, c1,c2 are the acceleration constants, V is the velocity of the particle, and X is the position of the particle.The optimization object is to obtain the minimum metal conductor usage and loss of windings; according to the structural characteristics of the windings, the size of the conductorWcoil, the height of horizontal channels Hcoil, and the thickness of insulation layer Lins are selected as the optimization variables. The object functions of metal conductor usage Mcu and loss Ploss can be written as(13)minMcu=πρm∑i=1ndiscdiWcoilHcoil,minPloss=πI2γ∑i=1ndiscdiWcoilHcoil,where ρ is the density of metal conductor, m is the number of discs, ndisc is the number of conductors in each disc, di is the average diameter of the ith conductor, I is the current of the conductor, and γ is the resistance of the conductor.The inequality constraints are(14)THot_max≤Tmax,Lmin≤L≤Lmax,Wcoil_min≤Wcoil≤Wcoil_max,Hcoil_min≤Hcoil≤Hcoil_max,Hduct_min≤Hduct≤Hduct_max,where THot_max is the maximum temperature of discs by calculation methods, Tmax is the maximum temperature constraints of the discs, L is the inductance of the reactor, and Wcoil_min, Wcoil_max, Hcoil_min, Hcoil_max, Hduct_min, and Hduct_max are the minimum and the maximum values of variables.


### 5.2. The Optimization Results

In particle swarm optimization algorithm, the variable parameters are set as follows: the radial width of the conductor is [1.6, 2.6] mm, the axial height of the conductor is [12, 14.6] mm, the thickness of the insulation layer is [0.25, 0.5] mm, and the height of horizontal channels is [3.9, 5.4] mm. The expected inductance difference is less than 2%, and the maximum temperature of winding is 75°C. The settings of particles in particle swarm optimization are as follows: the acceleration constantsc1 and c2 are 0.8, the iteration time is 200, and the inertia weight value is 0.8–1.2. The Pareto optimum results are shown in Figure 12,Figure 12
Pareto Frontier of particle swarm optimization.The change rate of metal conductor usage and loss,KM and Kloss, is defined as the following equation:(15)Kx=xx0,where x0 is the initial design parameter and x is the design parameter with the optimization method. In the initial parameters, the metal conductor usage and loss of windings are 90.98 kg and 11.9 kW before optimization, respectively.In Figure12, it can be seen that there exists an opposite trend between the metal conductor usage and loss of windings. Combined with the particle swarm optimization algorithm and thermal network model, the representative results can be obtained, as shown in Table 4. The inductance change rate is less than 1% compared with the initial parameters.Table 4
The pareto optimum solutions.
CaseWcoil (mm)Hcoil (mm)Lins (mm)Hduct (mm)Loss (kW)Weight (kg)11.613.60.354.515.0475.5722.612.10.405.111.06109.2531.613.10.454.515.6174.4541.812.10.454.815.2175.6452.414.60.404.29.81121.6861.614.10.503.914.5178.3471.813.60.403.913.5485.0181.813.10.354.814.0581.8991.812.60.353.914.6178.76101.813.10.504.214.0581.89111.612.60.504.216.2370.01122.012.10.504.213.8784.04132.413.10.404.210.94109.18142.614.10.254.29.49127.31151.813.10.504.514.0681.89162.414.60.403.99.81121.68171.812.60.354.214.6178.76182.012.10.255.113.8684.04192.012.60.305.113.3278.75201.613.10.303.915.6172.79In order to verify the correctness of the optimization results, taking the groups 2, 9, 17, and 20 in Table4 as an example, the maximum temperatures of windings are 54.91°C, 67.74°C, 68.56°C, and 73.40°C, respectively, based on the thermal network model. Meanwhile, the simulation results based on the finite element method are given in Figure 13; the maximum temperatures of windings are 55.10°C, 68.20°C, 70.00°C, and 72.63°C respectively; it can be seen that the temperature distribution laws of the four groups are basically the same, and the error coefficient is less than 3%; thus, the correctness of the optimization method has been verified.Figure 13
Simulation results: (a) Case 2; (b) Case 9; (c) Case 17; (d) Case 20.
(a)(b)(c)(d)In Table4, it can be seen that the metal conductor usage of windings is reduced by 23.05%, and the loss of the windings is reduced by 20.25% compared with the initial design parameters. As the purpose of low metal conductor usage and loss of winding are conflicted and cannot be optimized simultaneously, therefore, in the actual design process of the iron reactor, it is necessary to select the suitable structural parameter to realize the optimization of one objective with the other objectives which are not deteriorating.


## 5.1. The Particle Swarm Optimization Algorithm

The particle swarm optimization is one of the most influential multiobjective algorithms, including the prominent characteristics: the optimal solution of multiobjective can be easily obtained for its efficient search ability, and the algorithm has good multipeak search ability, which means it has a good ability to search for a global optimal value and it converges quickly. The flow chart of particle swarm optimization method is given in Figure11.Figure 11
The flow chart of particle swarm optimization algorithm.The formulas of particle update in particle swarm optimization are as follows [27]:(12)Vk+1=wVk+c1r1Pik−Xk+c2r2Pgk−Xk,Xk+1=Xk+Vk+1,where w is the inertia weight, r1,r2 are the random numbers distributed in [0, 1], k is the current iteration number, Pi is the individual optimal particle position, Pg is the global particle position, c1,c2 are the acceleration constants, V is the velocity of the particle, and X is the position of the particle.The optimization object is to obtain the minimum metal conductor usage and loss of windings; according to the structural characteristics of the windings, the size of the conductorWcoil, the height of horizontal channels Hcoil, and the thickness of insulation layer Lins are selected as the optimization variables. The object functions of metal conductor usage Mcu and loss Ploss can be written as(13)minMcu=πρm∑i=1ndiscdiWcoilHcoil,minPloss=πI2γ∑i=1ndiscdiWcoilHcoil,where ρ is the density of metal conductor, m is the number of discs, ndisc is the number of conductors in each disc, di is the average diameter of the ith conductor, I is the current of the conductor, and γ is the resistance of the conductor.The inequality constraints are(14)THot_max≤Tmax,Lmin≤L≤Lmax,Wcoil_min≤Wcoil≤Wcoil_max,Hcoil_min≤Hcoil≤Hcoil_max,Hduct_min≤Hduct≤Hduct_max,where THot_max is the maximum temperature of discs by calculation methods, Tmax is the maximum temperature constraints of the discs, L is the inductance of the reactor, and Wcoil_min, Wcoil_max, Hcoil_min, Hcoil_max, Hduct_min, and Hduct_max are the minimum and the maximum values of variables.


## 5.2. The Optimization Results

In particle swarm optimization algorithm, the variable parameters are set as follows: the radial width of the conductor is [1.6, 2.6] mm, the axial height of the conductor is [12, 14.6] mm, the thickness of the insulation layer is [0.25, 0.5] mm, and the height of horizontal channels is [3.9, 5.4] mm. The expected inductance difference is less than 2%, and the maximum temperature of winding is 75°C. The settings of particles in particle swarm optimization are as follows: the acceleration constantsc1 and c2 are 0.8, the iteration time is 200, and the inertia weight value is 0.8–1.2. The Pareto optimum results are shown in Figure 12,Figure 12
Pareto Frontier of particle swarm optimization.The change rate of metal conductor usage and loss,KM and Kloss, is defined as the following equation:(15)Kx=xx0,where x0 is the initial design parameter and x is the design parameter with the optimization method. In the initial parameters, the metal conductor usage and loss of windings are 90.98 kg and 11.9 kW before optimization, respectively.In Figure12, it can be seen that there exists an opposite trend between the metal conductor usage and loss of windings. Combined with the particle swarm optimization algorithm and thermal network model, the representative results can be obtained, as shown in Table 4. The inductance change rate is less than 1% compared with the initial parameters.Table 4
The pareto optimum solutions.
CaseWcoil (mm)Hcoil (mm)Lins (mm)Hduct (mm)Loss (kW)Weight (kg)11.613.60.354.515.0475.5722.612.10.405.111.06109.2531.613.10.454.515.6174.4541.812.10.454.815.2175.6452.414.60.404.29.81121.6861.614.10.503.914.5178.3471.813.60.403.913.5485.0181.813.10.354.814.0581.8991.812.60.353.914.6178.76101.813.10.504.214.0581.89111.612.60.504.216.2370.01122.012.10.504.213.8784.04132.413.10.404.210.94109.18142.614.10.254.29.49127.31151.813.10.504.514.0681.89162.414.60.403.99.81121.68171.812.60.354.214.6178.76182.012.10.255.113.8684.04192.012.60.305.113.3278.75201.613.10.303.915.6172.79In order to verify the correctness of the optimization results, taking the groups 2, 9, 17, and 20 in Table4 as an example, the maximum temperatures of windings are 54.91°C, 67.74°C, 68.56°C, and 73.40°C, respectively, based on the thermal network model. Meanwhile, the simulation results based on the finite element method are given in Figure 13; the maximum temperatures of windings are 55.10°C, 68.20°C, 70.00°C, and 72.63°C respectively; it can be seen that the temperature distribution laws of the four groups are basically the same, and the error coefficient is less than 3%; thus, the correctness of the optimization method has been verified.Figure 13
Simulation results: (a) Case 2; (b) Case 9; (c) Case 17; (d) Case 20.
(a)(b)(c)(d)In Table4, it can be seen that the metal conductor usage of windings is reduced by 23.05%, and the loss of the windings is reduced by 20.25% compared with the initial design parameters. As the purpose of low metal conductor usage and loss of winding are conflicted and cannot be optimized simultaneously, therefore, in the actual design process of the iron reactor, it is necessary to select the suitable structural parameter to realize the optimization of one objective with the other objectives which are not deteriorating.


## 6. Conclusions

According to the structure and parameters of oil-immersed iron core reactor, a multiobjective optimization method based on particle swarm optimization and thermal network model is proposed, and the following conclusions can be obtained:(1)
The maximum temperature error of discs is less than 3% between the thermal network model and finite element method, and the position of the temperature is basically the same; the correctness of thermal network model has been verified(2)
The maximum temperature region is mainly concentration of the middle of windings; the reasons are given by analyzing the fluid velocity distribution of oil, so this area is the main focus to monitoring temperature for the windings(3)
The optimization design of windings based on particle swarm algorithm and thermal network model is proposed, the metal conductor usage is decreased by 23.05%, and the loss of the windings is reduced by 20.25% compared with the initial design parameters, so the optimization method has an important guiding significance for the design of oil-immersed iron core



---

*Source: 6642620-2021-02-17.xml*