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

Abstract

During switching in electrical systems, transient electromagnetic processes occur. The resulting dangerous current surges are best studied by computer simulation. However, the time required for computer simulation of such processes is significant for complex electromagnetic devices, which is undesirable. The use of spectral methods can significantly speed up the calculation of transient processes and ensure high accuracy. At present, we are not aware of publications showing the use of spectral methods for calculating transient processes in electromagnetic devices containing ferromagnetic cores. The purpose of the work: The objective of this work is to develop a highly effective method for calculating electromagnetic transient processes in a coil with a ferromagnetic magnetic core connected to a voltage source. The method involves the use of nonlinear magnetoelectric substitution circuits for electromagnetic devices and a spectral method for representing solution functions using orthogonal polynomials. Additionally, a schematic model for applying the spectral method is developed. Obtained Results: A method for calculating transients in magnetoelectric circuits based on approximating solution functions with algebraic orthogonal polynomial series is proposed and studied. This helps to transform integro-differential state equations into linear algebraic equations for the representations of the solution functions. The developed schematic model simplifies the use of the calculation method. Representations of true electric and magnetic current functions are interpreted as direct currents in the proposed substitution circuit. Based on these methods, a computer program is created to simulate transient processes in a magnetoelectric circuit. Comparing the application of various polynomials enables the selection of the optimal polynomial type. The proposed method has advantages over other known methods. These advantages include reducing the simulation time for electromagnetic transient processes (in the examples considered, by more than 12 times than calculations using the implicit Euler method) while ensuring the same level of accuracy. The simulation of processes over a long time interval demonstrate error reduction and stabilization. This indicates the potential of the proposed method for simulating processes in more complex electromagnetic devices, (for example, transformers).

1. Introduction

Switching in electromagnetic systems causes overcurrent and overvoltage which are dangerous for equipment. Current jumps are the reason for 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.

The calculation of electromagnetic fields in transformers and reactors during transient modes is possible using very expensive software packages like ANSYS [1], COMSOL [2], and other similar software. However, not all users can afford a costly license. These packages commonly use the Finite Element Analysis method. Such software packages allow for the calculation of three-dimensional fields in transient modes when the electromagnetic devices are connected to an electrical circuit containing only a few elements. Computer simulation of dynamic electromagnetic fields in transformers included in complex electrical systems is generally not feasible.

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 [3]. Modeling electrical and magnetic circuits with lumped parameters requires significantly fewer computer resources than modeling fields. The combined magnetic and electrical circuits constitute a so-called magnetoelectric equivalent circuit. A single hybrid circuit is convenient for forming a common system of equations describing electrical and magnetic processes as well as their mutual relationship, and it 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 in [4].

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.

The 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 is 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 the 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 magnetic field intensity.

Let us 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\mathrm{\Phi }}{dt}}=N_w{\displaystyle \frac{di}{dt}},$$

(3)

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

Let us denote the derivatives with respect to time by a stroke. Equation (3) is represented as

$$R_d{\mathrm{\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 magnetic flux Φ (see Figure 2). The concept of a magnetic displacement current by analogy with an electric displacement current with density $dD/dt$ (D is a vector of electrical induction, historically called electric displacement) was introduced by Heaviside [3].

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 control coefficient N_w. Based on Kirchhoff’s voltage rule, taking into account (4), we compose the following equations of state:

$$\left.\begin{array}{l}L{i}^{\prime }+Ri+N_w{\mathrm{\Phi }}^{\prime }=u(t);\\ R_d{\mathrm{\Phi }}^{\prime }-N_w{i}^{\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 [5,6], etc. To solve System (6), it is necessary to attach to it the equations of calculation by the numerical method for state variables. Coefficient R_d depends on magnetic flux Φ and 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 accurate methods for approximating curves given values at reference points. It is the method of approximation by splines [7]. 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 equivalent circuits [8].

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

-

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

-

Transients in these transformers require a very long time and 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

In book [6], the foundations of the spectral method for solving differential equations are provided, highlighting the advantages of this method over the finite difference method. The author utilizes Chebyshev polynomials. In work [10], the prospects for using orthogonal polynomials to solve differential equations are demonstrated. This book aims to teach spectral methods for boundary value problems, eigenvalue problems, and time-dependent problems. It discusses Hermite, Laguerre, the rational of Chebyshev, sinc, and spherical harmonic functions. The method of decomposing functions into series of orthogonal polynomials is called by the author the spectral method, by analogy with the decomposition of functions into Fourier series using trigonometric functions. In this case, the collection of decomposition coefficients forms a spectrum. Boyd’s work also examines the convergence rate of series and includes a rich bibliography on the application of orthogonal polynomials for solving differential equations.

In [11], spectral methods for solving differential equations are detailed. The author notes the exceptional accuracy of spectral methods and their intensive study over the last few decades. The author uses the MATLAB program for numerical experiments.

Ref. [12] provides a brief introduction to spectral methods, with special attention to practical examples. For simplicity, the author limits the discussion to one-dimensional solvers. It is shown that spectral methods are more advanced than finite difference schemes. They can solve a wide range of problems. These methods can achieve high accuracy with moderate computational resources.

In [13], the foundations of spectral methods for solving partial differential equations are presented. The author demonstrates that spectral methods are not overly complex.

Ref. [14] focuses on spectral methods implemented in MATLAB programs, asserting that spectral methods are usually the best tools for solving ordinary differential equations with smooth functions. They require less computer memory and are more efficient than finite difference methods.

In work [15], a method for calculating transients in a simple electrical circuit based on representing solution functions as series of Chebyshev, Hermite, Legendre, and algebraic polynomials is presented. This work demonstrates the transformation of itegro-differential state equations into algebraic equations for current representations. Substitution circuits using the representations of the investigated currents are proposed. Kirchhoff’s laws for voltage and current representations are proven. This significantly increases the speed of computer simulation for transient processes. In this study, the method is proposed to be further developed and extended for studying transient processes in electromagnetic devices using magnetoelectric substitution circuits.

Let us briefly describe the essence of the spectral method used in [15]. 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)+\dots c_{N-1}{P}_{N-1}(t),$$

(7)

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

According to the Weierstrass theorem, “A function continuous on a finite closed interval can be approximated with any desired precision by a polynomial” (Bateman, [16]).

The domain of definition for the arguments of Chebyshev and Legendre polynomials is x ∈ [−1, 1]. When using these polynomials, the time argument t ∈ [a, b] must be replaced with a dimensionless argument x. Algebraic polynomials are defined on t ∈ [0, ∞], while Hermite polynomials are defined on t ∈ [−∞, ∞]. For calculations involving real-time argument values, normalization is required: x = t·K_n. The normalization coefficient is $K_n=2/\tau$, τ = b − a. Thus, the decomposition of the function describing current variation over time into a series (7) takes the following form:

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

(8)

Generation of Matrix Equations

Over the interval of the time argument [a, b], we define a set of collocation reference points t_m (m = 0, 1, 2, … N − 1). These correspond to points x_m of variable x on the normalized interval [−1, 1]. Equation (8) can then be written for each collocation point. As a result, a system of linear algebraic equations is obtained:

$$\left.\begin{array}{l}{c}_0{P}_0(x_0)+c_1{P}_1(x_0)+c_2{P}_2(x_0)+\dots 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)+\dots c_{N-1}{P}_{N-1}(x_1)=i(x_1)\\ \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \\ c_0{P}_0(x_{N-1})+c_1{P}_1(x_{N-1})+c_2{P}_2(x_{N-1})+\dots c_N{P}_{N-1}(x_{N-1})=i(x_{N-1})\end{array}\right\},$$

(9)

where symbol P denotes an orthogonal polynomial of any type.

• From all the equations in System (9), we subtract the first equation. As a result, we obtain a 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)]+\dots 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)]+\dots c_{N-1}[P_{N-1}(x_2)-P_{N-1}(x_0)]=i(x_2)-i_0\\ \dots \dots \dots \dots \dots \dots \\ c_1[P_1(x_{N-1})-P_1(x_0)]+c_2[P_2(x_{N-1})-P_2(x_0)]+\dots 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[P_1(x)\hspace{1em}{P}_2(x)\hspace{1em}\dots \hspace{1em}{P}_{N-1}(x)\right]}^T.$$

(11)

Let us denote

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

(12)

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

V·C = I − I₀,

(13)

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

I = [i(x₁) i(x₂)… i(x_N−1)]^T is a 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 the 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\\ \dots \end{array}\\ {i^{\prime }}_{N-1}\end{array}\right]=\left[\begin{array}{cccc}{P^{\prime }}_1(x_1)& {P^{\prime }}_2(x_1)& \dots & {P^{\prime }}_{N-1}(x_1)\\ \begin{array}{l}\hspace{0.17em}{P^{\prime }}_2(x_1)\\ \dots \end{array}& \begin{array}{l}{P^{\prime }}_2(x_2)\\ \dots \end{array}& \begin{array}{l}\dots \\ \dots \end{array}& \begin{array}{l}{P^{\prime }}_{N-1}(x_2)\\ \dots \end{array}\\ {P^{\prime }}_{N-1}(x_1)& {P^{\prime }}_2(x_{N-1})& \dots & {P^{\prime }}_{N-1}(x_{N-1})\end{array}\right]\cdot \left[\begin{array}{c}{c}_1\\ \begin{array}{l}{c}_2\\ \dots \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 a vector of values of current derivatives at reference points k = 1, 2,… N – 1.

In [16], it was shown that the integrals of Legendre polynomials are determined by the following formulas:

$$\begin{gathered} {\displaystyle \int P_n(x)dx}={\displaystyle \frac{P_{n+1}(x)-P_{n-1}(x)}{2n+1}};{\displaystyle \int P_1(x)dx}=P_2(x)/3+d; \\ {\displaystyle \int P_0(x)dx}=P_1(x)+d. \end{gathered}$$

(16)

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

$$J\left(x_m\right)=\underset{x_0}{\overset{x_m}{\mathrm{\int }}}\left(c_0{P}_0(x)+c_1{P}_1(x)+c_2{P}_2(x)+\dots c_{N-1}{P}_{N-1}(x)\right)dx=c_0\underset{x_0}{\overset{x_m}{\int }}{P}_0(x)dx+\underset{x_0}{\overset{x_m}{\int }}(\mathbf{P}·\mathbf{C})dx$$

(17)

From Expression (17), considering Formula (16), we obtain

$$\begin{array}{l}J(x_m)={\displaystyle \underset{x_0}{\overset{x_m}{\int }}p(x)dx}=\left[P_1(x_m)-P_1(x_0)\right]c_0+\left\{{\displaystyle \frac{1}{3}}\left[P_2(x_m)-P_2(x_0)\right]\right\}{c}_1+\mathrm{L}\\ +\left\{\left[{\displaystyle \frac{P_{k+1}(x_m)-P_{k+1}(x_0)}{2k+1}}-{\displaystyle \frac{P_{k-1}(x_m)-P_{k-1}(x_0)}{2k+1}}\right]\right\}{c}_k+\mathrm{L}\\ +\left\{\left[{\displaystyle \frac{P_N(x_m)-P_N(x_0)}{2N-1}}-{\displaystyle \frac{P_{N-2}(x_m)-P_{N-2}(x_0)}{2N-1}}\right]\right\}{c}_{N-1}.\end{array}.$$

(18)

Coefficient c₀ is determined from the first equation of System (8) as follows:

$$c_0=i_0-\left(c_1{P}_1(x_0)+c_2{P}_2(x_0)+\dots \hspace{0.17em}{c}_{N-1}{P}_{N-1}(x_0)\right)=i_0-P(x_0)\cdot C.$$

(19)

In (18), we substitute Expression (19) for c₀. We obtain

$$\begin{array}{l}J(x_m)={\displaystyle \underset{x_0}{\overset{x_m}{\int }}p(x)dx}=\left[P_1(x_m)-P_1(x_0)\right]c_0+\left\{{\displaystyle \frac{1}{3}}\left[P_2(x_m)-P_2(x_0)\right]\right\}{c}_1+\dots \\ +\left\{\left[{\displaystyle \frac{P_{k+1}(x_m)-P_{k+1}(x_0)}{2k+1}}-{\displaystyle \frac{P_{k-1}(x_m)-P_{k-1}(x_0)}{2k+1}}\right]\right\}{c}_k+\dots \\ +\left\{\left[{\displaystyle \frac{P_N(x_m)-P_N(x_0)}{2N-1}}-{\displaystyle \frac{P_{N-2}(x_m)-P_{N-2}(x_0)}{2N-1}}\right]\right\}{c}_{N-1}\hspace{0.17em}.\end{array}$$

(20)

Let us express Formula (20) through vector V:

$$\begin{array}{l}J(x_m)={\delta }_m{i}_0+\left\{{\displaystyle \frac{1}{3}}\left[V_2(x_m)\right]-P_1(x_0){\delta }_m\right\}{c}_1+\dots \\ +\left\{\left[{\displaystyle \frac{V_{k+1}(x_m)}{2k+1}}-{\displaystyle \frac{V_{k-1}(x_m)}{2k+1}}\right]-P_k(x_0){\delta }_m\right\}{c}_k+\dots \\ +\left\{\left[{\displaystyle \frac{V_N(x_m)}{2N-1}}-{\displaystyle \frac{V_{N-2}(x_m)}{2N-1}}\right]-P_{N-1}(x_0){\delta }_m\right\}{c}_{N-1}\hspace{0.17em}.\end{array}$$

(21)

The row vector of all terms (except the first) in Expression (21), as a function of argument x, has the form

$$S_x(x_m)=\left[\left\{{\displaystyle \frac{1}{3}}{V}_2(x_m)\right\}\dots \left\{{\displaystyle \frac{V_{k+1}(x_m)}{2k+1}}-{\displaystyle \frac{V_{k-1}(x_m)}{2k+1}}\right\}\dots \left\{{\displaystyle \frac{V_N(x_m)}{2N-1}}-{\displaystyle \frac{V_{N-2}(x_m)}{2N-1}}\right\}\right].$$

(22)

Expression (21) takes the form

$$J(x_m)={\displaystyle \underset{x_0}{\overset{x_m}{\int }}p(x)dx}={\delta }_m{i}_0+\left(S_{x\hspace{0.17em}m}-{\delta }_m\cdot P(x_0)\right)\cdot C$$

(23)

where

$${\delta }_m=P_1(x_m)-P_1(x_0)=x_m-x_0$$

(24)

For collocation points m = 1,2, …N − 1, Integral (23) takes the form

$$J(x_1)={\delta }_1{i}_0+\left(S_{x\hspace{0.17em}10}-{\delta }_1\cdot P_1(x_0\right)\cdot C$$

$$J(x_n)={\delta }_n{i}_0+\left(S_{x\hspace{0.17em}n0}-{\delta }_n\cdot P_n(x_0\right)\cdot C$$

$$J(x_{N-1})={\delta }_{N-1}{i}_0+\left(S_{x\hspace{0.17em}N-10}-{\delta }_{N-1}\cdot P_{N-1}(x_0\right)\cdot C.$$

These equations in the matrix form take the following form:

$$\begin{gathered} J=\left[\begin{array}{c}\begin{array}{l}{J}_1\\ \dots \end{array}\\ \begin{array}{l}{J}_n\\ \dots \end{array}\\ J_{N-1}\end{array}\right]=\left[\begin{array}{c}{S}_{x\hspace{0.17em}1,0}-{\delta }_1{P}_1(x_0) \\ \dots \\ S_{x\hspace{0.17em}n,0}-{\delta }_2{P}_n(x_0)\\ \dots \\ S_{x\hspace{0.17em}N-1,0})-{\delta }_{N-}{P}_{N-1}(x_0)\end{array}\right]\cdot \left[\begin{array}{c}\begin{array}{l}{c}_1\\ \dots \end{array}\\ \begin{array}{l}{c}_n\\ \dots \end{array}\\ c_{N-1}\end{array}\right]+\Delta \cdot I_0, \end{gathered}$$

or

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

(25)

where

$$\begin{gathered} S=\left[\begin{array}{c}{S}_{x\hspace{0.17em}1,0}-{\delta }_{1}P(x_0) \\ \dots \\ S_{x\hspace{0.17em}n,0}-{\delta }_{2}P(x_0)\\ \dots \\ S_{x\hspace{0.17em}N-1,0}-{\delta }_{N-}P(x_0)\end{array}\right], \end{gathered}$$

(26)

$$\Delta ={\left[{\delta }_1\hspace{1em}{\delta }_2\hspace{1em}\dots \hspace{1em}{\delta }_{N-1}\right]}^T$$

(27)

For algebraic polynomials, as well as for Chebyshev and Hermite polynomials, the calculations are similar.

In work [16], it is shown that the integro-differential state equations for a simple R-L-C electrical circuit are transformed 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,$$

(28)

where symbol B denotes the inverse capacitance of the capacitor;

U is the voltage matrix at the collocation points.

The solution of algebraic System (28) determines vector C, which represents the coefficients of the polynomial approximation for the current. The values of vector C, given the initial current i₀ and initial voltage across capacitor u_C0, allow for the calculation of the current values at any arbitrary points on the time segment. In [16], the interpolation error of the time function for current, its derivative, and the integral over time was estimated. In [16], it was shown that the interpolation error of the current function i(t) is minimal when the reference points are selected at the zeros of Chebyshev polynomials. This work also showed that the greatest approximation error occurs not in the function itself but in its derivative. To reduce the approximation error of the derivative, it is advisable to select the reference points not at the zeros of Chebyshev polynomials but at the points of their maxima. The optimal placement of reference points for different tasks is not always clear ahead of time and may be determined 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 of (6). The basics of using orthogonal polynomials for integrating of differential equations described in Section 2 are applied. To describe transients in the form of orthogonal polynomials (determination of polynomial coefficients), it is necessary to use Equations (13) and (15) for electric and magnetic currents. For the magnetic current (magnetic flux derivative), Equation (13) is as follows:

Φ′ = V·C_Φ + Φ′₀

(29)

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 initial point t₀.

According to the method described in Section 2, a number of collocation reference points t_m (m = 0, 1, 2, … N − 1) is selected at 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 (29), derivative of Function (15), we transform the system of Equations (6) into a system of equations for images. As a result, we obtain a system of algebraic equations:

$$\left.\begin{array}{l}L\cdot D\cdot C_i+R\cdot V\cdot C_i+N_w\cdot V\cdot C_{\Phi }=e-R{I}_0-N_w\cdot {{\Phi }^{\prime }}_0;\\ R_d\cdot V\cdot C_{\Phi }-N_w\cdot D\cdot C_i=-R_d\cdot \Phi {\prime }_0\end{array}\right\}$$

(30)

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

D is the matrix of Equation (14);

I₀ is the 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

Equations (30) 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 current C_i instead of current i(t) in the equivalent circuit. Image C_i is a vector of decomposition coefficients in series by the orthogonal polynomials of current function i(t).

Having obtained this vector, we can reproduce function i(t) with sufficient accuracy according to Expression (7).

In the original circuit (Figure 2), the voltage across the resistive element, according to Ohm’s law, is R·i(t). In the substitution circuit, the analog of the voltage across the resistive element is R V, with an additional constant voltage source of R i₀. Coefficient R V in Equation (30), along with additional source R i₀, is referred to as the representation of the resistance of the resistive element. A similar approach is applied to the inductive element.

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 some resistance L·K_n·D. We have image C_Φ instead of magnetic current Φ′(t) in the magnetic branch of the equivalent circuit. Image C_Φ is a vector of decomposition coefficients of the magnetic current function in series by orthogonal polynomials. In the electrical branch, there is direct voltage source controlled by the image of magnetic current N_w (V C_Φ + Φ′(t₀)).

Differential resistor R_d has 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 value R_d·Φ′(t₀) is connected. This branch has direct voltage source controlled by the image of electric current derivative N_w·K_n·D.

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 [1], COMSOL [2], and the method of contour magnetic fluxes [17] in static fields. Calculations show that this inductance is weakly dependent on the current and can be taken as a constant.

The matrix form of System (30) is

$$Z\cdot C=F,$$

(31)

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],$$

(32)

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

(33)

C = [C_i, C_Φ]^T

(34)

The result of solving Equation (31) is vector C. Using vectors C_i, C_Φ, and also the value of i₀, Φ′₀ at the initial point t₀, we can obtain the electric and magnetic current values in all the reference points in the time segment [a, b] according to (13) and (29). Taking into account (8), we obtain 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 calculation error. Let us perform calculations on each segment using the cyclic run method increasing the current time by τ each time.

The calculation should take into account that Matrix Z (32) contains the value of 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 following formula:

$$\mathit{\Phi }={\mathbf{\Phi }}_0+(\mathit{S}\cdot {\mathit{C}}_{\Phi }+{{\mathbf{\Phi }}^{\prime }}_0\mathit{\Delta })\cdot \tau /2$$

(35)

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 curve approximation given by values at reference points. It is the method of approximation by splines. The MATLAB system [18] has software support for spline functions. If the magnetization curve is set using two vectors MB and MH, which determine coordinates of values of the reference points of magnetic induction and magnetic field strength, then, using the standard MATLAB csapi-function, the so-called p-form can be obtained: 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 magnetic field strength H:H = fnval (B, p) for any value of magnetic induction B. Using the standard fnder-function and applying the p-form for approximation of function H (B), we can obtain a p-form that approximates the derivative of 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 transient modeling, computer program coil_VDS in the MATLAB system is developed. The special feature of this program is that Matrix Z, composed according to System (31), contains the value of magnetic resistance R_d, which depends on 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.
• is the following are set: the number N of reference points on a segment without a zero point (seven points are recommended for each segment), the length of the time segment, the segment number N_u calculated over the entire time investigated interval.
• To calculate the coefficients of orthogonal polynomials, the positions of 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 paper [16], 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 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 reference points x_k on [−1, 1], reference points t_k on time interval [a, b] are calculated using the following formula:

  t_k = a + τ (x_k + 1)/2.
• The initial differential magnetic resistance of the R_d0 of the magnetic branch is calculated from the initial value of magnetic flux Φ₀ of the ferromagnetic.
• The values of matrices V, D, S, ∆ for the left part of Equation (31) are calculated.
• Matrix Z is filled, in which the initial differential magnetic resistance R_d0 is used.
• Current calculations for each segment on time interval [a, b] are cyclically performed (Items 8–18). Segment numbers are cyclically changed from one up to N_u. When the ku cycle parameter changes, the following acts are performed.
• The values of the source voltage vector U at all points of the segment are calculated.
• The 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 differential magnetic resistance R_d, since this value depends on the value of the magnetic flux.
• Matrix Z containing the value of magnetic resistance R_d is clarified. The vector of the right parts F of System (33) is calculated.
• The system of algebraic equations is solved and the vector of polynomial coefficients C = Z⁻¹·F is determined.
• The vectors of polynomial coefficients C_i, C_Φ for each current are distinguished.
• Vectors of values of all currents (including magnetic currents) are calculated based on values of vectors C_i, C_Φ according to (13).
• A vector of magnetic flux values at reference points is calculated using Formula (35) for calculation of the integral.
• 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.
• The last current and magnetic flux values on the current segment become the initial current and magnetic flux values for the next segment when the iterative cycle is terminated.
• 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 current i₀, the initial value of magnetic flux Φ₀, and the start time of segment t₀ as the last values obtained in calculating the 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.
• 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 is also developed to perform the same task with the help of a 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 quite accurate.

The research was conducted for sinusoidal and trapezoidal voltage sources. The current and magnetic flux calculation errors were compared, and the simulation time was assessed for solution functions represented by Chebyshev, Hermite, Legendre, and algebraic polynomials. Different methods for distributing collocation points were applied.

Collocation points were applied at zeros and maxima of Chebyshev polynomials, as well as uniformly distributed points. To verify the adequacy of the proposed method for transient process modeling, a computer program, coil_VDS, was developed in the MATLAB system. The text of this program is provided in [19].

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. Characteristic of the main magnetization curve for steel 3408-03.

H (A/m)	0	1	7.58	10.8	15.2	20.8	23.2	26.2	31.9	51.4	97.3	520.7	1218	1.25 × 105
B (T)	0	0.0096	0.1	0.2	0.4	0.8	1.0	1.2	1.4	1.6	1.7	1.8	1.84	2

and shown in the graph (Figure 5).

To increase the accuracy and speed of calculation of approximation of the magnetization curve, splines were 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 a 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 s.

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. The initial magnetization of ferromagnetic core, B₀ = 0.7 T, and the initial value of current, i₀ = 0, were set.

In Figure 7 and Figure 8 (similarly, in Figures 11 and 12), 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 [20].

To assess the accuracy of calculations, the relative mean quadratic error evasions over the entire simulation interval were 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 was calculated:

$$\Delta i_k={(i_k-i_{eul\hspace{0.17em}k})}^2.$$

The mean quadratic deviation was calculated for the entire simulation interval,

$$\Delta i=\sqrt{{\displaystyle \frac{1}{\left(N-1\right)\cdot Nu}}{\displaystyle \sum _{k=1}^{\left(N-1\right)NU}\Delta i_k}},$$

as well as the relative mean square deviation in percent,

$$\Delta i_{rel}\%={\displaystyle \frac{\Delta i}{{\left|I\right|}_{\max }}}$$

Errors of the magnetic induction function were 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.

The relative mean square deviations for current function Δi and magnetic induction function ΔB were calculated. These values are presented in

Table 2. Relative mean quadratic deviations of simulation errors in percent under sinusoidal action of voltage source.

ID	∆ %	Algebraic Polynomials			Chebyshev Polynomials			Legendre Polynomials			Hermite Polynomials
B0, T		Zeros	Max	Uniform	Zeros	Max	Uniform	Zeros	Max	Uniform	Zeros	MAX	Uniform
0	Δi	0.0036	0.0036	0.0035	0.0036	0.0022	0.0033	0.0036	0.0309	0.0789	0.0480	0.0414	0.0397
	ΔB	0.0034	0.00346	0.00347	0.0034	0.0035	0.0035	0.0034	0.0039	0.0049	0.0035	0.0034	0.0035
0.7	Δi	0.1018	0.1032	0.1087	0.1018	0.1031	0.1085	0.1018	0.1209	0.1964	0.111	0.1117	0.1159
	ΔB	0.0092	0.0093	0.0092	0.0092	0.0093	0.0092	0.0092	0.0085	0.0203	0.0099	0.0093	0.0092

. Various methods for distributing collocation points were used: at the zeros and maxima of Chebyshev polynomials, as well as uniformly distributed points. The simulation time for the interval was set at 0.3 s. The number of collocation points per segment was consistently chosen as N = 7 for all calculations. The processor time for the calculations was evaluated using operators tic, toc; Elapsed time was 13.884 s.

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 [21].

Figure 10. Source voltage as function of time.

Figure 10. Source voltage as function of time.

Figure 11. Current in coil as function of time under trapezoidal action.

Figure 11. Current in coil as function of time under trapezoidal action.

Figure 12. Magnetic induction in coil as function of time under trapezoidal action.

Figure 12. Magnetic induction in coil as function of time under trapezoidal action.

Table 3. Relative mean square deviations of simulation errors in percent under the action of source trapezoidal voltage.

∆ %	Algebraic Polynomials			Chebyshev Polynomials			Legendre Polynomials			Hermite Polynomials
	Zeros	Max	Uniform	Zeros	Max	Uniform	Zeros	Max	Uniform	Zeros	Max	Uniform
Δi	1.0325	1.0458	1.1012	1.0324	1.0457	1.1012	1.0323	1.0451	1.1007	1.0553	1.0539	1.583
ΔB	0.280	0.2841	0.294	0.2807	0.284	0.2944	0.2792	0.2436	0.29632	1.0371	0.7169	2.962

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

Calculations showed that the values obtained by the proposed method and the Euler method coincide well and the errors do not increase with increasing time.

5. Discussion

To assess the efficiency of the proposed method, a comparison of processor time was performed using the developed programs. Calculations were conducted using a program based on the proposed method and a program based on the implicit Euler numerical integration method for differential equations. The latter program can be freely downloaded via the link provided in [20] and tested in the MATLAB system. The processor time for this program in the MATLAB system (measured using the built-in functions tic, toc) is 12.2 times longer than that required for calculations using orthogonal polynomials. This can be explained as follows: well-known numerical integration methods for differential equations calculate the solution step by step, based on values from the current or previous steps. When using polynomials, the solution function is predefined as a polynomial. On a single current segment, several values at reference points are computed simultaneously, significantly reducing 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 it 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 calculation.

Improved calculation accuracy can be achieved by reducing segment dimension. However, the magnitude of segment dimension decrease is directly proportional to calculation time increase. 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

• In electromagnetic devices, the calculation of transient processes is effectively performed using magnetoelectric substitution circuits. The developed spectral method for calculating electromagnetic transients in magnetoelectric circuits, based on representing solution functions with orthogonal polynomials, significantly reduces processor time compared to known methods.
• The proposed schematic model of the method is convenient for composing state equations.
• Computational studies of the method showed that when using different polynomials (algebraic, Chebyshev, Hermite, and Legendre), the accuracy of the calculation changes insignificantly. The method of distributing collocation points has a slight effect on calculation accuracy. The processor time for calculating transients using the proposed method, compared to the implicit Euler method, showed a reduction in processor time by a factor of 12.2. This fact is reflected in the title of the article: “Accelerated Modeling”, etc.
• Comparison of the calculation time by the proposed method and the implicit Euler method showed that calculation time decreased by 12.2 times.