Mathematical Model of a Nonlinear Electromagnetic Circuit Based on the Modified Hamilton–Ostrogradsky Principle

Abstract

This paper presents a mathematical model of a typical lumped-parameter electromagnetic assembly, which consists of two subassemblies: one includes a magnetic circuit and the other with selected elements of electric circuits. An interdisciplinary research approach is used, which assumes the use of a modified integral method based on the variational Hamilton–Ostrogradsky principle. The modification of the method is the extension of the Lagrange function by two components. The first one reflects the dissipation of electromagnetic energy in the system, while the second one reflects the effect of external non-potential forces acting on the electromagnetic system. This approach allows for the avoidance of the inconvenience of the classical theory, which assumes the decomposition of the entire integrated system into individual electrical subsystems. The state equations of the electromagnetic subassembly are presented solely on the basis of the energy approach, which in turn allows taking into account various latent motions in the system, because the equations are derived based on non-stationary constraints between subsystems. The adopted theory allows for the formulation of the model of the system in a vector form, which gives much more possibilities for the analysis of higher-order electromagnetic circuits. Another important advantage is that the state equations of the considered electrical object are given in Cauchy normal form. In this way, the equations can be integrated both explicitly and implicitly. The results of computer simulations are presented in graphical form, analysed, and discussed.

1. Introduction

Equations of nonlinear electromagnetic circuits play a very important role in applied electrical engineering. Most electrical devices contain magnetic circuits, usually made in the form of cores with windings. The formulation of a mathematical model encounters the problem of taking into account physical processes in the magnetic cores of electrical machines and devices supplied with alternating voltage. It should be noted that in the magnetisation and demagnetisation processes of ferromagnetic elements, deformed induced voltages are generated in the windings due to the effect of magnetic hysteresis and saturation, i.e., phenomena represented by a nonlinear magnetisation curve. All nonlinear processes should be represented in the model of electrical devices.

Mathematical modelling of transient processes in complex dynamic systems is important for both the economy and science. Production of machines or devices is impossible without prior analysis of various processes in power systems. Transient processes are analysed at the stages of both design and operation.

There are many approaches to modelling technical objects. Two are the most common. In the first, the mathematical model is based on classical approaches [1,2,3] and the law of conservation of energy; the second is based on variational approaches [1,2,4,5] and minimisation of the action functional, according to the Hamilton–Ostrogradsky integral variational principle. Each of them has its advantages and disadvantages. In the classical approach, a single integrated system is decomposed into subsystems assigned to individual scientific fields, interconnected and described with stationary equations. The variational approach is clearly different. The so-called action functional is created for the single integrated system by forming elements of the Lagrange function; and then, this functional is minimised. The result is the Euler–Lagrange equations, which create a mathematical model of the dynamic object.

The classical application of variational approaches, using a conservative Lagrangian, in which there is a difference between kinetic and potential energy, does not allow for the formulation of a mathematical model. This Lagrangian ignores the work of dissipative forces and external non-potential forces. Such approaches usually describe conservative systems. A mathematical model of a dynamical system based on these assumptions will be practically useless in the analysis of dynamical processes; variational approaches applied in a classical way are not able to match classical methods. Variational approaches should be somewhat modified to be used effectively. In this case, the idea of modifying variational approaches was used. For this purpose, a modified variational method [2] was used, based on the well-known Hamilton–Ostrogradsky integral principle and fully reflecting the Maupertuis’s least action principle. The concept of modifying the Hamilton–Ostrogradsky principle will be briefly presented in the third chapter.

Kirchhoff’s laws, applicable to both electric and magnetic circuits, are among the fundamental principles of applied electric engineering. MMFs in the magnetic circuits of machinery and apparatus are generated by currents in the electric wiring, which makes us conclude there is a complex mutual dependence between electric and magnetic circuits. A single electromagnetic circuit arises. This paper presents a method of obtaining a mathematical model for nonlinear electromagnetic circuits not via classical approaches, but in a different way using variational approaches [2,6] based on the modified Hamilton–Ostrogradsky principle.

This paper aims to derive equations of nonlinear electromagnetic circuits based on variational approaches and introduce the concept of mathematical modelling of transient processes in nonlinear electromagnetic circuits. For this purpose, the modified Hamilton–Ostrogradsky principle is used.

There are many ways to approximate and model the magnetic hysteresis loop, including the magnetisation curve [7]. The model introduced by Jiles and Atherton in 1986 [8] is widely used as a tool for modelling magnetic and electromagnetic circuits. This model can be used to describe the hysteresis loop of any ferromagnetic material. Many dynamic models based on the original Jiles–Atherton approach have been described, but they show some uncertainty in their energetic considerations [9]. Two reverse dynamic models based on the Jiles–Atherton approach and field separation theory are presented in [9]. The energetic aspects of the model are explained, and an alternative approach based on field separation is proposed to account for dynamic losses. Ref. [10] presents a modeling of dynamic hysteresis loops. It is based on the method proposed in [8]. The dynamic effects of eddy currents [10] are considered by means of the approach discussed in [11]. A satisfactory agreement between the measured and modelled dynamic hysteresis loops and their corresponding loss densities is achieved.

Yang, W. et al. [12] suggests a mathematical model of a transformer which the authors use to analyse transient electromagnetic processes, in particular investigating the effects of the transformer’s magnetic field on the power grid. The field part of the model relies on the finite element method employed in numerical calculations. The results of the analysis of voltage and current variations across short-circuited transformer winding and of the dynamic power control in external networks are presented.

Jang, S.M. et al. [13] compare and analyse the characteristics of high-speed permanent magnet synchronous machines (PMSMs) with diametrically magnetized rotors in line with the conditions of electromagnetic circuit design. Analyses of magnetic field distribution and electric parameter calculations are utilized to this end. Three types of diametrical magnet PMSMs are additionally suggested. The models allow not only for analysing the magnetic field distribution but also assessing and comparing electric parameters. The reliability of simulation results is corroborated with experimental data that show the model is fit for use.

Ramirez-Núñez, J. A. et al. [14] describes the theoretical aspects and experimental analysis of magnetic processes in asynchronous motors with respect to fault finding and later operation. The simplicity, low costs, and potential of the method make it promising for additional diagnostics provided by other well-known methods. In particular, an advanced algorithm is presented that improves the visualisation of higher harmonics associated with a range of motor defects. An optimised and improved design of a permanent magnet electric machine stator is presented in [15]. A mathematical model is developed for testing magnetic processes and magnetic field distribution, which allows for realising the magnetic field distribution in an optimised electric machine. The mathematical model relies on the electromagnetic field theory and uses the finite element method.

A mathematical model of a four-pole magnetic bearing is introduced in [16]. The model is based on the magnetic circuits theory, with the finite element method serving to visualise the field processes. The model is implemented in a MatLab/Simulink ver. R2022b software complex in order to present the computer simulation results. They are shown as the waveforms of the magnetic bearing’s electric and magnetic parameters.

In [17], transient magnetic fields are analysed, generated via a lightning current in the presence of 3D conducting plates. An equivalent circuit is created for the plates for the period of simulation. Several methods are proposed as a result in order to mitigate the calculation complexity as large plates are simulated. The numerical and experimental test results for the method are described as well.

A transformer’s mathematical model is built in [18] for the purpose of analysing its electric and magnetic processes. For numerical modelling, the model is implemented in MatLab/Simulink, and a series of computer simulations are undertaken. Processes in the transformer working with the power grid are examined in particular. Normal operation modes when loaded, idle, and short-circuited are tested. The waveforms of the transformer’s electric and magnetic quantities are presented. Some theoretical approaches to process modelling in magnetic circuits of any configuration are developed in [19]. The mathematical magnetic core models were based on magnetic and electric equivalent circuits of appropriate configurations; then, the equations of magnetic and electric equilibrium are formulated. Some instances of magnetic core modelling of diverse configurations and computer simulation results are presented.

Transient magnetic field and voltage induced in a photovoltaic system via lightning strikes are calculated in [20]. An effective method of computing the transient magnetic field and the voltage induced is suggested. This paper presents a description of the equivalent circuit model of a support system for photovoltaic panels. The model for analysing the magnetic field is derived from the vectoral potential for sloping, vertical, and horizontal branches of the photovoltaic panel frame system. Magnetic field distribution and voltages induced via lightning strike are calculated.

Roppert, K. et al. [21] suggests modelling nonlinear stationary magnetic processes in induction heating. An original approach is advanced for modelling of a fully coupled, nonlinear magneto-thermodynamic system. The modelling uses a harmonic balancing circuit, with a special emphasis on the selection of ferromagnetic materials. The model is implemented by means of the finite difference method and a computer simulation is carried out. An investigation of magnetic processes in a cylinder and a flat plate to improve and optimise induction heating is suggested in [22]. A mathematical model and numerical equation solutions are presented based on the finite element method. The computer modelling results are illustrated graphically and analysed.

The article [23] concerns the determination of the influence of saturation and hysteresis on transient processes in electrical steel plates of electrical machines. The authors developed an analytical method for calculating the transient process in a magnetic circuit with an approximated magnetisation curve. One part of the magnetic circuit takes into account saturation, and the other is linear. The operational calculus was used to solve the Maxwell equations characterising the transient process. The applied method was verified based on the simulation results using an adapted part of the MATLAB/Simulink software complex. The results of computer simulation of transient electromagnetic processes in a transformer were presented and compared with the results obtained via other methods.

The Lagrangian and Eulerian–Lagrangian formalisms are used by many authors [24,25,26,27,28,29]. Molina-Santana, E. et al. [24] model a buck–boost DC–DC converter including power losses and a strongly nonlinear inductor. Alternative Hamiltonian and Lagrangian formalisms are used in [25,26], where in [25] they are related to relativistic mechanics. The Lagrangian and Hamiltonian formalisms for a hypothetical charged and uncharged particle are considered on the basis of the relativistic equation of dynamics and the integrals of the particle’s motion. The relativistic invariance law for the energy and momentum conservation is given in the Lorentz representation. To select a variety of generalised coordinates and momenta, the Lagrange equation of the second kind can be modified with regard to the relativist laws of energy and momentum conservation. In [26], the Lagrangian formalism is established for differential equations with special mathematical physics functions as solutions. The formalism is based on standard or non-standard Lagrangians. It was also proved that the latter require a modification of the Euler–Lagrange equation by auxiliary conditions, which is a new development of the calculus of variations.

The real form of transient processes is best represented by applying the mathematical modelling apparatus which includes the fundamental principles of applied physics, in particular Hamilton’s principle [30,31,32,33,34]. Under the classic interpretation, the principle does not apply to electrical power systems, since dissipation and non-potential external forces cannot be addressed as part of the conservative Lagrange’s function [30]. The Hamilton–Ostrogradsky’s principle is modified to solve this problem by expanding the Lagrange function [2].

Fiori, S. [31] gives a theoretical description of the variational approach to formulating fourth-order dynamic systems including differentiable manifolds based on the Hamilton–d’Alembert principle of analytical mechanics. The approach introduces the Lagrange function, dependent on the kinetic energy and the covariant acceleration energy, and the potential energy function that addresses conservative forces. Mavroeidis, C.P. et al. [32] provide an original application of the variational method. The differential equations governing the motion of volume fluids are derived by means of the standard methods of variational calculus.

Necessary and sufficient conditions for the validity of Hamilton’s variational principle are derived in [33] for circuits consisting of Chua’s periodical table. In linear systems, the motion equations include only even-order derivatives of independent variables. The Lagrange’s formalism known from classic mechanics can be used to describe phenomena across electric circuits.

This review of the available literature certainly suggests a very wide applicability of the electric, magnetic, and electromagnetic circuit theory. The methods presented can be applied to far more than just electrical machines and devices. This application of the modified Hamilton–Ostrogradsky’s principle shows an energetic approach to mathematical formulas.

2. The Concept of Modifying the Hamilton–Ostrogradsky Principle

In the classical interpretation, Hamilton’s principle is formulated as follows [2]: “Among all possible virtual motions of the holonomic system, only for the real motion the Hamiltonian action functional takes a stationary value over the entire integration interval” and is mathematically justified with the following formula:

$$S={\displaystyle \underset{t_1}{\overset{t_2}{\int }}{L}^{*}dt }, L^{*}=L^{*}(q_k,{\dot{q}}_k,t)$$

(1)

where S is the Hamiltonian action functional; L* is a non-conservative Lagrangian function of the second kind [1], $q_k,{\dot{q}}_k$ are generalised coordinates and velocities, k is the number of degrees of freedom, and t is time.

An important fact follows from the formulation of the above-mentioned principle and expression (1): the variation of the Hamiltonian action functional for real motion is always equal to zero [1,2].

$$\mathrm{\delta }S=\mathrm{\delta }{\displaystyle \underset{t_1}{\overset{t_2}{\int }}{L}^{*}dt=0} , L^{*}=T-P$$

(2)

where T is kinetic energy, and P is potential energy.

It is obvious that for mathematical modelling of real systems, including electrical systems, Hamilton’s principle in versions (1) and (2) cannot be used, because in the non-conservative Lagrangian, external and internal energy dissipation in the system is not taken into account, and the influence of external non-potential forces is not taken into account. It is precisely the reason indicated above that became the basis for modifying the well-known Hamilton’s principle.

One of the more successful attempts to modify Hamilton’s principle was presented in the monograph [1], where the authors proposed extending the non-conservative Lagrange’s function with two additional terms: the first term was added to the kinetic energy (coenergy) and took into account the dissipation of energy in the analysed object, and the second term was added to the potential energy and took into account the influence of external non-potential forces. It is worth noting that the aforementioned extension of the function by two terms was carried out empirically, i.e., without mathematical justification. Thanks to such an extension, the conservative Lagrange’s function has acquired a non-conservative meaning. The authors gave it the title of the non-conservative Lagrange’s function.

In contrast to the authors of [1], the author of [2] proposed to obtain the non-conservative Lagrange function exclusively in a mathematical way. For this purpose, the well-known Lagrange equations of the second kind were used. Although in both cases the analytical image of the non-conservative function has the same appearance, the results presented in [2] mathematically justify the correctness of the idea of formulating the extended (non-conservative) Lagrange function in [1]. The extended (non-conservative) Lagrange function is as follows [1,2]:

$$L=\tilde{T}-P+\mathrm{\Phi }-D , {\mathrm{\Phi }}^{*}={\displaystyle \underset{0}{\overset{t}{\int }}{\mathrm{\Phi }}_R{({\dot{q}}_k,t)}_{|t=\tau }d\tau }$$

(3)

where L is the extended non-conservative Lagrange function; $\tilde{T}$ is the kinetic co-energy; D is the work of external non-potential forces (in [1] called the non-conservative potential); Φ is the function of external and internal dissipation of energy (in [1] called the non-conservative kinetic energy (co-energy)); Φ_R is the dissipative function (usually the Rayleigh function); τ is an additional variable of integration.

Analysing Equations (1) and (2) and the work [1], it can be concluded that Hamilton’s principle can be applied only to lumped parameter systems. For distributed parameter systems, the above approach is not applicable. Therefore, the applicability of the well-known Hamilton’s principle was mathematically extended by the Ukrainian mathematician Mykhail Ostrogradsky to distributed parameter systems without changing the conservative Lagrange’s function (density of the Lagrange’s function). Since then, Hamilton’s principle is called the Hamilton–Ostrogradsky’s integral variational principle. Note that the idea of extending the classical (conservative) Lagrange function is also valid for distributed parameter systems. In such a case, the Lagrange’s function is transformed into the density of the nonconservative Lagrange’s function. In [2], the modified Hamilton–Ostrogradsky’s principle was used to analyse both lumped- and distributed-parameter systems.

3. Mathematical Model of the System

Electromagnetic circuits are analysed as holonomic lumped parameter systems with a finite number of degrees of freedom. This means the number of generalised coordinates is equal to the number of the electromagnetic system’s degrees of freedom. In this case, the generalised coordinates are as follows: $q\equiv {(q_1, q_2,\dots , q_n)}^T$; generalised velocities: $\dot{q}\equiv {({\dot{q}}_1, {\dot{q}}_2,\dots , {\dot{q}}_n)}^T$. The following are used in the case of generalised coordinates: the flux linkages for circuit elements and charges in the windings of each ferromagnetic core element: $q_k={\mathrm{\Phi }}_k, q_{k,j}=Q_{k,j}$; in the case of generalised velocities, ${\displaystyle \frac{d{q}_k}{dt}}\equiv {\dot{q}}_k={\displaystyle \frac{d{\mathrm{\Phi }}_k}{dt}}=-e_k, {\dot{q}}_{k,j}={\displaystyle \frac{d{Q}_{k,j}}{dt}}=i_{k,j}$ are self-induction EMFs (electromotive forces) generated via a magnetic circuit (without regard to leakage flux linkages) and currents across all the windings of each element, $k=1, 2,\dots , N$, N is the number of magnetic elements in a magnetic system, $j=1, 2,\dots , M_k$, and $M_k$ is the number of coils around the k-th magnetic system element.

A jth magnetic system circuit included in a general electromagnetic circuit will be discussed as an example. The diagram of an electric circuit including a lumped parameter electric winding is illustrated in Figure 1. For the purpose of analysing physical processes, the computational diagram of the electromagnetic device is shown in Figure 2.

The key to Figure 2 is as follows: V is a quantity analogous to the voltage drop equal to the product of reluctance and magnetic flux, ρ_1,2,3 are nonlinear reluctances of the ferromagnetic branches of the magnetic circuit, ρ₀ is the linear reluctance of the air gap, and F is the magnetomotive force.

Based on Ampère’s circuital law, the magnetomotive force MMF in the k-th magnetic device is formulated as follows [2]:

$$F_k={\displaystyle \sum _{j=1}^{M_k}{w}_{k,j}{i}_{k,j}} , j=1,\dots ,M_k$$

(4)

where M_k is the number of windings across the k-th circuit element, and $w_{k,j}$ is the number of turns of the j-th winding on the k-th circuit element.

The expanded action functional will be expressed as per Hamilton–Ostrogradsky (see Equations (2) and (3)).

$$S={\displaystyle \underset{t_1}{\overset{t_2}{\int }}L}dt={\displaystyle \underset{t_1}{\overset{t_2}{\int }}\left(\tilde{T}-P+\sout{\mathrm{\Phi }}-D\right)}dt , L=\tilde{T}-P+\sout{\mathrm{\Phi }}-D,$$

(5)

where L is the expanded Lagrange function, $\tilde{T}$ is magnetic (co)energy, P is potential energy, $\sout{\mathrm{\Phi }}$ is dissipation energy, and D is the work of external non-potential forces.

To derive the electromagnetic state equations, the functional variation S, cf. (5), is compared to zero. In this way, the Euler–Lagrange’s equation can finally be generated, which describes the electromagnetic circuit’s mathematical model [2] (see Equations (2) and (3)):

$$\mathrm{\delta }S=\mathrm{\delta }{\displaystyle \underset{t_1}{\overset{t_2}{\int }}L}dt={\displaystyle \underset{t_1}{\overset{t_2}{\int }}\mathrm{\delta }\left(\tilde{T}-P+\sout{\mathrm{\Phi }}-D\right)}dt=0.$$

(6)

where δ is variation symbol.

The Euler–Lagrange equation is arrived at on defining the functional variation (5):

$${\displaystyle \frac{d}{dt}}{\displaystyle \frac{\partial L}{\partial {\dot{q}}_k}}-{\displaystyle \frac{\partial L}{\partial q_k}}=0$$

(7)

All the elements of the expanded Lagrangian will be expressed to define the electromagnetic circuit’s mathematical model:

$$\tilde{T}={\displaystyle \sum _{k=1}^N{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{i_{k,j}}{\int }}{\mathrm{\Psi }}_{k,j}d{i}_{k,j}}}}, P={\displaystyle \sum _{k=1}^N{\displaystyle \underset{0}{\overset{{\mathrm{\Phi }}_k}{\int }}\left({\mathrm{\rho }}_{0k}+{\mathrm{\rho }}_k\right){\mathrm{\Phi }}_{k}d{\mathrm{\Phi }}_k}};$$

(8)

$$\mathrm{\Phi }={\displaystyle \sum _{k=1}^N{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{t}{\int }}\frac{r_{k,j}^2{i}_{k,j}}{2}d\tau }}}, D={\displaystyle \sum _{k=1}^N{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{t}{\int }}{\left(u_{k,j}{i}_{k,j}\right)}_{|t=\tau }}d\tau }}.$$

(9)

The expanded Lagrange function will be as follow:

$$L={\displaystyle \sum _{k=1}^N[{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{i_{k,j}}{\int }}{\mathrm{\Psi }}_{k,j}d{i}_{k,j}}}}-{\displaystyle \underset{0}{\overset{{\mathrm{\Phi }}_k}{\int }}\left({\mathrm{\rho }}_{0k}+{\mathrm{\rho }}_k\right){\mathrm{\Phi }}_{k}d{\mathrm{\Phi }}_k}+{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{t}{\int }}\frac{r_{k,j}^2{i}_{k,j}}{2}}d\tau }-{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{t}{\int }}{\left(u_{k,j}{i}_{k,j}\right)}_{|t=\tau }}d\tau }].$$

(10)

Additionally, the reasons for using the theorem on the function differentiation change for cyclical motions are provided. In Helmholtz’s terminology, cyclical motions are motions for which energy dissipation depends only on the generalised coordinates: $f_k=f({\dot{q}}_k,t)$. This assumption holds true in our case.

Theorem 1.

The following relationship applies to cyclical motions (the function dependent only on generalised velocities and time: $f({\dot{q}}_k,t)$):

$${\displaystyle \frac{\partial }{\partial {\dot{q}}_k}}{\displaystyle \frac{df({\dot{q}}_k,t)}{dt}}={\displaystyle \frac{d}{dt}}{\displaystyle \frac{\partial f({\dot{q}}_k,t)}{\partial {\dot{q}}_k}}.$$

(11)

Proof of Theorem 1.

In the general case, the equation for non-cyclical motions is considered, namely $f=f(q_k,{\dot{q}}_k,t)$. Assuming k = 1, the full derivative of the composite function is derived according to the definition [2]:

$${\displaystyle \frac{df(q(t),\dot{q}(t),t)}{dt}}={\displaystyle \frac{\partial f(q,\dot{q},t)}{\partial q}}{\displaystyle \frac{dq}{dt}}+{\displaystyle \frac{\partial f(q,\dot{q},t)}{\partial \dot{q}}}{\displaystyle \frac{d\dot{q}}{dt}}+{\displaystyle \frac{\partial f(q,\dot{q},t)}{\partial t}}={\displaystyle \frac{\partial f}{\partial q}}\dot{q}+{\displaystyle \frac{\partial f}{\partial \dot{q}}}\ddot{q}+{\displaystyle \frac{\partial f}{\partial t}}.$$

(12)

Following from (12), the full total derivative of $\mathrm{\phi }={\displaystyle \frac{\partial f(q(t),\dot{q}(t),t)}{\partial \dot{q}(t)}}$ will be calculated:

$${\displaystyle \frac{d}{dt}}{\displaystyle \frac{\partial f(q,\dot{q},t)}{\partial \dot{q}}}={\displaystyle \frac{{\partial }^{2}f}{\partial q\partial \dot{q}}}{\displaystyle \frac{dq}{dt}}+{\displaystyle \frac{{\partial }^{2}f}{\partial {\dot{q}}^2}}{\displaystyle \frac{d\dot{q}}{dt}}+{\displaystyle \frac{{\partial }^{2}f}{\partial t\partial \dot{q}}}={\displaystyle \frac{{\partial }^{2}f}{\partial q\partial \dot{q}}}\dot{q}+{\displaystyle \frac{{\partial }^{2}f}{\partial {\dot{q}}^2}}\ddot{q}+{\displaystyle \frac{{\partial }^{2}f}{\partial t\partial \dot{q}}}.$$

(13)

The Function (13) is converted into the following:

$$\begin{gathered} {\displaystyle \frac{d}{dt}}{\displaystyle \frac{\partial f(q,\dot{q},t)}{\partial \dot{q}}}={\displaystyle \frac{{\partial }^{2}f}{\partial \dot{q}\partial q}}\dot{q}+{\displaystyle \frac{{\partial }^{2}f}{\partial {\dot{q}}^2}}\ddot{q}+{\displaystyle \frac{{\partial }^{2}f}{\partial \dot{q}\partial t}}={\displaystyle \frac{\partial }{\partial \dot{q}}}\left({\displaystyle \frac{\partial f}{\partial q}}{\displaystyle \frac{dq}{dt}}+{\displaystyle \frac{\partial f}{\partial \dot{q}}}{\displaystyle \frac{d\dot{q}}{dt}}+{\displaystyle \frac{\partial f}{\partial t}}\right)-{\displaystyle \frac{\partial f}{\partial q}}= \\ {\displaystyle \frac{{\partial }^{2}f}{\partial \dot{q}\partial q}}{\displaystyle \frac{dq}{dt}}+{\displaystyle \frac{\partial f}{\partial q}}{\displaystyle \frac{\partial \dot{q}}{\partial \dot{q}}}+{\displaystyle \frac{{\partial }^{2}f}{\partial {\dot{q}}^2}}{\displaystyle \frac{d\dot{q}}{dt}}+{\displaystyle \frac{\partial f}{\partial \dot{q}}}{\displaystyle \frac{\partial \ddot{q}}{\partial \dot{q}}}+{\displaystyle \frac{{\partial }^{2}f}{\partial \dot{q}\partial t}}-{\displaystyle \frac{\partial f}{\partial q}}= \\ {\displaystyle \frac{{\partial }^{2}f}{\partial \dot{q}\partial q}}{\displaystyle \frac{dq}{dt}}+{\displaystyle \frac{\partial f}{\partial q}}\cdot 1+{\displaystyle \frac{{\partial }^{2}f}{\partial {\dot{q}}^2}}{\displaystyle \frac{d\dot{q}}{dt}}+{\displaystyle \frac{\partial f}{\partial \dot{q}}}\cdot 0+{\displaystyle \frac{{\partial }^{2}f}{\partial \dot{q}\partial t}}-{\displaystyle \frac{\partial f}{\partial q}}={\displaystyle \frac{{\partial }^{2}f}{\partial \dot{q}\partial q}}\dot{q}+{\displaystyle \frac{{\partial }^{2}f}{\partial {\dot{q}}^2}}\ddot{q}+{\displaystyle \frac{{\partial }^{2}f}{\partial \dot{q}\partial t}}={\displaystyle \frac{d}{dt}}{\displaystyle \frac{\partial f(q,\dot{q},t)}{\partial \dot{q}}}. \end{gathered}$$

(14)

The bracketed expression in (12) is the total derivative of $f(q,\dot{q},t)$, see (9).

We will then write the following:

$${\displaystyle \frac{d}{dt}}{\displaystyle \frac{\partial f(q,\dot{q},t)}{\partial \dot{q}}}={\displaystyle \frac{\partial }{\partial \dot{q}}}{\displaystyle \frac{df(q,\dot{q},t)}{dt}}-{\displaystyle \frac{df(q,\dot{q},t)}{\partial q}}.$$

(15)

Since the motions are cyclical, $f=f(\dot{q},t)$, and the second term on the right-hand side of (15) will be equal to zero. □

If k = 1, …, N, the relationship (11) will of course be true, too. The same holds for a changed integration sequence:

Theorem 2.

$${\displaystyle \frac{d}{dt}}{\displaystyle \underset{0}{\overset{\dot{q}}{\int }}f(\dot{q},t)d\dot{q}}={\displaystyle \underset{0}{\overset{\dot{q}}{\int }}\frac{df(\dot{q},t)}{dt}d\dot{q}}.$$

(16)

Proof of Theorem 2.

The left and right sides of (16) will be differentiated against the coordinate $\dot{q}(t)$ considering the theorem on the derivative of the integral relative to its upper boundary, see (11) for cyclical motions:

$${\displaystyle \frac{\partial }{\partial \dot{q}}}{\displaystyle \frac{d}{dt}}{\displaystyle \underset{0}{\overset{\dot{q}}{\int }}f(\dot{q},t)d\dot{q}}={\displaystyle \frac{d}{dt}}{\displaystyle \frac{\partial }{\partial \dot{q}}}{\displaystyle \underset{0}{\overset{\dot{q}}{\int }}f(\dot{q},t)d\dot{q}}={\displaystyle \frac{df(\dot{q},t)}{dt}};$$

(17)

$${\displaystyle \frac{\partial }{\partial \dot{q}}}{\displaystyle \underset{0}{\overset{\dot{q}}{\int }}\frac{df(\dot{q},t)}{dt}d\dot{q}}={\displaystyle \frac{df(\dot{q},t)}{dt}}.$$

(18)

The results are identical. The differentiation addresses boundary conditions.

The expanded Lagrange function will be substituted into the Euler–Lagrange equation considering (11) and (17):

$$\begin{gathered} {\displaystyle \frac{d}{dt}}{\displaystyle \frac{\partial }{\partial i_{k,j}}}{\displaystyle \sum _{k=1}^N[{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{i_{k,j}}{\int }}{\mathrm{\Psi }}_{k,j}d{i}_{k,j}}}}-{\displaystyle \underset{0}{\overset{{\mathrm{\Phi }}_k}{\int }}\left({\mathrm{\rho }}_{0k}+{\mathrm{\rho }}_k\right){\mathrm{\Phi }}_{k}d{\mathrm{\Phi }}_k}+{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{t}{\int }}\frac{r_{k,j}^2{i}_{k,j}}{2}}d\tau }-{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{t}{\int }}{\left(u_{k,j}{i}_{k,j}\right)}_{|t=\tau }}d\tau }]- \\ {\displaystyle \frac{\partial }{\partial Q_{k,j}}}{\displaystyle \sum _{k=1}^N[{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{i_{k,j}}{\int }}{\mathrm{\Psi }}_{k,j}d{i}_{k,j}}}}-{\displaystyle \underset{0}{\overset{{\mathrm{\Phi }}_k}{\int }}\left({\mathrm{\rho }}_{0k}+{\mathrm{\rho }}_k\right){\mathrm{\Phi }}_{k}d{\mathrm{\Phi }}_k}+{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{t}{\int }}\frac{r_{k,j}^2{i}_{k,j}}{2}}d\tau }-{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{t}{\int }}{\left(u_{k,j}{i}_{k,j}\right)}_{|t=\tau }}d\tau }]= \\ {\displaystyle \frac{d}{dt}}{\displaystyle \frac{\partial }{\partial i_{k,j}}}{\displaystyle \sum _{k=1}^N{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{i_{k,j}}{\int }}{\mathrm{\Psi }}_{k,j}d{i}_{k,j}}}+}{\displaystyle \frac{\partial }{\partial i_{k,j}}}{\displaystyle \frac{d}{dt}}{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{t}{\int }}(\frac{r_{k,j}^2{i}_{k,j}}{2}}}-{\left(u_{k,j}{i}_{k,j}\right))}_{|t=\tau }d\tau -0= \\ {\displaystyle \frac{d{\mathrm{\Psi }}_{k,j}}{dt}}+r_{k,j}{i}_{k,j}-u_{k,j}=0\Rightarrow {\displaystyle \frac{d{\mathrm{\Psi }}_{k,j}}{dt}}=u_{k,j}-r_{k,j}{i}_{k,j}. \end{gathered}$$

(19)

We assume that the total flux linkages are the sum of all leakage and main (operational) flux linkages. They will be expressed as follows:

$${\mathrm{\Psi }}_{k,j}=L_{\mathrm{\sigma }k,j}{i}_{k,j}+{\mathrm{\psi }}_{k,j},$$

(20)

where the main fluxes linkages are formulated as follows:

$${\mathrm{\psi }}_{k,j}=w_{k,j}{\mathrm{\Phi }}_k.$$

(21)

In the event, the Euler–Lagrange’s equations will be written as follows, considering (11) and (17):

$$\begin{gathered} {\displaystyle \frac{d}{dt}}{\displaystyle \frac{\partial }{\partial (-e_k)}}{\displaystyle \sum _{k=1}^N[{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{i_{k,j}}{\int }}\left(L_{\mathrm{\sigma }k,j}{i}_{k,j}+w_{k,j}{\mathrm{\Phi }}_k\right)d{i}_{k,j}}}}-{\displaystyle \underset{0}{\overset{{\mathrm{\Phi }}_k}{\int }}\left({\mathrm{\rho }}_{0k}+{\mathrm{\rho }}_k\right){\mathrm{\Phi }}_{k}d{\mathrm{\Phi }}_k}+{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{t}{\int }}\frac{r_{k,j}^2{i}_{k,j}}{2}}d\tau }-{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{t}{\int }}{\left(u_{k,j}{i}_{k,j}\right)}_{|t=\tau }}d\tau }]- \\ {\displaystyle \frac{\partial }{\partial ({\mathrm{\Phi }}_k)}}{\displaystyle \sum _{k=1}^N[{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{i_{k,j}}{\int }}\left(L_{\mathrm{\sigma }k,j}{i}_{k,j}+w_{k,j}{\mathrm{\Phi }}_k\right)d{i}_{k,j}}}}-{\displaystyle \underset{0}{\overset{{\mathrm{\Phi }}_k}{\int }}\left({\mathrm{\rho }}_{0k}+{\mathrm{\rho }}_k\right){\mathrm{\Phi }}_{k}d{\mathrm{\Phi }}_k}+{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{t}{\int }}\frac{r_{k,j}^2{i}_{k,j}}{2}}d\tau }-{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{t}{\int }}{\left(u_{k,j}{i}_{k,j}\right)}_{|t=\tau }}d\tau }]= \\ {\displaystyle \frac{\partial }{\partial (-e_k)}}{\displaystyle \sum _{k=1}^N{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{i_{k,j}}{\int }}{w}_{k,j}\frac{d{\mathrm{\Phi }}_k}{dt}d{i}_{k,j}}}-}{\displaystyle \frac{\partial }{\partial ({\mathrm{\Phi }}_k)}}{\displaystyle \sum _{k=1}^N\left(-{\displaystyle \underset{0}{\overset{{\mathrm{\Phi }}_k}{\int }}\left({\mathrm{\rho }}_{0k}+{\mathrm{\rho }}_k\right){\mathrm{\Phi }}_{k}d{\mathrm{\Phi }}_k}\right)}= \\ {\displaystyle \frac{\partial }{\partial (-e_k)}}{\displaystyle \sum _{k=1}^N{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \underset{0}{\overset{i_{k,j}}{\int }}{w}_{k,j}(-e_k)d{i}_{k,j}}}}-\left(-\left({\mathrm{\rho }}_{0k}+{\mathrm{\rho }}_k\right){\mathrm{\Phi }}_k\right)={\displaystyle \sum _{j=1}^{M_k}{w}_{k,j}{i}_{k,j}}+\left({\mathrm{\rho }}_{0k}+{\mathrm{\rho }}_k\right){\mathrm{\Phi }}_k=0. \end{gathered}$$

(22)

The Equation (22) will be converted into the following:

$$\left({\mathrm{\rho }}_{0k}+{\mathrm{\rho }}_k\right){\mathrm{\Phi }}_k+{\displaystyle \sum _{j=1}^{M_k}{w}_{k,j}{i}_{k,j}}=V_k+F_k=0,$$

(23)

where

$$V_k=({\mathrm{\rho }}_{0k}+{\mathrm{\rho }}_k){\mathrm{\Phi }}_k, F_k={\displaystyle \sum _{j=1}^{M_k}{w}_{k,j}{i}_{k,j}},$$

(24)

where $V_k$ is the product of reluctance and magnetic flux (voltage-drop-like quantity) (calculation value), and $F_k$ is the magnetomotive force (MMF).

Equations (22) and (24) will be expressed in matrix-vectoral format:

$$({{\rm P}}_0+{\rm P}){\mathbf{\Phi }}_K+F=0 \Rightarrow V+F=0;$$

(25)

$${\displaystyle \frac{d{\mathbf{\Psi }}_j}{dt}}=u_j-r_j{i}_j, {\mathbf{\Psi }}_j=L_{\mathrm{\sigma }j}{i}_j+w_j\mathbf{\Phi },$$

(26)

where

$${{\rm P}}_0\equiv \mathrm{diag}({\mathrm{\rho }}_{01}, {\mathrm{\rho }}_{02}, \dots , {\mathrm{\rho }}_{0n}), {\rm P}\equiv \mathrm{diag}({\mathrm{\rho }}_1, {\mathrm{\rho }}_2, \dots , {\mathrm{\rho }}_n) , F\equiv {\left(w_j{i}_j\right)}^{\mathrm{T}},$$

(27)

$$w_j\equiv (w_1^{(k)}, w_2^{(k)}, \dots , w_{M_k}^{(k)}) , i_j\equiv (i_1^{(k)}, i_2^{(k)}, \dots , i_{M_k}^{(k)}),$$

(28)

where M_k is the number of turns of the jth electric circuit (winding).

The system of Equations (25)–(28) implies that topological methods of analysing magnetic circuits are needed to calculate the real magnetic fluxes.

The system (25)–(28) contains the matrix-vectoral equations of a nonlinear magnetic circuit that are solved jointly. □

4. Practical Implementation of the Mathematical Model

Equations (25)–(28) take into account the geometrical dimensions of the magnetic core according to Figure 3, with the core thickness h being 0.15 m. The parameters of the electromagnetic circuit are as follows: the coil supply voltage has the same analytical form in all cases:

$u(t)=U_m\sin (\omega t+\mathrm{\psi })$ V, $U_{m1}=300$ V, $U_{m3}^{\langle 1\rangle }=200$ V, $U_{m3}^{\langle 2\rangle }=250$ V, $w_1=100$, $w_3^{\langle 1\rangle }=200$, $w_3^{\langle 2\rangle }=150$, $L_{\mathrm{\sigma }1}=40$ mH, $L_{\mathrm{\sigma }3}^{\langle 1\rangle }=30$ mH, $L_{\mathrm{\sigma }3}^{\langle 2\rangle }=40$ mH, $r_1=1$ $\Omega$, $r_3^{\langle 1\rangle }=1.5$ $\Omega$, $r_3^{\langle 2\rangle }=2$ $\Omega$, ${\mathrm{\psi }}_1={\displaystyle \raisebox{1ex}{$\mathrm{\pi }$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$, ${\mathrm{\psi }}_3^{\langle 1\rangle }={\displaystyle \raisebox{1ex}{$\mathrm{\pi }$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}$, ${\mathrm{\psi }}_3^{\langle 2\rangle }=-{\displaystyle \raisebox{1ex}{$\mathrm{\pi }$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}$. The steel magnetisation curve: $H(B)=\mathrm{\nu }(B)B=65B+270B^5$.

Based on the second Maxwell’s equation $\left(\nabla \times E=-{\displaystyle \frac{\partial B}{\partial t}}\right)$, the magnetic induction vector B is utilised in the electromagnetic field equations. The finite difference method can serve to solve the second Maxwell’s equation. In the event, the equation related to B is presented in the normal Cauchy format, integrated over time with explicit or implicit numerical methods. Magnetic induction vector B, not the magnetic field strength vector H, is used for the purposes of the calculations. The coefficient of reverse magnetic permeability $\mathrm{\nu }=1/\mathrm{\mu }$ is introduced to this end instead of the usual coefficient of magnetic permeability $\mathrm{\mu }$. This has one more advantage: the field equations including B and $\mathrm{\nu }$ are less rigid than the usual form, which involves H and $\mathrm{\mu }$. The traditional magnetisation curve $B(H)$ will thus be converted into a reverse magnetisation curve $H(B)$, symmetrical to the line $y=x$ [2].

The equations for the magnetic part of the electromagnetic device in Figure 3 and Figure 4 will be expressed in the matrix-vectoral format considering circuit directions and Equation (25):

$$F=({{\rm P}}_0+{\rm P}){\mathbf{\Phi }}_K.$$

(29)

The coordinate format of (29), cf. 4., is as follows:

$$\left[\begin{array}{c}{w}_1{i}_1\\ w_3^{\langle 1\rangle }{i}_3^{\langle 1\rangle }-w_3^{\langle 2\rangle }{i}_3^{\langle 2\rangle }\\ 0\end{array}\right]=\left[\begin{array}{ccc}{\mathrm{\rho }}_{01}+{\mathrm{\rho }}_1(B_1)+{\mathrm{\rho }}_2(B_2)& -{\mathrm{\rho }}_2(B_2)& -{\mathrm{\rho }}_2(B_2)\\ -{\mathrm{\rho }}_2(B_2)& {\mathrm{\rho }}_3(B_3)+{\mathrm{\rho }}_2(B_2)& {\mathrm{\rho }}_2(B_2)\\ -{\mathrm{\rho }}_2(B_2)& {\mathrm{\rho }}_2(B_2)& {\mathrm{\rho }}_2(B_2)+{\mathrm{\rho }}_0\end{array}\right]\left[\begin{array}{c}{\mathrm{\Phi }}_{k1}\\ {\mathrm{\Phi }}_{k2}\\ {\mathrm{\Phi }}_{k3}\end{array}\right] , \left[\begin{array}{c}{\mathrm{\Phi }}_1\\ {\mathrm{\Phi }}_2\\ {\mathrm{\Phi }}_3\\ {\mathrm{\Phi }}_0\end{array}\right]=\left[\begin{array}{ccc}1& & \\ -1& 1& 1\\ & 1& \\ & & 1\end{array}\right]\left[\begin{array}{c}{\mathrm{\Phi }}_{k1}\\ {\mathrm{\Phi }}_{k2}\\ {\mathrm{\Phi }}_{k3}\end{array}\right],$$

(30)

where ${\mathrm{\Phi }}_j$ are real flux linkages in the system j = 1, 2, 3, $\langle 1\rangle , \langle 2\rangle$ represents the first and second coils in the third magnetic circuit branch, see Figure 3 and Figure 4, ${\mathrm{\Phi }}_0$ is the magnetic flux non-compensated in a core centre, and ${\mathrm{\rho }}_0$ is the calculation air gap reluctance associated with ${\mathrm{\Phi }}_0$.

Equations for the electric part of the device will be expressed classically, cf. (26):

$${\displaystyle \frac{d{\mathrm{\Psi }}_1}{dt}}=u_1-r_1{i}_1 , {\displaystyle \frac{d{\mathrm{\Psi }}_3^{\langle 1\rangle }}{dt}}=u_3^{\langle 1\rangle }-r_3^{\langle 1\rangle }{i}_3^{\langle 1\rangle } , {\displaystyle \frac{d{\mathrm{\Psi }}_3^{\langle 2\rangle }}{dt}}=u_3^{\langle 2\rangle }-r_3^{\langle 2\rangle }{i}_3^{\langle 2\rangle },$$

(31)

$${\mathrm{\Psi }}_1=L_{\mathrm{\sigma }1}{i}_1+w_1{\mathrm{\Phi }}_1 , {\mathrm{\Psi }}_3^{\langle 1\rangle }=L_{\mathrm{\sigma }3}^{\langle 1\rangle }{i}_3^{\langle 1\rangle }+w_3^{\langle 1\rangle }{\mathrm{\Phi }}_3 , {\mathrm{\Psi }}_3^{\langle 2\rangle }=L_{\mathrm{\sigma }3}^{\langle 2\rangle }{i}_3^{\langle 2\rangle }+w_3^{\langle 2\rangle }{\mathrm{\Phi }}_3.$$

(32)

The reluctances of the core parts are computed classically:

$$\mathrm{\rho }(B)={\displaystyle \frac{\mathrm{\nu }(B)l}{S}}, {\mathrm{\rho }}_0={\displaystyle \frac{l}{{\mathrm{\mu }}_{0}S}}={\displaystyle \frac{{\mathrm{\nu }}_{0}l}{S}} , \mathrm{\nu }(B)={\displaystyle \frac{1}{\mathrm{\mu }(B)}}={\displaystyle \frac{H(B)}{B}} , {\mathrm{\nu }}^{\partial }(B)={\displaystyle \frac{\partial H(B)}{\partial B}}.$$

(33)

The magnetic induction of the core limbs will be calculated based on the definitions, see Figure 3.

$$\mathrm{\Phi }={\displaystyle \underset{S}{\iint }B\cdot dS} \Rightarrow {\mathrm{\Phi }}_j=B_j{S}_j\Rightarrow B_j={\displaystyle \frac{{\mathrm{\Phi }}_j}{S_j}},$$

(34)

where S_j is the magnetic core cross-sectional area.

Equation (30) become as follows in coordinate format:

$$\begin{gathered} w_1{i}_1=\left({\mathrm{\rho }}_{01}+{\mathrm{\rho }}_1(B_1)+{\mathrm{\rho }}_2(B_2)\right){\mathrm{\Phi }}_{k1}-{\mathrm{\rho }}_2(B_2){\mathrm{\Phi }}_{k2}-{\mathrm{\rho }}_2(B_2){\mathrm{\Phi }}_{k3}=\left({\mathrm{\rho }}_{01}+{\mathrm{\rho }}_1(B_1)\right){\mathrm{\Phi }}_{k1}- \\ {\mathrm{\rho }}_2(B_2)(-{\mathrm{\Phi }}_{k1}+{\mathrm{\Phi }}_{k2}+{\mathrm{\Phi }}_{k3})=\left({\mathrm{\rho }}_{01}+{\mathrm{\rho }}_1(B_1)\right){\mathrm{\Phi }}_1-{\mathrm{\rho }}_2(B_2){\mathrm{\Phi }}_2; \end{gathered}$$

(35)

$$\begin{gathered} w_3^{\langle 1\rangle }{i}_3^{\langle 1\rangle }-w_3^{\langle 2\rangle }{i}_3^{\langle 2\rangle }=-{\mathrm{\rho }}_2(B_2){\mathrm{\Phi }}_{k1}+({\mathrm{\rho }}_3(B_3)+{\mathrm{\rho }}_2(B_2)){\mathrm{\Phi }}_{k2}+{\mathrm{\rho }}_2(B_2){\mathrm{\Phi }}_{k3}={\mathrm{\rho }}_3(B_3){\mathrm{\Phi }}_{k2}+ \\ {\mathrm{\rho }}_2(B_2)(-{\mathrm{\Phi }}_{k1}+{\mathrm{\Phi }}_{k2}+{\mathrm{\Phi }}_{k3})={\mathrm{\rho }}_2(B_2){\mathrm{\Phi }}_2+{\mathrm{\rho }}_3(B_3){\mathrm{\Phi }}_3; \end{gathered}$$

(36)

$$\begin{gathered} 0=-{\mathrm{\rho }}_2(B_2){\mathrm{\Phi }}_{k1}+{\mathrm{\rho }}_2(B_2){\mathrm{\Phi }}_{k2}+\left({\mathrm{\rho }}_0+{\mathrm{\rho }}_2(B_2)\right){\mathrm{\Phi }}_{k3}={\mathrm{\rho }}_0{\mathrm{\Phi }}_{k3}+{\mathrm{\rho }}_2(B_2)(-{\mathrm{\Phi }}_{k1}+{\mathrm{\Phi }}_{k2}+{\mathrm{\Phi }}_{k3})= \\ {\mathrm{\rho }}_0{\mathrm{\Phi }}_0+{\mathrm{\rho }}_2(B_2){\mathrm{\Phi }}_2\Rightarrow {\mathrm{\Phi }}_0=-{\displaystyle \frac{{\mathrm{\rho }}_2(B_2)}{{\mathrm{\rho }}_0}}{\mathrm{\Phi }}_2. \end{gathered}$$

(37)

Differentiating (32) with respect to time and considering the initial conditions will result in the following:

$${\displaystyle \frac{d{\mathrm{\Psi }}_1}{dt}}=L_{\mathrm{\sigma }1}{\displaystyle \frac{d{i}_1}{dt}}+w_1{\displaystyle \frac{d{\mathrm{\Phi }}_1}{dt}} , {\displaystyle \frac{d{\mathrm{\Psi }}_3^{\langle 1\rangle }}{dt}}=L_{\mathrm{\sigma }3}^{\langle 1\rangle }{\displaystyle \frac{d{i}_3^{\langle 1\rangle }}{dt}}+w_3^{\langle 1\rangle }{\displaystyle \frac{d{\mathrm{\Phi }}_3}{dt}} , {\displaystyle \frac{d{\mathrm{\Psi }}_3^{\langle 2\rangle }}{dt}}=L_{\mathrm{\sigma }3}^{\langle 2\rangle }{\displaystyle \frac{d{i}_3^{\langle 2\rangle }}{dt}}+w_3^{\langle 2\rangle }{\displaystyle \frac{d{\mathrm{\Phi }}_3}{dt}}.$$

(38)

Differentiating (35) and (36) with respect to time and considering the initial conditions will result in the following:

$$w_1{\displaystyle \frac{d{i}_1}{dt}}={\displaystyle \frac{d}{dt}}\left[\left({\mathrm{\rho }}_{01}+{\mathrm{\rho }}_1(B_1)\right){\mathrm{\Phi }}_1\right]-{\displaystyle \frac{d}{dt}}\left({\mathrm{\rho }}_2(B_2){\mathrm{\Phi }}_2\right)={\displaystyle \frac{d{\mathrm{\rho }}_1(B_1)}{dt}}{\mathrm{\Phi }}_1+\left({\mathrm{\rho }}_{01}+{\mathrm{\rho }}_1(B_1)\right){\displaystyle \frac{d{\mathrm{\Phi }}_1}{dt}}-\left({\displaystyle \frac{d{\mathrm{\rho }}_2(B_2)}{dt}}{\mathrm{\Phi }}_2+{\mathrm{\rho }}_2(B_2){\displaystyle \frac{d{\mathrm{\Phi }}_2}{dt}}\right),$$

(39)

$$w_3^{\langle 1\rangle }{\displaystyle \frac{d{i}_3^{\langle 1\rangle }}{dt}}-w_3^{\langle 2\rangle }{\displaystyle \frac{d{i}_3^{\langle 2\rangle }}{dt}}={\displaystyle \frac{d{\mathrm{\rho }}_2(B_2)}{dt}}{\mathrm{\Phi }}_2+{\mathrm{\rho }}_2(B_3){\displaystyle \frac{d{\mathrm{\Phi }}_2}{dt}}+{\displaystyle \frac{d{\mathrm{\rho }}_3(B_3)}{dt}}{\mathrm{\Phi }}_3+{\mathrm{\rho }}_3(B_3){\displaystyle \frac{d{\mathrm{\Phi }}_3}{dt}},$$

(40)

with the total derivatives of nonlinear reluctances defined as follows:

$${\displaystyle \frac{d{\mathrm{\rho }}_j(B_j)}{dt}}={\displaystyle \frac{\partial {\mathrm{\rho }}_j(B_j)}{\partial B_j}}{\displaystyle \frac{d{B}_j}{dt}}={\displaystyle \frac{1}{S_j}}{\displaystyle \frac{l_1}{S_j}}{\displaystyle \frac{\partial {\mathrm{\nu }}_j(B_j)}{\partial B_j}}{\displaystyle \frac{d{\mathrm{\Phi }}_j}{dt}}={\displaystyle \frac{l_j}{S_j^2}}{\mathrm{\nu }}_j^{\partial }(B_j){\displaystyle \frac{d{\mathrm{\Phi }}_j}{dt}} , j=1, 2, 3;$$

(41)

Considering (41), we arrive at the following:

$$w_1{\displaystyle \frac{d{i}_1}{dt}}={\displaystyle \frac{l_1}{S_1^2}}{\mathrm{\nu }}_1^{\partial }(B_1){\displaystyle \frac{d{\mathrm{\Phi }}_1}{dt}}{\mathrm{\Phi }}_1+\left({\mathrm{\rho }}_{01}+{\mathrm{\rho }}_1(B_1)\right){\displaystyle \frac{d{\mathrm{\Phi }}_1}{dt}}-{\displaystyle \frac{l_2}{S_2^2}}{\mathrm{\nu }}_2^{\partial }(B_2){\displaystyle \frac{d{\mathrm{\Phi }}_2}{dt}}{\mathrm{\Phi }}_2-{\mathrm{\rho }}_2(B_2){\displaystyle \frac{d{\mathrm{\Phi }}_2}{dt}};$$

(42)

$$w_3^{\langle 1\rangle }{\displaystyle \frac{d{i}_3^{\langle 1\rangle }}{dt}}-w_3^{\langle 2\rangle }{\displaystyle \frac{d{i}_3^{\langle 2\rangle }}{dt}}={\displaystyle \frac{l_2}{S_2^2}}{\mathrm{\nu }}_2^{\partial }(B_2){\displaystyle \frac{d{\mathrm{\Phi }}_2}{dt}}{\mathrm{\Phi }}_2+{\mathrm{\rho }}_2(B_2){\displaystyle \frac{d{\mathrm{\Phi }}_2}{dt}}+{\displaystyle \frac{l_3}{S_3^2}}{\mathrm{\nu }}_3^{\partial }(B_3){\displaystyle \frac{d{\mathrm{\Phi }}_3}{dt}}{\mathrm{\Phi }}_3+{\mathrm{\rho }}_3(B_3){\displaystyle \frac{d{\mathrm{\Phi }}_3}{dt}}.$$

(43)

Taking (31) and (32) into account, we formulate the following:

$$L_{\mathrm{\sigma }1}{\displaystyle \frac{d{i}_1}{dt}}+w_1{\displaystyle \frac{d{\mathrm{\Phi }}_1 }{dt}}=u_1-r_1{i}_1;$$

(44)

$$L_{\mathrm{\sigma }3}^{\langle 1\rangle }{\displaystyle \frac{d{i}_3^{\langle 1\rangle }}{dt}}+w_3^{\langle 1\rangle }{\displaystyle \frac{d{\mathrm{\Phi }}_3}{dt}}=u_3^{\langle 1\rangle }-r_3^{\langle 1\rangle }{i}_3^{\langle 1\rangle };$$

(45)

$$L_{\mathrm{\sigma }3}^{\langle 2\rangle }{\displaystyle \frac{d{i}_3^{\langle 2\rangle }}{dt}}+w_3^{\langle 2\rangle }{\displaystyle \frac{d{\mathrm{\Phi }}_3}{dt}}=u_3^{\langle 2\rangle }-r_3^{\langle 2\rangle }{i}_3^{\langle 2\rangle }.$$

(46)

Applying Kirchhoff’s first law to for a magnetic circuit node in consideration of (37) and Figure 4 leads to the following:

$${\mathrm{\Phi }}_1+{\mathrm{\Phi }}_2-{\mathrm{\Phi }}_3+{\mathrm{\Phi }}_0=0 \Rightarrow {\mathrm{\Phi }}_1+{\displaystyle \frac{1}{{\mathrm{\rho }}_0}}\left({\mathrm{\rho }}_0+{\mathrm{\rho }}_2(B_2)\right){\mathrm{\Phi }}_2-{\mathrm{\Phi }}_3=0.$$

(47)

If the magnetic fluxes are not compensated, i.e., ${\mathrm{\rho }}_0\to \infty$, then ${\mathrm{\Phi }}_0\to 0$ and (47) is simplified as follows:

$${\mathrm{\Phi }}_1+{\mathrm{\Phi }}_2-{\mathrm{\Phi }}_3=0.$$

(48)

Differentiating (47) with respect to time and considering (41) and the initial conditions will result in the following:

$${\displaystyle \frac{d{\mathrm{\Phi }}_1}{dt}}+{\displaystyle \frac{1}{{\mathrm{\rho }}_0}}\left({\displaystyle \frac{l_2}{S_2^2}}{\mathrm{\nu }}_2^{\partial }(B_2){\mathrm{\Phi }}_2+\left({\mathrm{\rho }}_0+{\mathrm{\rho }}_2(B_2)\right)\right){\displaystyle \frac{d{\mathrm{\Phi }}_2}{dt}}-{\displaystyle \frac{d{\mathrm{\Phi }}_3}{dt}}=0,.$$

(49)

Six algebraic differential equations are solved jointly. There are six unknowns: three winding currents and three magnetic fluxes in core parts.

$$w_1{\displaystyle \frac{d{i}_1}{dt}}-\left[{\displaystyle \frac{l_1}{S_1^2}}{\mathrm{\nu }}_1^{\partial }(B_1){\mathrm{\Phi }}_1+\left({\mathrm{\rho }}_{01}+{\mathrm{\rho }}_1(B_1)\right)\right]{\displaystyle \frac{d{\mathrm{\Phi }}_1}{dt}}+\left[{\displaystyle \frac{l_2}{S_2^2}}{\mathrm{\nu }}_2^{\partial }(B_2){\mathrm{\Phi }}_2-{\mathrm{\rho }}_2(B_2)\right]{\displaystyle \frac{d{\mathrm{\Phi }}_2}{dt}}=0,.$$

(50)

$$w_3^{\langle 1\rangle }{\displaystyle \frac{d{i}_3^{\langle 1\rangle }}{dt}}-w_3^{\langle 2\rangle }{\displaystyle \frac{d{i}_3^{\langle 2\rangle }}{dt}}-\left[{\displaystyle \frac{l_2}{S_2^2}}{\mathrm{\nu }}_2^{\partial }(B_2){\mathrm{\Phi }}_2+{\mathrm{\rho }}_2(B_2)\right]{\displaystyle \frac{d{\mathrm{\Phi }}_2}{dt}}-\left[{\displaystyle \frac{l_3}{S_3^2}}{\mathrm{\nu }}_3^{\partial }(B_3){\mathrm{\Phi }}_3+{\mathrm{\rho }}_3(B_3)\right]{\displaystyle \frac{d{\mathrm{\Phi }}_3}{dt}}=0,$$

(51)

where

$${\mathrm{\nu }}_j(B_j)={\displaystyle \frac{H(B_j)}{B_j}}=65+270{B_j}^4 , {\mathrm{\nu }}_j^{\partial }(B_j)={\displaystyle \frac{\partial H(B_j)}{\partial B_j}}=65+1350{B_j}^4;$$

(52)

$${\mathrm{\rho }}_j(B_j)={\displaystyle \frac{{\mathrm{\nu }}_j(B_j)l_j}{S_j}} , B_J={\displaystyle \frac{{\mathrm{\Phi }}_j}{S_j}}, j=1, 2, 3..$$

(53)

The system (44)–(46) and (49)–(51) will be expressed in the matrix-vectoral format:

$$A\cdot {\displaystyle \frac{dX}{dt}}=B\Rightarrow {\displaystyle \frac{dX}{dt}}=A^{-1}\cdot B, \det A\ne 0,$$

(54)

where A is the matrix of nonlinear coefficients, X is the columnar vector of unknown functions of currents and magnetic fluxes, and B is the columnar vector of the right equation sides.

The input nonlinear differential state equations of the electromagnetic device (54) will finally be presented in coordinate format:

$$\begin{gathered} {\displaystyle \frac{d}{dt}}\left[\begin{array}{c}{i}_1\\ i_3^{\langle 1\rangle }\\ i_3^{\langle 2\rangle }\\ {\mathrm{\Phi }}_1\\ {\mathrm{\Phi }}_2\\ {\mathrm{\Phi }}_3\end{array}\right]={\left[\begin{array}{cccccc}{L}_{\mathrm{\sigma }1}& & & w_1& & \\ & L_{\mathrm{\sigma }3}^{\langle 1\rangle }& & & & w_3^{\langle 1\rangle }\\ & & L_{\mathrm{\sigma }3}^{\langle 2\rangle }& & & w_3^{\langle 2\rangle }\\ w_1& & & -\left[{\displaystyle \frac{l_1}{S_1^2}}{\mathrm{\nu }}_1^{\partial }(B_1){\mathrm{\Phi }}_1+\left({\mathrm{\rho }}_{01}+{\mathrm{\rho }}_1(B_1)\right)\right]& {\displaystyle \frac{l_2}{S_2^2}}{\mathrm{\nu }}_2^{\partial }(B_2){\mathrm{\Phi }}_2-{\mathrm{\rho }}_2(B_2)& \\ & w_3^{\langle 1\rangle }& -w_3^{\langle 2\rangle }& & -\left[{\displaystyle \frac{l_2}{S_2^2}}{\mathrm{\nu }}_2^{\partial }(B_2){\mathrm{\Phi }}_2+{\mathrm{\rho }}_2(B_2)\right]& -\left[{\displaystyle \frac{l_3}{S_3^2}}{\mathrm{\nu }}_3^{\partial }(B_3){\mathrm{\Phi }}_3+{\mathrm{\rho }}_3(B_3)\right]\\ & & & 1& {\displaystyle \frac{1}{{\mathrm{\rho }}_0}}\left[{\displaystyle \frac{l_2}{S_2^2}}{\mathrm{\nu }}_2^{\partial }(B_2){\mathrm{\Phi }}_2+\left({\mathrm{\rho }}_0+{\mathrm{\rho }}_2(B_2)\right)\right]& -1\end{array}\right]}^{-1}\times \\ \left[\begin{array}{c}{u}_1-r_1{i}_1\\ u_3^{\langle 1\rangle }-r_3^{\langle 1\rangle }{i}_3^{\langle 1\rangle }\\ u_3^{\langle 2\rangle }-r_3^{\langle 2\rangle }{i}_3^{\langle 2\rangle }\\ 0\\ 0\\ 0\end{array}\right]. \end{gathered}$$

(55)

Similar systems are commonly created including expanded, over-dimensioned magnetic circuits, where the magnetic circuit characteristic can be assumed to be linear. In the event, the magnetic circuit’s electromagnetic state equations are considerably simplified:

$$w_1{i}_1=\left({\mathrm{\rho }}_0+{\mathrm{\rho }}_1\right){\mathrm{\Phi }}_{k1}+{\mathrm{\rho }}_2({\mathrm{\Phi }}_{k1}-{\mathrm{\Phi }}_{k2})=\left({\mathrm{\rho }}_0+{\mathrm{\rho }}_1\right){\mathrm{\Phi }}_1-{\mathrm{\rho }}_2{\mathrm{\Phi }}_2;.$$

(56)

$$w_3^{\langle 1\rangle }{i}_3^{\langle 1\rangle }-w_3^{\langle 2\rangle }{i}_3^{\langle 2\rangle }={\mathrm{\rho }}_3{\mathrm{\Phi }}_{k2}+{\mathrm{\rho }}_2({\mathrm{\Phi }}_{k1}-{\mathrm{\Phi }}_{k2})={\mathrm{\rho }}_2{\mathrm{\Phi }}_2+{\mathrm{\rho }}_3{\mathrm{\Phi }}_3.$$

(57)

$$w_3^{\langle 1\rangle }{\displaystyle \frac{d{i}_3^{\langle 1\rangle }}{dt}}-w_3^{\langle 2\rangle }{\displaystyle \frac{d{i}_3^{\langle 2\rangle }}{dt}}-{\mathrm{\rho }}_2{\displaystyle \frac{d{\mathrm{\Phi }}_2}{dt}}-{\mathrm{\rho }}_3{\displaystyle \frac{d{\mathrm{\Phi }}_3}{dt}}=0,$$

(58)

Assuming the linearity of the B vs. H relationship, the Equation (47) will be reformulated as below:

$${\displaystyle \frac{d{\mathrm{\Phi }}_1}{dt}}+\left(1+{\displaystyle \frac{{\mathrm{\rho }}_2}{{\mathrm{\rho }}_0}}\right){\displaystyle \frac{d{\mathrm{\Phi }}_2}{dt}}-{\displaystyle \frac{d{\mathrm{\Phi }}_3}{dt}}=0.$$

(59)

Hence, the input linear differential electromagnetic state equations of the device will be expressed in the normal form:

$${\displaystyle \frac{d}{dt}}\left[\begin{array}{c}{i}_1\\ i_3^{\langle 1\rangle }\\ i_3^{\langle 2\rangle }\\ {\mathrm{\Phi }}_1\\ {\mathrm{\Phi }}_2\\ {\mathrm{\Phi }}_3\end{array}\right]={\left[\begin{array}{cccccc}{L}_{\mathrm{\sigma }1}& & & w_1& & \\ & L_{\mathrm{\sigma }3}^{\langle 1\rangle }& & & & w_3^{\langle 1\rangle }\\ & & L_{\mathrm{\sigma }3}^{\langle 2\rangle }& & & w_3^{\langle 2\rangle }\\ w_1& & & -\left({\mathrm{\rho }}_{01}+{\mathrm{\rho }}_1\right)& {\mathrm{\rho }}_2& \\ & w_3^{\langle 1\rangle }& -w_3^{\langle 2\rangle }& & -{\mathrm{\rho }}_2& -{\mathrm{\rho }}_3\\ & & & 1& \left(1+{\displaystyle \frac{{\mathrm{\rho }}_2}{{\mathrm{\rho }}_0}}\right)& -1\end{array}\right]}^{-1}\left[\begin{array}{c}{u}_1-r_1{i}_1\\ u_3^{\langle 1\rangle }-r_3^{\langle 1\rangle }{i}_3^{\langle 1\rangle }\\ u_3^{\langle 2\rangle }-r_3^{\langle 2\rangle }{i}_3^{\langle 2\rangle }\\ 0\\ 0\\ 0\end{array}\right].$$

(60)

$${\mathrm{\rho }}_j={\displaystyle \frac{l_j}{{\mathrm{\mu }}_j{S}_j}} , j=0, 1, 2, 3 , {\mathrm{\mu }}_k={\mathrm{\mu }}_k^{*}{\mathrm{\mu }}_0, {\mathrm{\mu }}_k^{*}=12000-15000.$$

(61)

A system of nonlinear differential equations will be solved jointly: (55) with (52) and (53) or a system of linear differential Equation (60) with (61).

5. Computer Simulation Results

The nonlinear variant, i.e., (52), (53), and (55), has been employed to obtain the current and magnetic flux waveforms. The solution relies on the Runge–Kutta method with an adapted integration step and the procedure of reversing the matrix of coefficients (A) at every step of Gauss integration. MATHCAD ver. 15/Prime 2 (Win 7) is utilised to this end.

A static dependence of electric steel magnetising is shown in Figure 5 as a reverse steel magnetising curve $H(B)=\mathrm{\nu }(B)B=65B+270B^5$. Figure 5 shows that the steel curve becomes deeply saturated at magnetic inductions above 1.5 T. Since the magnetic branch solid is generally an isotropic environment, the tensor of the reverse magnetic permeability becomes a scalar.

Static relationships for the dynamic and static reverse magnetic permeability as a function of magnetic induction are illustrated in Figure 6 and Figure 7. Dynamic reverse permeability, in contrast to static reverse permeability, exhibits a rapid field intensity growth $H(B)$ at lower magnetic inductions B. This is reasonable in both mathematical (cf. (52)) and physical senses. The dynamic relationship shows an element’s properties in the transient state, according to the definition of total derivative with respect to time; see (41).

In Figure 8, Figure 9 and Figure 10, transient currents across three electromagnetic circuit windings for the first and third magnetic branches are shown, cf. Figure 3 and Figure 4. The currents across winding electric circuits are similar. The current amplitude is maximum in the lower coils of the third magnetic circuit branch (Figure 3).

Transient magnetic flux waveforms for each of the three respective magnetic circuit branches are shown in Figure 11, Figure 12 and Figure 13.

Transient responses are similar for the second and third branches, but clearly distinct for the first. The amplitude of flux oscillations across it is about 8-times lower than in the other two branches. The cause is well known. There is an air gap in the first branch, which generates very large reluctance. If magnetic fluxes are compensated, the algebraic total across all of the integration zone will be zero (the Kirchhoff’s first law for magnetic circuits).

Figure 14 and Figure 15 show transient dynamic and static waveforms of reverse magnetic permeabilities for the second magnetic branch of the electromagnetic circuit.

The oscillations of these functions are peculiar. The functions depend on the polynomial approximations to the dynamic and static waveforms; see (52). Both the functions are even.

6. Recapitulation

The adoption of the classical approach to mathematical modelling of electromagnetic circuits in order to derive the circuit state equations results in the decomposition of the uniform electromagnetic system into two subsystems. The first is the magnetic subsystem, which contains the magnetic circuit of the analysed object, and the second is the electrical subsystem, which contains the electrical circuits. For this purpose, so-called formalised methods are quite often used. For the first subsystem, the structural magnetic state equations are formulated, which describe the topology of the circuit (these equations are derived based on the circuit graph), and for the second subsystem, the equations of the electric state are formulated, which describe all the elements linked with the branches of the magnetic circuit. In the case of variational approaches, the situation is different, i.e., the method of decomposing the uniform electromagnetic circuit is not used to formulate the mathematical model of the electromagnetic circuit. Instead, four energy functions are defined, which are terms of the non-conservative Lagrangian (see Equation (3)). As a result of equating the variation of the above-mentioned functional to zero, the well-known Euler–Lagrange’s equation is obtained, the solution of which allows for the formulation of equations describing the dynamic system, in our case the electromagnetic circuit. When considering the choice between variational and classical approaches, it is difficult to indicate the advantage of one of them. It is obvious that when both approaches are correctly applied, the equations of state of the analysed object must be identical. A certain advantage of variational methods over classical approaches is the possibility of taking into account some latent motions in a single integrated system composed of a significant number of subsystems assigned to different fields of science. Ultimately, the choice of the method is up to the potential researcher, depending on his knowledge.

7. Conclusions

The modification of the integral variational Hamilton–Ostrogradsky’s principle significantly expands the possibilities of interdisciplinary variational approaches and gives them the opportunity to compete with classical methods.

In order to formulate a mathematical model of an electromagnetic circuit based on variational approaches, it is necessary to correctly represent the extended non-conservative Lagrange’s function using analytical forms of all four energy terms of this function.

In the case of high-order electromagnetic circuit analysis, formalised calculation methods should be used to formulate structural equations of the magnetic part of the integrated electromagnetic circuit.

For mathematical modelling of transient processes in nonlinear electromagnetic systems, nonlinear functions of the inverse of magnetic permeabilities should be used. On the one hand, this simplifies the mathematical form of the system of differential dynamic state equations of the object, and on the other hand, it reduces their stiffness.