Real-Time Digital Twin of a Wound Rotor Induction Machine Based on Finite Element Method

: Monitoring and early fault prediction of large electrical machines is important to maintain a sustainable and safe power system. With the ever-increasing computational power of modern processors, real-time simulation based monitoring of electrical machines is becoming a topic of interest. This work describes the development of a real-time digital twin (RTDT) of a wound rotor induction machine (WRIM) using a precomputed ﬁnite element model fed with online measurements. It computes accurate outputs in real-time of electromagnetic quantities otherwise difﬁcult to measure such as local magnetic ﬂux, current in bars and torque. In addition, it considers space harmonics, magnetic imbalance and fault conditions. The development process of the RTDT is described thoroughly and outputs are compared in real-time to measurements taken from the actual machine in rotation. Results show that they are accurate with harmonic content respected.


Introduction
Real-time simulation based monitoring of physical systems through digital twins is getting more attention due to the ever-increasing computational power of modern processors [1,2]. The term "digital twin" holds many definitions, each of them being more or less different, but it is commonly known as a set of accurate models capable of representing a physical system throughout all of its life cycle, from design phase to operation [3,4]. Real-time monitoring with a digital twin fed with sensor data coming from the actual system has many benefits, such as fault detection, output optimization and downtime planning. The challenge lies in developing accurate multi-physical models able to interact between each others.
For electrical machines, the core model is the electromagnetic model, which computes electrical quantities such as currents, voltages and electromagnetic torque. As of now, dq models are used extensively for real-time simulation due to their relative simplicity and small computation time, mainly for hardware-in-the-loop testing purpose [5] and grid simulation [6]. A computationally efficient model is needed because real-time constraints impose that the simulation time step be larger than the computation time it takes in real-world clock [7]. However, dq models are based on numerous assumptions; phenomenons like space harmonics or magnetic imbalances are neglected. Furthermore, they only provide information about the quantities seen at the machine terminals.
That said, the only way to obtain a high order model capable of representing accurately any machine structure and computing local quantities as magnetic densities and current densities is by using the well-known finite element method (FEM). It computes locally the Maxwell's equation solution for any given machine geometry considering even local magnetic saturation. The downside of this technique is the enormous computation time it takes for magnetodynamic resolution which make it challenging for real-time implementation.
Several detailed models of electrical machines have been developed as an alternative to FEM in an effort to make compromises between accuracy and computational requirement for implementation on real-time simulators [8]. Among these, magnetic equivalent circuits (MEC) models and lumped circuit models using winding functions are the most prominent.
MEC [9] is a model based on a network of flux tubes capable of representing the unique geometry and winding distribution of a machine. In [10], an induction motor was implemented for real-time execution on a high clock speed FPGA with a time step of 500 µs. A time step of 150 µs is attained in [11] for a switched reluctance motor. While adaquate accuracy can be achieved at a fraction of the time required by FEM, the amount of calculations required in one time step remains high.
On another hand, lumped circuit models using winding function calculations [12] were also implemented for real-time. In [13], a synchronous machine model was implemented at 20 µs time step using a CPU-based real-time digital simulator. These models considers the geometry of the machine to a certain extent only, because simplifications are made to make analytical calculations manageable and efficient.
However, it was demonstrated in [14] that, by using a magnetically coupled circuit approach using precomputed inductance functions with a FEM software, it is possible to reach relatively small computation times while keeping the accuracy of FEM. This model is called CFE-CC-Combination of Finite Element with Coupled Circuits [15]. It has three major strong points: 1. It considers space harmonics and magnetic imbalances since any machine geometry can be modeled in a FEM software; 2. Any number of electrical circuits can be added to the models, including phase windings of course, but also bars, dampers and even search coils. This means the model can output actual currents flowing in bars and dampers and also local magnetic fluxes by using search coils. 3. It can compute distribution of power losses inside the machine [16], allowing calculations of local temperature elevation if used with a thermal model [17].
For these reasons, the CFE-CC model is very close to its physical counterpart. All it takes for building it is the FEM model, which is available from the design phase.
This work presents the methodology to implement for real-time execution the CFE-CC model fed with measurements taken from the real machine, thus allowing effective online monitoring of internal electromagnetic quantities. The proposed electromagnetic model can then be used alongside thermal and mechanical models to complement this real-time digital twin (RTDT). The machine used as case-study to realize the proposed RTDT is a wound rotor induction machine (WRIM), also known as doubly-fed induction machine. This particular type of machine was chosen because all of its electrical circuits are accessible for measurements, thus making it convenient for signal waveform validation. It is important to note however that the described method herein is applicable to every type of electrical machine. After a detailed explanation on the model's construction and implementation on digital real-time simulator (DRTS), results are validated by comparing in real-time the RTDT's outputs to the physical machine's measurements.

Model's Equations
The model used for the purpose of this work does not include magnetic saturation. Iron losses as well as capacitive effects are also neglected.
A wound rotor induction machine, as well as many other types of electrical machines, comprises n electrical circuits where current can flow, such as the phases and bars. The total magnetic flux φ at a given circuit is the sum of the magnetic fluxes received from all other circuits and itself. It is worth noting that every quantity referenced in this manuscript is instantaneous and no referential transformation is used.
The magnetic flux received at a given circuit from another one depends on the current i flowing in the sender and the magnetic coupling L between the two. The coupling is a function of the rotor angular position θ.
All aforementioned equations are included into a single matrix equation: Flux variation seen by any circuit creates a voltage v given by: By combining Equations (3) and (4), the dynamic system of Equation (5) is obtained.
As for electromagnetic torque T em , it can be computed using: For search coils or any electrical circuit always left open, Equation (7) is used instead to compute back-emf v w for better numerical stability of ordinary differential equations solvers. Matrix [L w ] defines the coupling between the search coils and the other circuits carrying a current.
To sum up, identification of the model requires to know the inductance matrix function [L(θ)] which can then be inverted and differentiated. What makes the CFE-CC model versatile is the identification process of [L(θ)] with an FEM software. It uses a lookup table to store the inverse and derivative inductance matrix for many discrete positions of the rotor.
As in [14], it is possible to consider magnetic saturation with the precomputed inductance functions by using a correction coefficient afterward on computed magnetic flux linkages. Performances and computing time differ depending on type of machine and chosen quantities to compute the coefficient. Table 1 gathers specifications of the WRIM used in this work. As depicted in Figure 1b, there are three windings on the stator (a s , b s , c s ) and three windings on the rotor (a r , b r , c r ). A small search coil is also added around one stator tooth (w s ) in order to measure the voltage across it during operation. These seven circuits make up the CFE-CC model.

Identification of Electrical Circuits
The resistance matrix [R] of the WRIM is diagonal as shown below. Each terminal is accessible so resistance values can easily be measured. The search coil is not included because it is left as an open circuit, thus Equation (7) is used instead of (5) for this particular case.
As for [L], for any given rotor position, it has the form shown below. Magnetic couplings of the search coil are also included in a separate matrix [L w ].
L asas L asbs L ascs L asar L asbr L ascr L bsas L bsbs L bscs L bsar L bsbr L bscr L csas L csbs L cscs L csar L csbr L cscr L aras L arbs L arcs L arar L arbr L arcr L bras L brbs L brcs L brar L brbr L brcr L cras L crbs L crcs L crar L crbr L crcr [L w ] = L wsas L wsbs L wscs L wsar L wsbr L wscr (10)

Computation of the Inductance Matrix Using FEM
Computation of [L(θ)] is performed with an FEM software and accurate plans of the machine geometry and windings. A 3D FEM software is an ideal choice as it consider geometries that are impossible to model in a 2D environment. However, a 2D FEM software is normally easier to use and less computationally expensive. It should be noted that computation time required to solve the finite element domain is not related to the CFE-CC model's computation time at each time step, because results of the FEM are used as precomputed lookup tables.
For the present work, the authors opted for a 2D software to quickly compute [L(θ)]. The geometry of the WRIM makes it convenient to use a 2D space, but a few issues cannot be accounted in 2D alone such as rotor skewing and coil ends. Their respective effects are added to the model by altering the computed inductance matrix. Details of this procedure are presented at Appendix A. Figure 2 shows the FEM model of the WRIM in FLUX2D, the FEM software used for this work.
The domain can be reduced to only half of the complete machine because the machine has two pairs of poles and presents a symmetry. All magnetic materials are assumed to be linear, so saturation is not considered. Meshing is performed using the software's automatic mesh tool. Mesh size was selected by comparing the solution obtained using a fine mesh to more and more coarse ones, until the error starts increasing noticeably. Once the FEM model is ready, the procedure for identifying magnetic couplings is launched. Inductance matrix [L] is computed for many rotor positions to build the lookup table. Inductance matrix [L(θ)] is obtained by first computing the inductance matrix of the machine without windings, i.e., the inductance matrix of the slots [L slot (θ)]. That way, one may generate any winding configuration without restarting the FEM process simply by applying the following matrix transformation: (11) [N] is the winding configuration matrix. It contains the number of turns N of each electrical circuit in each slot. Equation (12) shows how to build the matrix, where n c is the total number of circuits and n s is the total number of slot.
In total, 1440 positions were computed evenly for a 180 • rotation, i.e., 60 positions per stator slot pitch or an angular step of 0.125 • . It is plenty to have a good resolution. Impact of inductance discretization is described thoroughly in [18]. At this point, it is important to take into consideration the physical memory to store the lookup table. Since the WRIM has six circuits without the search coil, the total number of elements in the lookup table for the inverse inductance matrix is 51,840. Another lookup table of the same size is necessary for the derivative if torque computation is included.

Real-Time Simulator And RT-Lab
The processing unit of the digital twin is an OP4510 entry-level DRTS of OPAL-RT. Real-time computations are performed by a dedicated multi-core CPU capable of time step of 2 µs. For applications requiring smaller time step, one may rather opt for a higher-end DRTS or FPGA based platforms instead. Table 2 sums up the main technical specifications of the DRTS used in this work. The OP4510 needs to be connected via an ethernet cable to a personal computer (PC) with RT-LAB installed. RT-LAB is a software environment for real-time simulation to interface with OPAL-RT's DRTS. Using a PC and RT-LAB, user can: 1. Modify execution state of the program (run, stop, pause); 2. Change parameters in real-time of the running model; 3. Receive data from the running model.

Voltage and Current Measurement
Measured voltages are the inputs of the model. Currents are also measured, but they are used for validation purposes to compare with the RTDT's outputs. A printed circuit board (PCB) was designed for voltage and current acquisition from the machine. The signals are then sampled using the DRTS's analog input channels.

Angular Position Sensor
The angular position of the rotor is needed to retrieve the inductance matrix from the lookup table. It is measured using an absolute shaft encoder. The measured position is encoded in a 12 bits signal sent through a parallel communication interface. Each bit is acquired by a digital input channel of the DRTS and binary to decimal conversion is performed in software. Figure 3 shows an overview of the hardware setup. The CFE-CC model is executed in real-time along with the actual machine in operation. All estimated waveforms can be displayed on an oscilloscope using the DRTS's analog output channels.

Programming Method
Programming of the RTDT was realized entirely within the MATLAB/Simulink environment using RT-LAB. It generated C code from a Simulink block diagram and sent it to the DRTS for compilation and execution.

Position Tracking Scheme
The acquired position signal cannot be fed directly to the CFE-CC model after decimal conversion because it is a noisy discrete signal. The model would see speed impulsions that would impact its estimations. In order to filter the position signal, the control loop shown at Figure 4 was programmed in the DRTS. The resulting closed loop transfer function is: where K p and K i are the PI controller's proportional and integral gains. The measured position was transformed beforehand from a sawtooth waveform into a ramp. The difference between measured and estimated position was sent to a PI controller and an integrator. Hence, the error in steady state was theoretically reduced to zero. The controller's gains should be chosen high enough to keep track of the speed variations of the machine.

CFE-CC Model Implementation for Real-Time Execution
Resolution of the global system was achieved using the State-Space-Nodal (SSN) solver of OPAL-RT available in the ARTEMIS package [19]. It allowed us to solve power systems equations with a hybrid method between pure state-space and nodal approach [20]. It was built upon the SimPowerSystems (SPS) package, thus the Simulink diagram was designed the same way by using SPS blocks (resistances, sources, inductances, etc.).
A custom block made for SSN resolution was implemented for the CFE-CC model. Figure 5 shows the WRIM bloc containing the CFE-CC model's equations as well as external components connected to it. The stator's terminals were connected to controlled voltage sources, which were the measured stator voltages of the physical machine. The rotor's terminals were connected to external resistors adjustable in real-time from the PC. By placing the machine in a separate nodal group than the external components, it had a decoupling effect with SSN resolution, thus greater stability was achieved even with large external resistances [21]. In order to implement such a model for SSN resolution, one must first select which circuits are available to the user as external connections, i.e., the Simulink block's terminals of the machine. The WRIM had six circuits, all available for external connection, thus vectors [v in ] and [i out ] were each six units large. They were respectively the input voltages and output currents at the Simulink block's terminals. The search coil w s was not treated as a circuit because no current could flow through.
[φ] = φ as φ bs φ cs φ ar φ br φ cr T (14) [v in ] = v as v bs v cs v ar v br v cr T (15) [i out ] = i as i bs i cs i ar i br i cr T (16) Dynamical system of Equation (5) is rearranged as standard state-space form. To facilitate reading of the following equations, the θ dependency is omitted. (19) [B] = [I] (20) [ and [L(θ)] −1 were stored in a lookup table. At each time step, the exact matrices were computed as the result of an interpolation between two adjacent matrices of the lookup table by using the angular position measurement. The same type of storage and interpolation was used to retrieve the derivative inductance matrix for torque computation using Equation (6). Now regarding the search coil, since no current, i.e., no state, was associated with it, it was not included in the dynamical system. Instead, it was computed using Equation (7).

Time Step Selection
The execution time step of the program was chosen to be as small as possible without any risk of overrun. An overrun occurred if the DRTS did not finish all computations within the specified time step. In that case, a time step was skipped resulting in erroneous estimation. RT-LAB gives information regarding computational requirement of the program. Based on it, a time step of 6 µs was chosen. It was important to have a small time step particularly when a circuit of the CFE-CC model was connected to a large external resistance. A large external resistance created a stiff state-space system, thus requiring a smaller time step to keep the solver stable. The SSN solver increased stability compared to the traditional approach but accuracy decreases when the time step was too large.

Validation
Validation of the RTDT was performed by comparing its real-time outputs to measurements available on the physical WRIM during operation. Waveforms shown below are real-time raw data acquired by an oscilloscope. For comparison purpose, a dq model of the WRIM was also implemented on another core of the DRTS for parallel execution alongside the CFE-CC model.
For the first setup, shown in Figure 6a, the WRIM's stator windings were connected to a three-phase 60 Hz autotransformer. It was used as an induction motor, so rotor terminals were short-circuited. The WRIM was mechanically coupled to a synchronous generator loaded by resistors. Figures 7 and 8 show waveform comparisons between measurements, CFE-CC model and dq model. Only winding a is shown because b and c were only 120 • phase-shifted. Both models yielded good accuracy, but it was possible to draw more conclusions from the frequency decomposition. Tables 3-5 show the notables frequencies of the spectra. As expected, the CFE-CC model approximated well higher frequencies which were related to the geometry of the machine, whereas the dq model only computed harmonics related to input voltage.

Frequency (Hz) Measurements (V)
60 156.5 180 1.0 300 3.9 420 1.8 Table 4. Main frequency components in stator current.   Geometric irregularities of the machine had an impact on the torque as well, as seen from Figure 9. By analysing the electromagnetic torque of a machine, one can draw many conclusions about its condition.    Figure 11 shows the dynamic response of rotor currents when the mechanical load was suddenly disconnected from the shaft. Here, the position tracking control loop played an important role. It needed to be fast enough to minimize the error during speed variations. Figure 11. Rotor current in winding a r starting from load operation and suddenly disconnecting the load.

Frequency (Hz) Measurements (A) CFE-CC (A) dq (A)
The following Figures 12-14 show the RTDT's current outputs when a phase imbalance was introduced. The setup is depicted at Figure 6b. A resistor was connected to winding a r of the physical machine, and the same was done to the RTDT. While the frequency content was accurate, a small phase discrepancy was present at winding b r and c r . The cause could be related to the methodology employed to compute [L] in Appendix A, such as the way the coil end inductances were added. Many assumptions were made. Magnetic saturation could also play a role since it was not taken into account in the model.

Conclusions
The detailed CFE-CC model of a WRIM was implemented in real-time on an entry-level DRTS of OPAL-RT with a time step of 6 µs. This model considers the machine's particular geometry and winding configuration can be changed easily using a simple matrix transformation. Any number of search coils can be added to the model, as well as faults such as broken bar and eccentricity. Comparisons with real-time measurements have shown that the CFE-CC model is accurate and computationally efficient. This work paves the way to new real-time monitoring techniques for electrical machines using the CFE-CC model alongside others in a digital twin. However, the machine studied in this work has only seven circuits. Including more electrical circuits in the CFE-CC model results in more computations at each time step and more memory requirement for lookup tables. Future works include studying the feasibility of real-time execution of the CFE-CC model for large synchronous machines comprising many electrical circuits. and reduce torque oscillations. The skew angle θ sk is generally given by the manufacturer. Fortunately, it is possible to add this effect into the inductance matrix [L(θ)] previously identified in 2D [22]. By dividing an unskewed machine along its length into M several slices, or sub-machines, the inductance matrix [L m (θ)] of each sub-machine is simply the initial one [L(θ)] divided by M. As for the WRIM of this work, θ sk is equal to 7.5 • , which corresponds to a slot pitch. The machine was divided into 61 sub-machines. Figure A1 shows the impact of the added skew effect on rotor currents during load operation, with stator windings fed with perfectly sinusoidal voltages and rotor windings short-circuited. Current ripples are reduced as expected.

Appendix A.2. Coil Ends
Coil ends can be modeled in a 3D FEM software, whereas it cannot be in 2D. In the latter case such as in the present work, their inductive effect must be added to [L(θ)] by other means. For simplicity, we ignore the impact on mutual inductances. The technique used is to conduct a standstill frequency response (SSFR) test [14] and add a self-inductance to the windings in order to make the experimental response fit the simulated one. The test was conducted while rotor windings are short-circuited. Figure A2 compares the response of the FEM model to the experiment with and without coil ends. Inductances of 0.01 H and 0.00047 H were added to stator and rotor windings respectively, which is less than 10% of the magnetizing inductance. Figure A2.
Inductance seen from the stator versus frequency with and without coil ends, rotor terminals short-circuited.