Calculation of linear induction motor features by detailed equivalent circuit method taking into account non-linear electromagnetic and thermal properties

Computers and Mathematics with Applications 78 (2019) 3187–3199

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Calculation of linear induction motor features by detailed equivalent circuit method taking into account non-linear electromagnetic and thermal properties ∗

Ivan Smolyanov , Fedor Sarapulov, Fedor Tarasov Ural Federal University, Russian Federation

article

info

Article history: Available online 11 June 2019 Keywords: Detailed equivalent circuit CFD Thermal modes Linear induction motor Non-linear aerodynamic resistance Analytical coefficient

a b s t r a c t The main purpose of the work is to suggest the technique of solving multi-physics problems (calculation of magnetic, temperature and velocity fields) for the purpose of units with forced cooling system and complex shapes. The general problem is the presence of narrow region in which it is necessary to calculate the velocity field of a cooling liquid. The problem is solved by a simplified method based on a detailed equivalent circuit, which significantly reduces required computational resources in comparison with numerical methods (finite element and finite volume methods). The paper presents the transition from partial differential equations to algebraic equations, created on the basis of detailed equivalent circuits. The results of the model are analyzed on a real linear induction motor applied in the transport system. Thermal modes are being considered depending on cooling intensity and magnitude of power supply parameters of the motor. It is proposed to reduce the computation time by simplified calculation of the velocity field using an analytical coefficient obtained from the analysis of the equivalent circuits. The vectorization of cycles intended for building the stiffness matrix, storage of the matrix in sparse form, and rational sequences of matrix calculations in the course of the solution made it possible to reduce the computational time almost twice as well. The comparison of the accuracy of the results obtained is presented, depending on the discretization of the computational domain and chosen type of interpolation. Verification results of electromagnetic and thermal calculations obtained from experimental data are also presented. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Nowadays, there are many fields starting from high-speed surface transport to various automation and electromechanical systems where motors with open core are used. Analysis of thermal processes in them is a vital task in studies of such electric drives with linear induction motor. Systemic research of this field began in the 40-50 s of the last century [1–3]. And the first notes about devices operating with the linear inductor can be found in works from 1920 [4,5]. Although the research of linear motors is not new, a number of problems in this domain still remains unresolved. For example, there are problems associated with simulation of superconductivity phenomena for high-speed linear motors [6], building of simple and accurate models for optimization problems and parametric studies [7], the calculation of specific structures ∗ Corresponding author. E-mail address: [email protected] (I. Smolyanov). https://doi.org/10.1016/j.camwa.2019.05.015 0898-1221/© 2019 Elsevier Ltd. All rights reserved.

3188

I. Smolyanov, F. Sarapulov and F. Tarasov / Computers and Mathematics with Applications 78 (2019) 3187–3199

and a number of others. Moreover, the results can be used in related scientific fields, for instance, pumping [5,8] or mixing [9] liquid metals. One of such problems is analyzing thermal operating conditions of linear motor. The necessity of such analysis appears both at the stage of the motor design and operational stage (changing the operational mode, changing the other parts of electric drive etc.) [10–12]. For instance, very strict requirements are imposed on geometrical dimensions of traction motors in the design of transport systems. It may be proved that motor does not generate required driving force over the required time due to structural components (usually coil windings) overheating that can cause their damage. Research of transient thermal processes should be performed when any solution to the problem (usage of cooling plants, temperature control systems, changing the engine structure etc.) is applied. A great number of narrow (thin) areas in the structure represent the main problem of thermal processes calculation for mechanisms of such type. It requires a very fine model discretization when solving the problem by numerical methods, for instance, finite element method or finite volume method [13,14]. As a consequence, the computational cost of such models increases significantly. At the present time, most of works on calculation of multi-physics problems (magnetic and temperature fields) are based on finite element method and non-detailed equivalent circuit. The authors [15,16] suggest a model based on a low-detailed circuit. The comparisons show good agreement between the results obtained by the detailed equivalent circuit and by the 3D FEM model, but the results obtained by the low-detailed equivalent circuit method require accurate determination of a set of analytical coefficients and functions. In the works [17,18], the problems of optimized design with regard to the temperature conditions are solved. A bunch of temperature and magnetic fields is analyzed using conventional equivalent circuits, which can significantly reduce the time of solving the optimization problems. Analysis of linear motors with cooling system using an analytical model was proposed by the authors [19]. The calculation of temperature fields for simple constructions based on the 2D finite element method was carried out in [20] and on the basis of 3D FEM (without taking into account the convection) in [21]. All these approaches have a number of advantages but are poorly focused on the calculation of installations with narrow simulated domains. In the article, an application of a simplified technique of numerical simulation based on detailed equivalent circuits used for calculation of velocity field of cooling medium, thermal processes and electromagnetic field induced by the stator is considered in solution of the coupled problem. The approach makes it possible to avoid building a mesh due to describing the behavior of physical phenomena based on the circuit theory. The article serves as an extension of a series of studies [22–24] where an essential simplification of calculations without loss of results accuracy is proposed. 2. Mathematical model description A model of coupled electromechanical, thermal and aerodynamic processes is implemented in Matlab software by detailed equivalent circuit method. The main model of the operation algorithm is depicted in Fig. 1. Automated control system unit generates a power supply signal (current dependence of voltage U(I) and frequency supply f ) by means of task for velocity for electromagnetic field calculation in Detailed magnetic equivalent circuit unit. Velocity (Vmov ) of the moving element and heat power (P) in different motor parts is calculated taking into account load forces (Ff ) in Detailed magnetic equivalent circuit unit. More detailed mathematical description of detailed magnetic equivalent circuit unit is given in Section 2.1 section and in papers [13,14,24–26]. Thermal processes are calculated in the detailed aerodynamic equivalent circuit unit (Section 2.2) by means of detailed thermal equivalent circuit taking into account cooling air movement in the area of inductor. The value of input cooling air velocity (Vc ,i ) is set at the input of the unit. Air velocity distribution in each branch is a solution. It should be noted that the magnetic and electrical properties depend on temperature (T ). 2.1. Detailed magnetic equivalent circuit The transition from partial differential equations to magnetic circuit calculation can be performed with the application of circuit analysis [24,26] and finite-difference approximation of derivatives [27,28]. On assumption that plane-parallel electromagnetic problem is being solved let us write the differential equation for phasor A:

∂ 2A ∂ 2A ∂A + 2 = ȷ · ωγ A − µJs + µγ Vx , ∂ x2 ∂y ∂x

(1)

where ω is the magnetic field angular frequency, γ denotes the electric conductivity of medium, µ stands for the magnetic permeability of medium and Js represents the extraneous current density, Vx is the moving velocity along x- axis. If the domain of magnetic field of thickness L is divided into right angled grid cells (Fig. 2(b), a) with a pitch of L1, L2, L3, L4, magnetic vector potential can be expanded in a Taylor series or well-known finite-difference approximation of derivatives with respect to coordinates can be used. Then, provided that L1 = L2 = L3 = L4 = h we obtain (2) [27,28] Ai+1,j − 2Ai,j + Ai−1,j

−

h2 ȷ · ωγ Ai,j

µ·h

+

Ai,j+1 − 2Ai,j + Ai,j−1 h2

− γ Vx Φi+1,j + Js = 0

.

(2)

I. Smolyanov, F. Sarapulov and F. Tarasov / Computers and Mathematics with Applications 78 (2019) 3187–3199

3189

Fig. 1. Flow chart of mathematical model.

By Stokes theorem for node 0 node (as shown in 2(b)) magnetic flows will be written as follows:

Φi+1,j = (Ai+1,j − Ai,j ) · L; Φi−1,j = (Ai−1,j − Ai,j ) · L Φi,j+1 = (Ai,j+1 − Ai,j ) · L; Φi,j−1 = (Ai,j−1 − Ai,j ) · L ,

(3)

and then Eq. (2) in view of (3) can be written as follows:

Φi+1,j + Φi−1,j Φi,j+1 + Φi,j+1 ωγ γ Vx Φi+1,j + −ȷ Φi,j + Js − =0 2h2 µL 2h2 µL µL L

(4)

or

Φi+1,j

(

1 2 µL 2h(

+Φi,j−1

+

γ Vx

1

2h2 µL

)L

)

+ Φi−1,j

− Φi,j

ȷωγ

µL

(

1

)

2h2 µL

(5)

= Fi,j

and in simplified form:

Φi+1,j Ri+1,j + Φi−1,j Ri−1,j + Φi,j+1 Ri,j+1 + Φi,j−1 Ri,j−1 − Φi,j Ri,j = Fi,j .

(6)

It follows that computational domain shown in Fig. 2(a) can be presented in the form of magnetic equivalent circuit (Fig. 2(b)). Thus, magnetic fluxes in branches of magnetic circuit can be calculated by Kirchhoff’s laws for magnetic circuits. Transformation using transition to magnetic circuits allows move away from solving partial differential equations, turning to solving ordinary differential equations or linear algebraic equations a system depending on formulation of the problem. With this approach, it is possible to significantly reduce the number of unknown unknowns by compiling the system of equations on methods for calculating circuits (loop current method, nodal potential method, etc.). In planeparallel problems based on the detailed equivalent circuit, it is possible to take into account the third coordinate only for a part of the design under study with inclusion of additional resistances. Therefore, one can create both three-dimensional model and their hybrids (mixed two- and three-dimensional tasks) using detailed magnetic equivalent circuit, it makes possible to reduce the number of unknowns. In the article electromagnetic processes are calculated by means of magnetic circuits (Fig. 3), where RtQ ,N , RnQ ,N - the magnetic resistance of tangential (along the inductor) and normal component of magnetic flux above the Q-slot in the Nth layer and FQ is a magneto motive force of Q-slot. The circuit can be expressed mathematically in terms of the following matrix equation: Z M Φ = Ic Is , M

(7) c

where Z is the matrix of magnetic impedance, Φ is the magnetic flux matrix, I is the column vector of currents in layers of computational domain, Is = KI IΦ is the column vector of stator currents, KI is the linear transformation matrix of phase currents into slot currents, IΦ is the column vector of phase currents.

3190

I. Smolyanov, F. Sarapulov and F. Tarasov / Computers and Mathematics with Applications 78 (2019) 3187–3199

Fig. 2. Arbitrary computational domain (a), equivalent circuit for arbitrary computational domain (b).

Fig. 3. Detailed magnetic equivalent circuit.

2.2. Aerodynamic equivalent circuit The flow continuity equation and the equation describing velocity distribution of fluid (gas) for the steady-state case

∇ · ρu = 0 (u · ∇ )u = f −

1

ρ

1

∇ p + v ∆ u + v ∆ u. 3

(8)

Here, ρ denotes the mass density, ν denotes the kinematic viscosity, f stands the vector of tension of mass forces, p is the pressure. The cooling flow can be considered incompressible. In the simplest case, for one-dimensional incompressible flow, the system will take the form

ρ

dQ

=0 dl dQ dp ρu = ρ fL − . dl dl

(9)

I. Smolyanov, F. Sarapulov and F. Tarasov / Computers and Mathematics with Applications 78 (2019) 3187–3199

3191

Considering a motor as a hydraulic system, they limit to two first equations, whose integrals have the form and u2 /2 + U + p/ρ = const (U is depicted potential of mass forces). When a magnitude of viscous medium flows in a channel with a variable cross-sectional area F these relations (9) can be reduced to the Bernoulli equation (10).

ρ uF = const pi + ρ u2i /2 + ρ Ui = pi+1 + ρ u2i+1 /2 + ρ Ui+1 + ∆p

(10)

where ∆p means the pressure loss between points i and (i + 1) due to friction. According to Newton’s hypothesis about the similarity of hydromechanical flows, the friction can approximately be expressed by the formula F = kρ Su2 . If one assumes that the dependence of resistance for all influencing factors should be taken into account by the coefficient k then the equation for pressure drop caused by resistance is

∆p =

F S

= kρ u2

(11)

In order to express the pressure drop ∆p in terms of the dynamic pressure ρ u2i /2, we introduce a coefficient ∆p = ξ ρ u2 /2 then ξ will be equal to 2k and we will call the drag coefficient. Observing the similarity criterion for hydro mechanics, it is possible to write the pressure drop for the same geometric objects as

∆p = ξ (Re)

ρ

u2 (12) 2 Often it makes sense to introduce a total coefficient of local resistances of the non-rectifying channel that is defined as the sum of the coefficients of resistances of thorns, turns and narrowings (or expansions). In this paper, the problem is calculated with respect to the flow rate Q , to which it is easy to proceed from the principle of conservation of mass (13) as applied to a channel with a different cross section S along the length. u1 S1 = u1 S1 = · · · = un Sn = Q = const

(13)

Calculation of aerodynamic ventilation reduces to determining the cooling medium velocities in channels, cooling system hydraulic resistance and pressure losses in the channel from input of cooling agent into the machine to the output. Differences in heating of the cooling agent in different parts of the cooling system is not taken into consideration and set as equal for all parts. This assumption is usually reasonable. Moreover, consideration of differences in heating of cooling agent requires additional iterative process that leads to a reduction of model performance and impossibility of its application as a part of the full dynamic thermal model. Let us compose a simplified model [22,29,30] (Fig. 4(a)) and determine colors for the following parts: inductor teeth (gray color), inductor slots (yellow color) and ends winding (orange color). The cooling air (or special gas) circulates in slot and end winding only. That is why the aerodynamic resistances and flow rates are calculated just in these areas of npk trp ep computational domain. The following notations Zout , Zin , Znkol , Zntro , Zn , Zn , Zn used in Fig. 4 corresponds to aerodynamic resistances of input, output, knee, T-joint in a branch, slot, T-joint in a direct channel and end winding for nth section respectively [22]. Transform now the equivalent circuit to an easy-to-calculate view (Fig. 4(b)), then the aerodynamic resistances can be written (depending on flow rate) as follows [22]: Znd Znc

Qns

(

=

p

Qn−1 Qns

(

Qnc−1

( Zns

)

Qns Qnd−1

)

Zntrp

= Zntrp )

, Resn

Qns

(

)

p

Qn−1 Qns

(

Qnc−1

= 2 · Zntro

+ Znep

)

+ Znep ) ( Qns

Qnd−1

(14)

( ) + Znep Resn ,

where Zns , Znc , Znd are the aerodynamic resistances of the nth section of channel, the aerodynamic resistances of collecting and the distributing branches respectively; Qns , Qnd , Qnc is the flow rate in the nth section of channel, the flow rate of collecting and the distributing branches respectively. Reynold’s number dependence of track resistance will be neglected. Eq. (15) are written in accordance with Kirchhoff’s first law for distributing branch, number of equations being (N − 1), i.e. n = 1 . . . N − 1. Qns + Qnd + Qnc = 0.

(15)

Qnd−1 ,

when n = N: there is no equation as there is no zero node. When n = 1 : Equations by Kirchhoff’s first law for collecting branch, number of equations being (N − 2), i.e. n = 1 . . . N − 2: Qnc + Qns + Qnc−1 = 0.

(16)

When n = 1: there is no equation as there is no node, when n = 2: , when n = N: we do not write an equation. Equations by Kirchhoff’s second law, number of equations is (N − 1) i.e. n = 1 . . . N − 1: Zns · Qns

(

)2

( )2 ( )2 ( )2 + Znc · Qnc − Znk+1 · Qns+1 − Znd · Qnd = 0.

(17)

3192

I. Smolyanov, F. Sarapulov and F. Tarasov / Computers and Mathematics with Applications 78 (2019) 3187–3199

Fig. 4. General aerodynamic equivalent circuit of cooling system (a), transformed detailed aerodynamic equivalent circuit (b).. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

It should be noted that in this equation aerodynamic resistances are flow rate-dependent variables according to Eqs. (15). When the flow rate is determined and the section is known, we can calculate a medium velocity for each section. Solution procedure of simultaneous equations is as follows: • 1st stage. Calculate the flow rate in branches provided that local resistances do not depend on flow rates ratio. Resistance and flow rate values obtained will be used as initial values. • 2nd stage. Calculate T-joint resistances depending on flow rates obtained in the 1 stage. • 3rd stage. Solve simultaneous equations to determine flow rates and velocities in branches. Flow rates obtained in the previous iteration will be set as initial values for solving the simultaneous equations, resistances calculated in the previous iteration will be used as coefficients. • 4th stage. Stages 2 and 3 should be repeated until calculation results in two subsequent iterations agree with each other enough.

The calculation of velocity of fluid based on aerodynamic resistance in most cases is applicable to a pipe-shaped and similar structures. But with a correct choice of aerodynamic resistance based on the similarity criterion for other similarity groups, this method can be generalized. 2.3. Detailed thermal equivalent circuits Thermal conductivity equation in case of mass transfer along x axis will be written as follows [31]:

∂T ∂ 2T Q ∂T +V =a 2 + , ∂t ∂x ∂x cp d

(18)

where V = dx is the moving velocity; p is the specific material heat capacity; d is the material density; Q is the dt total power of heat release; a is the thermal diffusivity coefficient, T is the temperature. Using the finite-difference approximation, Eq. (12) will be written as [32–34]: cpi · Mi

dTi dt

=

λi Si hi

(Ti+1 − Ti ) +

+cpi · Mi

Vx hi

λi Si hi

(Ti−1 − Ti )

(Ti−1 − Ti ) + Pi ,

(19)

where cpi denotes the heat capacity, Mi is the mass, Ti is the temperature, lambdai stands for the heat conductivity, Si is the surface area, hi represents the length, Pi is the power for the ith section. By analogy with 2.1 and 2.2 sections let us compose detailed thermal equivalent circuit (Fig. 5) taking into account all thermal fluxes. According to this circuit,

I. Smolyanov, F. Sarapulov and F. Tarasov / Computers and Mathematics with Applications 78 (2019) 3187–3199

3193

Fig. 5. Detailed thermal equivalent circuits for frontal (a), and profile place (b).

the system of ordinary differential equations is compiled, which is solved using the explicit multi-step Adams–Bashforth method with a variable order. 3. Results Let us consider possibilities of detailed equivalent circuits method based on linear induction motor calculation for monorail transport put into operation in Moscow. Strict requirements were imposed on width of inductor in the design of traction motor. In this context winding diagram where all 50 slots are half-filled was taken. The other half of each slot section is the channel for cooling air. Basic geometrical dimensions are shown in Fig. 6. The motor operates within a voltage range from 40 to 240 V, a frequency range from 5 to 40 Hz, phase current range from 70–270 A with max generated force of 7500 N. Inductor’s width is equal to 240 mm, secondary element width – 400 mm. The coils and the conductive layer are made from cooper, the core and the back iron consist of laminated iron and entire part of space is filled with cooling liquid. The kinematic viscosity of the coolant obeys the law (20), and the thermal conductivity depending on the average temperature (21). The heat capacity of the fluid is constant and equal to 3000 J/(kg K), and the mass density is 1100 kg/m3 . The dependences of electrical and thermal conductivity of the copper is described by the analytical expression (22). The heat capacity is 430 J/(kg K). The thermal conductivity and the heat capacity of steel are set at 20 W/(kg K) and 500 J/(kg K) and the saturation curve taken from [35]. Let us consider possibilities of detailed equivalent circuits method based on calculation of the linear induction motor for monorail transport put into operation in Moscow. Strict requirements were imposed on width of the inductor in the design of traction motor. In this context, the winding diagram where all 50 slots are half-filled, was considered. The other half of each slot section is the channel for cooling air. Basic geometrical dimensions are shown in Fig. 6. The motor operates within the voltage range from 40 to 240 V, the frequency range from 5 to 40 Hz, phase current range from 70 to 270 A with maximal generated force of 7500 N. The inductor’s width is equal to 240 mm, the secondary element width is 400 mm. The coils and conductive layer are made of cooper, the core and back iron consist of laminated iron and entire part of space is filled with cooling liquid. The kinematic viscosity of the coolant obeys the law (20), and thermal conductivity depends on the average temperature (21). The heat capacity of the fluid is constant and equal to 3000 J/(kg K) and mass density is 1100 kg/m3 . The dependencies of electrical and thermal conductivities of copper are described by analytical expressions (22). The heat capacity is 430 J/(kg K). The thermal conductivity and heat capacity of steel are set at 20 W/(kg K) and 500 J/(kg K), respectively, and the saturation curve is taken from [33].

νm (Tm ) = 7, 33 · 10−11 · Tm2 + 9, 737 · 10−8 · Tm + 1, 356 · 10−5 ,

(20)

λm (Tm ) = 5, 586 · 10−5 · Tm + 0, 025,

(21)

λCu (T ) = 393, 22 − 0, 0712 · T [ ( )]−1 σ = 1.667 · 10−8 1 + 3.862 · 10−3 (T − T0 ) ,

(22)

3194

I. Smolyanov, F. Sarapulov and F. Tarasov / Computers and Mathematics with Applications 78 (2019) 3187–3199

Fig. 6. Mathematical model geometry.

Fig. 7. Cooling agent velocity in distributing (a) and collecting (b) branches and in slots.

Let us perform calculations for cooling medium velocity equal to 10 m/s at the entrance to the system. The velocity curves (Fig. 7) will be calculated for the first and the last iterations in order to demonstrate the influence of flow rate ratios on their aerodynamic resistance and final result respectively. Flow rate curves and pressure drops curves will be shown for the last iteration only. An error is too high particularly for some regions in calculating with constant resistances. A reasonable number of iterations must be at least 10 (error is less than 1%) or 5–6 (error is less than 5%–7%). Calculation with 15–20 iterations gives an adequate accuracy. Influence of flow rate dependence of aerodynamic resistance on each section can be estimated quantitatively by comparison of total pressure drop and aerodynamic resistance in the system. The calculation has shown that with due regard to flow rates dependence, the total pressure drop is 3858.8 Pa and resistance is equal to 79 727 Pa s/m6 and 4129.9 Pa and 85 329 Pa s/m6 without due regard to its linear dependence, respectively. Resistances of channels with low flow rate (middle section of inductor) increase in the course of iteration processes and resistances of channels with high flow rate decrease. It leads to velocity (and flow rates) redistribution as seen in Figs. 7(b) and 8(a). Almost full absence of ventilation in several central channels is observed that causes deterioration of winding cooling conditions in this section. Dramatic dips to zero point shown in Figs.7b and 8 correspond to fully filled slots (2 and 49). It should be noted that some amount of cooling air enters these slots as seen in Fig. 8(a) even if velocity in the filled channels is very small. Significant velocity difference of cooling medium in collecting and distributing branches for the first and the last sections is observed. This fact corresponds to physical content of processes in the plant that are connected with the fact that cooling air enters in the distributing branch and leaves through the collecting branch. Calculation of velocity distribution takes 70% of total computation time necessary for solving the problem. Computation time can be reduced significantly when introducing the analytical coefficients based on the model. The velocity in each circuit section can be determined by means of these coefficients according to formulas Vns (Vin ) = Ks (n) · Vin , Vnc (Vin ) = Kc (n) · Vin , Vnd (Vin ) = Kd (n) · Vin .

(23)

I. Smolyanov, F. Sarapulov and F. Tarasov / Computers and Mathematics with Applications 78 (2019) 3187–3199

3195

Fig. 8. Flow rate of cooling liquid and gas (a) and pressure drop in different circuit sections (b).

where Vns , Vnc , Vnd are the velocity in inductor’s slots, collecting and distributing branches of the nth section, respectively, Vin is velocity at the entrance to cooling system, Ks , Kc , Kd are the analytical coefficients depending on the order number of branch for slot, collecting and distributing branches, respectively. Calculation results for various air velocities at the entrance (from 1 to 10 m/s) are shown. When the velocity value is far less than 1 m/s Reynold’s numbers in all sections will be sufficiently small i.e. natural convection equations will be used for calculation. Let us decrease number of zones from 50 to 13, put it differently we reduce model’s discretization where the first zone contains the first three slots and 13th zone represents the last three slots, other zones are enlarged sections containing 4 slots. The results of model calculation are shown in Fig. 9. If we compare these results (Fig. 9) with velocity distribution for input velocity of 10 m/s (Fig. 7) we can conclude that reduction of model discretization influences on results accuracy insignificantly. Calculation results for various air velocities at the entrance (from 1 to 10 m/s) are shown. When the velocity value is much less than 1 m/s, Reynold’s numbers in all sections will be sufficiently small, i.e., natural convection equations will be used for calculation. Let us now decrease the number of zones from 50 to 13 putting them differently. In this way we reduce the model’s discretization where the first zone contains the first three slots and 13th zone represents the last three slots, other zones are enlarged sections containing 4 slots. The results of model calculation are shown in Fig. 9. If we compare these results (Fig. 9) with velocity distribution for input velocity of 10 m/s (Fig. 7) we can conclude that the reduction of model discretization influences the accuracy of results only insignificantly. Based on dependencies (Fig. 9), analytical functions (Fig. 10) that allow determining the velocity in any section of the cooling system by formulas (23) (if the input velocity is a known value) were derived. The fact makes it possible to exclude aerodynamic processes calculation each time when the temperature field is calculated. Linear induction motor running at frequencies of 5, 10, 15, 23.2 and 40 Hz4 was considered. Certain currents were given for each frequency: 135, 165, 190, 215, 240 and 270 A, respectively. The idle mode under natural ventilation conditions (cooler velocity at the entrance to the system is equal to train speed) was simulated. Based on inductor length-averaged temperature values of winding (Fig. 11) we can say that steady-state temperature value is higher than 180–200 ◦ C when current value is more than 135 A. It should be noted that in steady-state thermal process the maximum temperature of winding is 1.3–1.5 times higher than inductor length-averaged temperature of winding shown in the curve (Fig. 11). In transition thermal processes, the temperature of winding can be 1.3–1.8 times higher than inductor length-averaged temperature (due to the velocity of temperature rise in different construction elements it varies). Heating up of the secondary element was not considered in the article because it is insignificant in comparison with heating up of inductor winding. Current value of 135 A is within the range of operating currents of the machine. That is why the induced ventilation (cooling air or liquid) is required. Let us consider the most difficult mode in the context of thermal behavior with a frequency of 5 Hz and load of 8000 N. The model allows us constructing the dependence of temperature on flow rate of input cooling agent (Fig. 12). As it can be seen in the curve, an adequate flow rate of cooling agent ensures a safe mode of linear motor inductor (the winding temperature is lower than 200 ◦ C). Results given Fig. 12 are in agreement with experimental data provided by TEMP Engineering and Scientific center. Besides, electromagnetic calculations were verified by comparison of experimental data with simulation data by forces (Table 1). The most expensive mathematical operation is to solve a system of algebraic equations. The study of various methods for solving systems of linear algebraic equations in that work was not carried out. The compilation of the stiffness matrix

3196

I. Smolyanov, F. Sarapulov and F. Tarasov / Computers and Mathematics with Applications 78 (2019) 3187–3199

Fig. 9. Velocity distribution in slots (a), collecting branch (b), distributing branch (c).

Fig. 10. Analytical functions for equations.

I. Smolyanov, F. Sarapulov and F. Tarasov / Computers and Mathematics with Applications 78 (2019) 3187–3199

3197

Fig. 11. Dependencies of inductor length-averaged temperature of winding at various frequencies.

Fig. 12. Dependence of inductor length-averaged temperature of winding on flow rate of cooling medium. Table 1 Comparison of experimental data with calculation results. Phase current, A

135

165

190

215

240

270

Ft Ft Ft Ft Ft Ft

1950 2095 1650 1780 1550 1421

2750 3129 2150 2659 2200 2122

3800 4149 2900 3426 2950 2814

5200 5313 4100 4515 4050 3604

6750 6621 5300 5626 5200 4490

8400 8379 6800 7121 – 5883

(exp. 5 Hz), N (theor. 5 Hz), N (exp. 10 Hz), N (theor. 10 Hz), N (exp. 15 Hz), N (theor. 15 Hz), N

was based on the concept of the theory of circuits depending on the number of nodes and branches of the circuit. The next most expensive computer time is the function of filling the matrices of the model. The vectorization of cycles, storage of matrices in a sparse form, and rational sequence of matrix calculation at the time of solution made it possible to reduce the computation time by almost two times. It is also possible to reduce the calculation time using different types of interpolation and different model of discretization step. Therefore, it should be estimated the error of calculations of the arrays of thrust, core losses and secondary element (SE) losses for different type of interpolation (linear and spline) and the number of degrees of freedom (DOF) (Table 2). The mean square relative error of interpolation over all points and the relative error of interpolation at each point (in order to identify the maximum or exceeding a specified level) were considered. Test points are selected randomly. The root-mean-square error of piecewise-linear (linear) interpolation for all arrays of characteristics does not

3198

I. Smolyanov, F. Sarapulov and F. Tarasov / Computers and Mathematics with Applications 78 (2019) 3187–3199 Table 2 Comparison of relative error using different type of interpolation and degrees of freedom. DOF

1 056 5 670 111 670

Thrust

Core losses

SE losses

Linear

Spline

Linear

Spline

Linear

Spline

0.0311 0.007 0.00255

0.0155 0.0027 0.0018

0.0625 0.0101 0.00275

0.0358 0.00435 0.0026

0.0822 0.0239 0.0113

0.03 0.00215 0.00084

exceed 8.25%, even in the case of the smallest DOF array. The non-convergence is much smaller for spline interpolation with a high number of DOFs and does not exceed 2.5%. 4. Conclusion The possibility of calculating magnetic, temperature and velocity fields using simplified method based on detailed equivalent circuit have been presented. The approach to solving significantly simplifies the procedure of transition from the physical laws to a system of algebraic or ordinary differential equations due to the specific nature of the dimensional discretization using equivalent resistances. The opportunity and correctness of representation of physical phenomena with help of detailed equivalent circuit have been shown using the derivation of equations compiled according the equivalent circuits from classical partial differential equations describing the physical laws. The results were analyzed and verified on a real linear motor for the purpose of a transport system. The question of accounting for nonlinearities in models was considered qualitatively and quantitatively. The need to take into account nonlinearity when calculating the velocity field especially depends on the values of the input velocity. One of the key elements discussed in the paper is an opportunity to reduce the calculation time for the model. The effect of the results accuracy depending on the degrees of freedom and type of interpolation has been considered quantitatively with the help of the root-mean-square relative error. The required discretization of the computational domain and the choice of the type of interpolation significantly depends on the degree of nonlinearity of the calculated functions. Using the analytical coefficient to determine the coolant velocity in various parts of the design can significantly reduce the computation time of the entire model. The analytical function can be calculated for any other device that is calculated using this model. The proposed approach to analyzing the velocity field is advisable to use for systems with narrow regions where it is difficult to build a finite element mesh (due to the high disproportion or complexity) or to solve the problem in this area (high non-convergence). The aerodynamic resistances in this paper are calculated based on the similarity criterion for tube-type structures. Therefore, the model will have high accuracy for problems with pipe-shaped geometry, but it should be noted that it is possible to expand the use of this method for any other design, taking into account the correct choice of the aerodynamic coefficient of friction for the shape of the design. Models based on detailed equivalent circuits can be considered as hybrid spatial problems, i.e., they consider part of the task as two-dimensional, and part as a three-dimensional formulation. In the following papers, the possibility of reducing the calculation time of such models using this representation will be considered. Acknowledgment The work was supported by Act 211 Government of the Russian Federation, contract no. 02.A03.21.0006. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]

E.R. Laithwaite, V. Duxbury, Electromagnetic shuttle-propelling devices, J. Textile Inst. Proc. 48 (1957) P214–P224. R.E. Graham, Linear servo theory, Bell Syst. Tech. J. 25 (1946) 616–651. Gradenwitz, Le Canon Elecromagnetique de Birkeland. Eclairage 26, 267. J. Hartmann, Apparat zur erzeugung eines ununterbrochenen liitenden, in: Flussigkeitsstrahles Mit Regelbarer Geschwindigkeit, 1920. V. Geze, B. Nacke, Numerical simulation of core-free design of a large electromagnetic pump with double stator, Magnetohydrodynamic 52 (2016) 417–431. L. Jin, Z. Deng, W. Lei, H. Li, J. Li, N. Qian, Dynamic characteristics of the HTS maglev vehicle running under a low-pressure environment, IEEE Trans. Appl. Supercond. 29 (2019) 1–4. A. Shiri, A. Shoulaie, Design optimization and analysis of single-sided linear induction motor, considering all phenomena, IEEE Trans. Energy Convers. 27 (2012) 516–525. F. Sarapulov, S. Sarapulov, V. Frizen, Use of detailed equivalent circuitmethod for investigation of electromagnetic, thermal and hydrodynamic processes in induction electric engineering units, Acta Tech. CSAV 60 (2015) 131–153. E. Shvydkiy, V. Zaharov, K. Bolotin, I. Molyanov, S. Sarapulov, Numerical modeling of the travelling magnetic field stirrer for liquid lithium, Magnetohydrodynamics 53 (2017) 707–713. J. Gieras, Linear Induction Drives, Oxford University Press, New York, 1994, p. 293. R. Hellinger, P. Mnich, Linear motor-powered transportation: History, present status, and future outlook, Proc. IEEE (2009) 1892–1900.

I. Smolyanov, F. Sarapulov and F. Tarasov / Computers and Mathematics with Applications 78 (2019) 3187–3199

3199

[12] J. Gieras, J. F. Mews, P. Splawski, Analytical Calculation of electrodynamic levitation forces in a special-purpose linear induction motor, IEEE Trans. Ind. Appl. 48 (2012) 106–116. [13] F. Sarapulov, I. S.M.olyanov, F. Tarasov, Study of the linear induction motor with bimetallic secondary element, Acta Tech. CSAV 63 (2018) 205–220. [14] F. Sarapulov, S. Sarapulov, I. Smolyanov, Research of thermal regimes of linear induction motor, in: 2017 18th International Conference on Computational Problems of Electrical Engineering (CPEE), 2017, pp. 1–4. [15] T. Tong, J. Gong, Y. Gao, N. Bracikowski, Lecture Notes in Electrical Engineering, 2018, pp. 563–570. [16] N. Bracikowski, M. Hecquet, P. Brochet, S.V. Shirinskii, Multiphysics modeling of a permanent magnet synchronous machine by using lumped models, IEEE Trans. Ind. Electron. 59 (2012) 2426–2437. [17] D. Pan, L. Li, M. Wang, Modeling and optimization of air-core monopole linear motor based on multiphysical fields, IEEE Trans. Ind. Electron. 65 (2018) 9814–9824. [18] S.R. Aleksandrov, T.T. Overboom, E.A. Lomonova, Design optimization and performance comparison of two linear motor topologies with PM-less tracks, IEEE Trans. Magn. 54 (2018) 1–8. [19] L. Zhang, B. Kou, Y. Chen, Y. Jin, Analysis and calculation of eddy current braking force for an ironless linear synchronous motor with cooling system, in: 9th International Conference on Electrical Machines and Systems, 2016, pp. 1–5. [20] Y.H. Wu, X.M. Liu, C.F. Ji, Research on the temperature field of special linear motor, Appl. Mech. Mater. 416–417 (2013) 169–174. [21] J.-J. Kim, Y.H. Jeong, D.-W. Cho, Thermal behavior of a machine tool equipped with linear motors, Int. J. Mach. Tools Manuf. 44 (2004) 749–758. [22] F.N. Sarapulov, V.V. Goman, Development of mathematical models of thermal processes in linear asynchronous motors, Russ. Electr. Eng. 80 (2009) 431–435. [23] F. Sarapulov, V. N. Frizen, E. E. Shvydkiy, I. L. Smol’yanov, A, Mathematical modeling of a linear-induction motor based on detailed equivalent circuits, Russ. Electr. Eng. 89 (2018) 270–274. [24] V. Dmitrievskii, V. Goman, F. Sarapulov, V. Prakht, S. Sarapulov, Choice of a numerical differentiation formula in detailed equivalent circuits models of linear induction motors, in: 2016 International Symposium on Power Electronics, Electrical Drives, Automation and Motion (SPEEDAM), 2016, pp. 458–463. [25] A. Nosrati, J. Nazarzadeh, Analysis of linear induction machines with internal fault by MEC, COMPEL - Int. J. Comput. Math. Electr. Electron. Eng. 36 (2017) COMPEL–08–2016–0360. [26] V. Begalov, F. Sarapulov, A. Gorlov, Propulsion and normal forces in linear induction drives, in: Proc. of ACED–04, 2004. [27] T. Wu, R. Xu, An optimal compact sixth-order finite difference scheme for the Helmholtz equation, Comput. Math. Appl. 75 (2018) 2520–2537. [28] A. Arakelyan, Convergence of the finite difference scheme for a general class of the spatial segregation of reaction–diffusion systems, Comput. Math. Appl. 75 (2018) 4232–4240. [29] A. Huminic, G. Huminic, A. Soica, Study of aerodynamics for a simplified car model with the underbody shaped as a venturi nozzle, Int. J. Veh. Des. 58 (2012) 15. [30] S. Arsene, I. Sebesan, G. Popa, V. Stefan, The aerodynamic resistances determined by the rolling equipment for the electric locomotive LE 060 EA 5100 KW. 2015. [31] M. Goffredo, M. Schmid, S. Conforto, F. Bilotti, C. Palma, L. Vegni, D. Alessio, M. Goffredo, M. Schmid, S. Conforto, F. Bilotti, C. Palma, L. Vegni, COMPEL : The International Journal for Computation and Mathematics in Electrical and Electronic Engineering Article information : To cite this document :. 2014. [32] P. Solin, K. Segeth, I. Dolezel, Higher-Order nite Element Methods, first ed., Chapman and Hall/CRC, 2004, p. 408. [33] C.O. Maidana, Thermo-magnetic systems for space nuclear reactors, in: SpringerBriefs in Applied Sciences and Technology, Springer International Publishing, Cham, 2014, p. 53. [34] A. Tessarolo, C. Bruzzese, Computationally efficient thermal analysis of a low-speed high-thrust linear electric actuator with a three-dimensional thermal network approach, IEEE Trans. Ind. Electron. 62 (2015) 1410–1420. [35] I.A. S.M.olyanov, P. Karban, Optimal design of MHD pump, in: 2018 ELEKTRO, 2018, pp. 1–4.