Accelerated Modeling of Transients in Electromagnetic Devices Based on Magnetoelectric Substitution Circuits

1. Introduction

Switching in electromagnetic systems causes overcurrent and overvoltage which are dangerous for equipment. Current jumps are reason of significant mechanical stresses in power transformers. Therefore, the study of transients in electromagnetic systems is relevant. Recently, scientific interest in computer calculations of transient electromagnetic processes in various components of electrical systems containing power transformers, reactors and other devices has been increasing. The use of modern computer research of electromagnetic processes can significantly reduce the financial costs and time for expensive physical modeling and solve those problems that are not available in analytical research. In this case, it becomes possible to study electromagnetic processes in individual elements of magnetic and electrical parts of objects of the electrical system, taking into account their real design features and interaction of objects in this system.

Calculation of dynamic electromagnetic fields in transformers incorporated into complex electrical systems requires a large computational resource and time. Therefore, to accelerate the modeling of dynamic electromagnetic fields, the task is often reduced to modeling electromagnetic processes in electrical and magnetic circuits, which are interconnected [1]. Modeling electrical and magnetic circuits with lumped parameters requires significantly less computer resources than modeling fields. The combined magnetic and electrical circuits constitute a so-called magnetoelectric equivalent circuit (MEEC). A single hybrid circuit is convenient for forming a common system of equations describing electrical and magnetic processes, their mutual relationship and is modeled as a single circuit. At the same time, there is no reason to carry out additional transformations that lead the equations to a purely electric or purely magnetic circuit, as was done in [2].

The increasing complexity of transformer designs, reactors and the electrical circuits in which they are incorporated, the greater detailing of design features in product models, the consideration of non-linearity, and the increased requirements for accuracy of calculations have led to the scientific task of improving the modeling of transient electromagnetic processes in complex electromagnetic devices with non-linear elements.

Purpose of the work is to develop a highly efficient method for calculating electromagnetic transients in magnetoelectric circuits for replacing a coil on a ferromagnetic core connected to a voltage source. The method is based on the use of nonlinear magnetoelectric schemes for replacing electromagnetic devices and a spectral method for representing solution functions by orthogonal polynomials. At the same time, a schematic model for the use of the spectral method was developed.

2. Basics of Transformations for Obtaining Magnetoelectric Equivalent Circuits

It is convenient to consider the basis of the transformations for obtaining magnetoelectric equivalent circuits with the help of example of calculation of processes in a circuit (Figure 1) consisting of a closed ferromagnetic core with a median line length l, on which a coil containing N_w turns is located. The main magnetization curve for the magnetic circuit without taking into account hysteresis in the ferromagnetic core is set.

If an alternating voltage source u(t) is connected to the coil, an electric current i(t) flows through the coil turns. According to Ampère’s law one can write:

$$H\cdot l=N_w\cdot i$$

(1)

where H is the magnetic field intensity.

Let’s differentiate expression (1) with respect to time:

$${\displaystyle \frac{dH}{dB}}{\displaystyle \frac{dB}{dt}}l=N_w{\displaystyle \frac{di}{dt}},$$

(2)

where B is magnetic induction.

We transform expression (2) using the designation of differential magnetic permeability ${\mu }_d=dB/dH$

$${\displaystyle \frac{l}{S{\mu }_d}}{\displaystyle \frac{d\Phi }{dt}}=N_w{\displaystyle \frac{di}{dt}},$$

(3)

where S is magnetic core cross-section area, Φ is magnetic flux.

Let’s denote the derivatives with respect to time by a stroke and Equation (3) is represented as:

$$R_d{\Phi }^{\prime }=N_w\cdot i^{\prime },$$

(4)

where differential magnetic resistance is used

$$R_d={\displaystyle \frac{l}{{\mu }_{d}S}}.$$

(5)

Equation (4) can be interpreted as follows. The so-called “magnetic current” flows through the magnetic circuit, the value of which is equal to the time derivative of the magnetic flux Φ (see Figure 2). The concept of a magnetic displacement current by analogy with an electric displacement current with a density $dD/dt$ (D is a vector of electrical induction, historically called as electric displacement) was introduced by Heaviside [1].

The right side (4) is interpreted as a magnetic voltage source controlled by the coil current derivative. The electrical circuit contains voltage source u(t), resistance R, inductance L (due to magnetic fluxes in air) and a voltage source controlled by a magnetic flux derivative Φ’ with a control coefficient N_w. Based on Kirchhoff’s voltage rule, taking into account (4), we compose the equations of state:

$$\left.\begin{array}{l}Li{\text{'}}^{\prime }+Ri+N_w\Phi \text{'}=u(t);\\ R_d\Phi \text{'}-N_{w}i{\text{'}}^{\prime }=0\end{array}\right\}$$

(6)

These equations can be solved by any numerical method for solving ordinary differential equations, for example, the implicit Euler method, Runge-Kutta, Gere methods [3,4], etc. To solve the system (6), it is necessary to attach to it the equations of calculation by the numerical method for state variables. The coefficient R_d depends of the magnetic flux Φ and it is calculated from the magnetization curve of the core. To approximate the nonlinearity of the magnetization curve, we use one of the most convenient and an accurate method for approximating curves given values at reference points. It is the method of approximation by splines [5]. Using splines to approximate nonlinearity allows to calculate derivatives faster than other numerical methods. The fact is that the derivative of the polynomial is known, and the differentiating operation of the function is reduced to changing the coefficients of the corresponding p-form.

Similarly, it is possible to create the equivalent magneto-electric circuit for various devices, including single-phase and three-phase transformers, and to work out differential equations of state according to the equivalent circuits [6].

This method was tested using the developed programs in the MATLAB system using the Gere method up to the fourth order [7]. Calculations with the help of these programs showed the following:

-

equivalent magneto-electric circuits of power transformers are very complicated;

-

transients in these transformers have very long time, and with rapidly changing components of processes;

-

the time of simulating transients is significant, which is undesirable.

-

To speed up the simulating process, it is proposed to use orthogonal polynomials.

3. Basics of Using Orthogonal Polynomials to Integrate Differential Equations

Papers [8] show the prospect of orthogonal polynomials using to integrate differential equations. The method of polynomials decomposing functions in series is called as the spectral method in these works. This is due to the analogy of the functions decomposition in Fourier series by trigonometric functions, where the set of decomposition functions is the frequency spectrum. Work [9] is devoted to spectral methods implemented on MATLAB programs. The author argues that spectral methods are usually the best tool for solving ordinary differential equations with smooth functions. Spectral methods require less computer memory than alternative methods. The spectral method is more efficient than finite difference methods.

The paper [10] proposes a spectral method of calculating transients in electric circuits based on the representation of solution functions by series of the Chebyshev, Hermite, Legendre polynomials, and also by algebraic polynomials. In this work, the spectral method was demonstrated and studied with the help of transients calculation in a simple electrical circuit. So far, the transformation of the integro-differential equation of state into the algebraic equation for current images has been performed. Equivalent circuits using images of studied currents are proposed for electric circuits. Kirchhoff’s rules for voltage and current images have been proven. As a result, a method allowing significantly increasing of the speed of transient’s simulation was obtained.

In this work it is proposed to develop and extend this method for study of transient processes in electromagnetic devices using magneto-electric equivalent circuits.

Let us briefly describe the essence of the spectral method used in [10]. The function of changing the current of time can be decomposed in series of orthogonal polynomials:

$$i(t)\approx p(t)=c_0{P}_0(t)+c_1{P}_1(t)+\mathrm{...}{c}_{N-1}{P}_{N-1}(t),$$

(7)

where P_n(t) is an orthogonal polynomial, argument t∈[a, b].

Weierstrass’s theorem states: “A function continuous on a finite closed segment can be approximated by a polynomial with a given accuracy” [11].

The Chebyshev and Legendre polynomials [11] have the range of the argument$x\in \left[-1,1\right]$. When using these polynomials, the time argument t∈[a, b] is replaced by the dimensionless argument $x\in \left[-1,1\right]$, where $x=\left(2t-\left(a+b\right)\right)/\left(b-a\right)$ . Algebraic polynomials have the range of the argument t∈[0, ∞], and the Hermite polynomials t∈[-∞, ∞]. To improve the accuracy of calculations in real time, it is advisable to normalize the values of the argument : x=t·K_n. Normalization factor is  ,

τ = b − a.

Then the expansion in series (7) takes the form:

$$i(x)\approx p(x)=c_0{P}_0(x)+c_1{P}_1(x)+c_2{P}_2(x)+\mathrm{...}{c}_{N-1}{P}_{N-1}(x)$$

Generation of Matrix Equations

We select at the time interval [a, b] a number of collocation reference points t_m (m = 0, 1, 2, …N-1), corresponding to the points x_m of the variable x on the normalized interval [-1, 1]. If condition (8) is written for each reference point, then we get a system of linear algebraic equations:

$$\left.\begin{array}{l}{c}_0{P}_0(x_0)+c_1{P}_1(x_0)+c_2{P}_2(x_0)+\mathrm{...}{c}_{N-1}{P}_{N-1}(x_0)=i_0\\ c_0{P}_0(x_1)+c_1{P}_1(x_1)+c_2{P}_2(x_1)+\mathrm{...}{c}_{N-1}{P}_{N-1}(x_1)=i(x_1)\\ \mathrm{................................}\\ c_0{P}_0(x_{N-1})+c_1{P}_1(x_{N-1})+c_2{P}_2(x_{N-1})+\mathrm{...}{c}_N{P}_{N-1}(x_{N-1})=i(x_{N-1})\end{array}\right\}$$

(9)

where the symbol P denotes an orthogonal polynomial of any type.

Subtract the first equation from the equations of system (9), then we get an reduced system:

$$\left.\begin{array}{l}{c}_1[P_1(x_1)-P_1(x_0)]+c_2[P_2(x_1)-P_2(x_0)]+\mathrm{...}{c}_{N-1}[P_{N-1}(x_1)-P_{N-1}(x_0)]=i(x_1)-i_0\\ c_1[P_1(x_2)-P_1(x_0)]+c_2[P_2(x_2)-P_2(x_0)]+\mathrm{...}{c}_{N-1}[P_{N-1}(x_2)-P_{N-1}(x_0)]=i(x_2)-i_0\\ \mathrm{....................}\\ c_1[P_1(x_{N-1})-P_1(x_0)]+c_2[P_2(x_{N-1})-P_2(x_0)]+\mathrm{...}{c}_{N-1}[P_{N-1}(x_{N-1})-P_{N-1}(x_0)]=i(x_{N-1})-i_0\end{array}\right\}.$$

(10)

We enter the column vector of polynomials as a function of x:

$$P(x)={\left[\begin{array}{cccc}{P}_1(x)& P_2(x)& \mathrm{...}& P_{N-1}(x)\end{array}\right]}^T$$

(11)

Let’s denote:

$$V=\left[\begin{array}{c}P(x_1)\\ P(x_2)\\ \mathrm{...}\\ P(x_{N-1})\end{array}\right]-\left[\begin{array}{c}P(x_0)\\ P(x_0)\\ \mathrm{...}\\ P(x_0)\end{array}\right]$$

(12)

The matrix form of system (10) takes the form as shown in [10]:

V·С = I – I₀ ,

(13)

where C = [c₁ c_2… c_N-1]^T is vector of values of coefficients of approximating polynomial (8) without coefficient c₀;

I = [i(x₁) i(x₂)… i(x_N-1)]^T is vector of current values at reference points 1, 2, ... N-1.

I₀ = [i₀ i_{0 …} i₀]^T is a vector of current values at point t₀ with dimension N-1.

Series decomposed in orthogonal polynomials can be differentiated and integrated.

Similarly, we can compose a system of linear algebraic equations for derivatives and integrals. System (9) for derivatives in matrix form, taking into account ${\displaystyle \frac{dP}{dt}}={\displaystyle \frac{dP}{dx}}\cdot {\displaystyle \frac{dx}{dt}}={\displaystyle \frac{dP}{dx}}\cdot K_n$, takes the form:

$$\left[\begin{array}{c}{i^{\prime }}_1\\ \begin{array}{l}{i^{\prime }}_2\\ \mathrm{...}\end{array}\\ {i^{\prime }}_{N-1}\end{array}\right]=\left[\begin{array}{cccc}{P^{\prime }}_1(x_1)& {P^{\prime }}_2(x_1)& \mathrm{...}& {P^{\prime }}_{N-1}(x_1)\\ \begin{array}{l}{P^{\prime }}_2(x_1)\\ \mathrm{...}\end{array}& \begin{array}{l}{P^{\prime }}_2(x_2)\\ \mathrm{...}\end{array}& \begin{array}{l}\mathrm{...}\\ \mathrm{...}\end{array}& \begin{array}{l}{P^{\prime }}_{N-1}(x_2)\\ \mathrm{...}\end{array}\\ {P^{\prime }}_{N-1}(x_1)& {P^{\prime }}_2(x_{N-1})& \mathrm{...}& {P^{\prime }}_{N-1}(x_{N-1})\end{array}\right]\cdot \left[\begin{array}{c}{c}_1\\ \begin{array}{l}{c}_2\\ \mathrm{...}\end{array}\\ c_{N-1}\end{array}\right],$$

(14)

Or

I’ = K_n·D· C,

(15)

where D is a matrix of a system of linear equations for derivatives (14),

I^′ = [i^′ (x₁) i^′ (x₂)… i^′ (x_N-1)]^T is vector of values of current derivatives at reference points k=1, 2,… N-1.

The integrals with the upper limit x_m of the function (8) in the form of a polynomial decomposition are following as:

$$J=S\cdot C+\Delta \cdot I_0.$$

(16)

Changing the numbers of reference points m from 1 up to N-1, we obtain a system of linear equations with respect to coefficients c₀ and c_k (k=1,.. m), which make up the vector C. In matrix form for all reference points, the system of linear equations for integrals has the form [10]:

$$J=S\cdot C+\Delta \cdot I_0,$$

(17)

where S, Δ are special matrices, the calculation of which for different polynomials is shown in [10]. This work shows that integro-differential equations of state for a simple electric circuit R-L-C-e are converted into algebraic equations for the image of current C:

$$(LD+RV+BS)C=U-u_{C0}-R{i}_0-B\cdot \Delta \cdot i_0,$$

(18)

where B is the inverse value of the capacitance of the condenser,

U is matrix-column of source voltage values at the range of time points..

As a result of solution (18), a vector C of coefficients of polynomial collocation of current is determined. Knowing the vector C at the known initial value of current i₀ and the initial value of voltage across the capacitor u_C0 it is possible to determine the values of current at all arbitrary points of the time segment. The interpolation error of the time function of the current, the derivative and the time integral of the current is estimated in [10].

It is shown that the interpolation error of the current time function is minimal when choosing the position of the reference points in the zeros of Chebyshev’s polynomials. This work shows that the derivative of the solution function has the greatest approximation error. When analyzing transients, we are also faced with an approximation of the current derivative. To reduce this error, it is advisable to choose the positions of the reference points not in the zeros of the Chebyshev’s polynomials, but in the maxima. Thus, the optimal choice of the position of the reference points for different tasks is not clear in advance, but can be clarified by practical calculations.

4. Using Orthogonal Polynomials to Calculate Transients in a Ferromagnetic Core Coil

4.1. General Equations

To calculate transients in a coil with a ferromagnetic core using orthogonal polynomials, we use the magnetoelectric equivalent circuit Figure 2. The equations of state are in the form (6). It is applied the basics of using orthogonal polynomials for integrating of differential equations described in section 2. To describe transients in the form of orthogonal polynomials (determination of polynomial coefficients), it is necessary to use Equations (13) and (14) for electric and magnetic currents. For the magnetic current (magnetic flux derivative), Equation (13) is as follows:

Φ΄ = V·С_Φ + Φ₀΄,

(19)

where C_Φ is the image vector of the magnetic current function;

V is matrix (12);

Φ′₀ is the magnetic current value vector (magnetic flux derivative) at the initial point t₀.

According to the method described in section 2, it is selected a number of collocation reference points t_m (m = 0, 1, 2, … N-1) at the time interval [a, b], which corresponds to points x_m of variable x at the normalized interval [-1, 1]. Using the image of the desired current function (13), magnetic current function (19), derivative of function (15), we transform the system of Equation (6) into a system of equations for images. As a result, we get a system of algebraic equations:

where C_i is a vector of decomposition coefficients of the electric current function;

D is the matrix of Equation (14);

I₀ is vector of electric current values at the initial point;

Φ′₀ is the vector of magnetic current values at the initial point.

4.2. Schematic Interpretation of the Method of Numerical Calculation of Transients in Magnetoelectric Circuits

Equation (20) can be interpreted using a special equivalent circuit in Figure 3 corresponding to the original scheme (Figure 2). Current i(t) flows in the electrical branch of the original circuit. We have the image of the current C_i instead of the current i(t) in the equivalent circuit. The image C_i is a vector of decomposition coefficients in series by the orthogonal polynomials of the current function i(t).

In the electrical branch of the equivalent circuit the resistive element has image R·V and a direct voltage source of value R·i₀ is connected into series with it. The inductive element has an image in the form of a some resistance L·K_n·D. We have the image C_Φ instead of the magnetic current Φ′(t) in the magnetic branch of equivalent circuit. The image C_Φ is a vector of decomposition coefficients of the magnetic current function in series by orthogonal polynomials. In the electrical branch there is the direct voltage source controlled by the image of the magnetic current N_w (V C_Φ +Φ′(t₀)).

The differential resistor R_d has an image R_d·V in the magnetic branch of the equivalent circuit and, in the opposite direction of the current in series to it a direct voltage source with a value R_d·Φ′( t₀) is connected. This branch has direct voltage source controlled by the image of the electric current derivative N_w ·K_n·D.

The inductance L is due to the magnetic flux passing outside of the magnetic core. The calculation of this inductance is a separate problem that could be solved by field methods using software complexes ANSYS [12], COMSOL [13], the method of contour magnetic fluxes [14] in static fields. Calculations show that this inductance is weakly dependent of current and can be taken as a constant.

Matrix form of the system (20) is:

$$Z\cdot C=F$$

(21)

where

$$Z=\left[\begin{array}{cc}-N_w\cdot D& R_d\cdot V\\ L\cdot D+R\cdot V& N_w\cdot V\end{array}\right],$$

(22)

$$F=\left[\begin{array}{c}-R_d\cdot \Phi {\text{'}}_0\\ U-R{I}_0-N_w\cdot \Phi {\text{'}}_0\end{array}\right]$$

(23)

C=[C_i· С_Φ]^T

The result of solving Equation (21) is a vector C. Using vectors C_i, C_Φ. and also the value of i₀, Φ′₀ at the initial point t₀, we can get the electric and magnetic current values in all the reference points in the time segment [a, b] according (13) and (19). Taking into account (8) we get the value of electric and magnetic currents in all arbitrary points of the segment.

The entire calculation time interval is divided by N_u segments with time step τ=b-a to reduce the calculation error. Let’s perform calculations on each segment using the cyclic run method increasing the current time by τ each time.

The calculation should take into account that the matrix Z (21) contains the value of the magnetic resistance R_d, which depends on the magnetic flux. Therefore, at the reference points on each time segment, the magnetic flux values of the branch containing the ferromagnetic core are calculated and, according to the dynamic magnetization curve, the differential magnetic resistance R_d is calculated. Since the magnetization curve is nonlinear, the calculation of the magnetic resistance on each time segment is performed in an iterative cycle. The magnetic flux values of the magnetic branch at all reference points can be calculated by the formula:

$$\Phi ={\Phi }_0+(S\cdot C_{\Phi }+\Phi {\text{'}}_0\Delta )\cdot \tau /2$$

(25)

4.3. Calculation Algorithm, Programs and Simulation Results

To approximate the nonlinear dependence of the ferromagnetic magnetization curve, we use one of the most convenient and accurate methods of curves approximation given by values at reference points. It is the method of approximation by splines. System MATLAB [5] has software support for spline functions. If you set the magnetization curve using two vectors MB, MH, which determine a coordinates of values of the reference points of magnetic induction and magnetic field strength, then using the standard MATLAB csapi-function, you can get the so-called p-form: p = csapi (MB, MH). P-form is a structure that stores all coefficients and other information about the spline approximation. Using the p-form, it is possible to calculate the value of the magnetic field strength H: H = fnval (B, p) for any value of magnetic induction B. Using the standard fnder-function and applying p-form for approximation of the function H (B), we can obtain p-form that approximates the derivative of the function dH/dB (B): dp = fnder (p). Then, for an arbitrary value of B, we can calculate the derivative value: dH/dB = fnval (B, dp).

To check the adequacy of the proposed method of transients modeling, a computer program coil_VDS in the MATLAB system has been developed. The special feature of this program is that the matrix Z, composed according to the system (21), contains the value of the magnetic resistance R_d, which depends on the magnetic flux. Therefore the magnetic flux values of the branch (magnetic current integral) containing the ferromagnetic core are calculated at the reference points on each time segment. Differential magnetic resistance is calculated according to the magnetization curve. Since the magnetization curve is nonlinear, the calculation of magnetic resistance is carried out in an iterative cycle.

Calculations are performed in the following sequence, according to the block diagram (Figure 4). Each action sequence item described below has its own block number on the block diagram

• Input of initial data (start time t_begin and end time of simulation t_end, time segment dimension τ, electric circuit parameters R, L, e(t), core coil parameters S, l, magnetization curve B (H), initial conditions Φ₀, i₀) is performed. The orthogonal polynomial type is selected.
• It is set: the number N of reference points on a segment without a zero point (7 points are recommended for each segment), the length of the time segment, the segments number N_u is calculated over the entire time investigated interval.
• To calculate the coefficients of orthogonal polynomials the positions of the reference points x_k on interval [-1, 1] of one segment (the position of the reference points on all segments is the same) are set. In the paper [10] it is shown that the reference points can be set uniformly, and also thickened to the ends of the segment or to the middle of the segment.
• The segment is specified in the interval [a, b] on the timeline, where a and b are the start and end times of the segment, with τ=b-a. Corresponding to the locations of the reference points x_k on [-1, 1] the reference points t_k on the time interval [a, b] are calculated using the formula

t_k = a +τ (x_k + 1)/2.

5.

The initial differential magnetic resistance of the R_d0 of magnetic branch is calculated from the initial value of the magnetic flux Φ₀ of the ferromagnetic.

6.

The values of matrices V, D, S, Δ for the left part of the Equation (21) are calculated.

7.

It is filled in matrix Z, in which the initial differential magnetic resistance R_d0 is used.

8.

Current calculations for each segment on time interval [a, b] are cyclically performed (items 8-18). Segment numbers are cyclically changed from 1 up to N_u. When the ku cycle parameter changes, it is performed the following acts :

9.

The values of the source voltage vector U at all points of the segment are calculated.

10.

Iterative cycle (items 10-16) of current calculation for reference points of the current segment is performed. The iterative cycle is necessary to clarify the value of the differential magnetic resistance R_d, since this value depends on the value of the magnetic flux.

11.

The matrix Z containing the value of the magnetic resistance R_d is clarified. The vector of the right parts F of the system (23) is calculated.

12.

System of algebraic equations is solved and vector of polynomial coefficients C=Z^-1·F is determined.

13.

The vectors of polynomial coefficients C_i , C_Φ for each current are distinguished.

14.

Vectors of values of all currents (including magnetic currents) are calculated based on values of vectors C_i, C_Φ according to (13).

15.

A vector of magnetic flux values at reference points is calculated using the formula (25) for calculation of integral.

16.

The differential magnetic resistance R_d of the branch is calculated from the values of the magnetic flux according to the magnetization curve of the ferromagnetic, using the function of approximation by splines of the magnetization curve. The end condition of the iteration loop is checked. The iteration cycle (items 10-16) ends if the values of magnetic resistance of adjacent iteration cycles do not exceed the specified error, otherwise we return to item 10.

17.

The last current and magnetic flux values on the current segment will become the initial current and magnetic flux values for the next segment when the iterative cycle is terminated.

18.

Current and magnetic flux values are stored in arrays for output at the end of calculation.

The condition of the end of the transient calculation is checked. If the termination condition is not met, the initial data for calculating the next segment is set, namely: the initial value of the current i₀, the initial value of the magnetic flux Φ₀ and the start time of the segment t₀ as the last values obtained in calculating of previous segment. The transient calculation ends (steps 8-18) after calculation of the last segment, when the current time of the simulating process reaches the set value t_end.

19.

The graphs of the transient in the time area are generated.

Due to the adopted simplifications of the representation of the magnetization curve (without taking into account hysteresis), the authors do not set the purpose of full compliance of modeling results with processes in a real object. The purpose of modeling is to compare the accuracy and speed of the proposed method and the well-known and widely used implicit Euler method. A computer program has also been developed to perform the same task with the help of numerical implicit Euler calculation. The integration step is chosen so that when it is reduced, the calculation results do not change. Therefore, the calculations by the Euler method can be considered as quite accurate.

The research was carried out for sinusoidal and trapezoidal source voltages. The calculation errors of the current and magnetic flux were compared and the simulation time was estimated when representing the solution functions in series by Chebyshev, Hermite, Legendre polynomials and also by algebraic polynomials. The calculations were carried out with a different choice of the method of distribution-division the reference points: in zeros of Chebyshev polynomials, in maxima of Chebyshev polynomials, and also with a uniform distribution of points. Computer program coil_VDS in the MATLAB system was developed to check the adequacy of the proposed method of transient modeling. The text of this program is given in [16].[M1]

The following object parameters are set for the simulation; u (t) = 1000√2*sin (314t) - in the first case, L = 0.0017 H, R = 1.8 Ohm. The main magnetization curve of the ferromagnetic core is given in the form of Table 1 and shown in the graph (Figure 5).

To increase the accuracy and speed of calculation of approximation of the magnetization curve, splines are used in this program. The maximum relative approximation error did not exceed 0.01%. The order of the used orthogonal polynomials was set by the number of reference points in segment: N=7.

Figure 6 shows the results of transient calculation for Chebyshev polynomials at initial currents and magnetic flux which are equal to zero. It is a fragment of the process. The transient lasts more than 10 sec.

The accuracy and calculation time for non-zero initial conditions (deeper core saturation) were studied.

The results of the calculation using the coil_VDS program under sinusoidal action for Chebyshev polynomials, and the calculation using the implicit Euler method, are shown in Figure 7 and Figure 8. Initial magnetization of ferromagnetic core is B₀=0,7 T, initial value of current i₀=0 are set.

In Figure 7 and Figure 8, the values calculated by the implicit Euler method are shown by asterisks, and the results obtained by the proposed method using the coil_VDS program are shown by solid lines. The text of this program is given in [17].

To assess the accuracy of calculations, the relative mean quadratic error evasions over the entire simulation interval are determined. For this, the quadratic deviation of the current function value calculated in the program from the value obtained by the Euler method at each reference point k is calculated:

The mean quadratic deviation is calculated for the entire simulation interval

and relative mean square deviation in percent

Errors of the magnetic induction function are similarly calculated.

The change in the percentage of relative error in the calculation process is shown in Figure 9. The reference points are distributed at the maxima of the Chebyshev polynomials. It can be seen from the figures that the values obtained by different methods coincide well and the error does not increase with increasing time value.

Table 2 shows the relative mean square deviations in percent for the current function Δi and the magnetic induction function ΔB. The calculations were carried out for a different choice of the method of distribution of reference points: in zeros of Chebyshev polynomials, in maxima of Chebyshev polynomials and also with a uniform distribution of points. Calculation time is set to 0.3 s, the number of reference points on the segment N = 7 and for all variants of calculations is the same one. Processor calculation time was estimated by tic, toc operators - Elapsed time is 13.884412 seconds.

Calculations were also carried out under the action of the trapezoidal voltage source (Figure 10). Calculation results for B₀=0.7 T, u_max=1000·1,41 V, u_min=-1000·1,41 V, R=6,8 Ohm are shown in Figure 11 and Figure 12. The text of this program is given in [18].

Table 3 shows the relative mean square deviations in percent for the current function Δi, the magnetic induction function ΔB under the action of source trapezoidal voltage.

Calculations have shown that the values obtained by the proposed method and the Euler method coincide well and the error do not increase with increasing time.

5. Discussion

The comparison of time of calculation was carried out using the program based on the proposed method and the program based on the implicit Euler method of numerical integration of differential equations. This program can be freely downloaded from the specified link [17] and tested in the MATLAB system. The processor calculation time for this program in the MATLAB system (estimated using the built-in function tic, toc) is 12.2 times longer than when using orthogonal polynomials. This can be explained as follows. Traditional methods of numerical integration of differential equations calculate sequentially current step values which are based on the results on the previous steps. When using polynomials, the polynomial function of solution is already known. Several values are calculated at the reference points of the current segment that is lead to reduction of the calculation time.

Comparison of the errors obtained in calculations by different methods showed the following. The accuracy of the methods is somewhat worse for trapezoidal action of voltage source than for sinusoidal action, but is quite acceptable for practical calculations.

The accuracy of methods using Chebyshev, Legendre, and algebraic polynomials is better than using Hermite polynomials. Therefore, the use of Hermite polynomials is not recommended for the calculation of magnetoelectric equivalent circuits. The method of distribution of reference points slightly affects the accuracy of the calculation.

Improved calculation accuracy can be achieved by reducing of the segment dimension. However, how many times the segment dimension decreases, the calculation time increases by the same number of times. Increasing the number of reference points N >7 does not improve accuracy, but increases calculation time.

In conclusion, we note that the proposed method for calculating nonlinear magneto-electric circuits of arbitrary complexity can be used to study transients in transformers included in complex electrical circuits.

6. Conclusions

1. The developed unified spectral method of calculation of electromagnetic transients in magneto-electric circuits based on representation of solution functions by orthogonal polynomials allows to significantly reduce simulation time as compared to known methods.

2. The completed development of the schematic model of the method is convenient for practice.

3. The study of the method showed that when using different polynomials (algebraic, Tchebishev, Hermite and Legendre) the accuracy of calculation changes slightly. The method of distribution of reference points slightly affects to the accuracy of the calculation.

4. Comparison of the calculation time by the proposed method and the implicit Euler method showed that the calculation time decreased by 12.2 times.