1 Introduction↩︎

Transformers and electric machines rely on laminated ferromagnetic cores to guide and concentrate magnetic flux. Accurate modeling of electromagnetic fields and magnetic losses within these laminated regions is crucial for the reliable prediction of induced voltages and overall device efficiency [1], [2]. However, the electromagnetic behavior of laminated cores is particularly challenging to model due to the inhomogeneous distribution of eddy currents, the nonlinear and history-dependent nature of magnetic hysteresis, and the coupling between these two phenomena.

Modeling hysteresis alone is already challenging. The Preisach [3] and Jiles-Atherton [4] models have long been considered standard hysteresis models, and have been widely implemented in two-dimensional (2D) finite element (FE) simulations [5]–[7]. However, both models were originally formulated for scalar hysteresis, and their extensions to vector hysteresis lack a clear physical foundation which limits their applicability to vector field problems. In addition, they may suffer from limitations related to accuracy, computational cost, or numerical convergence [8], [9]. More recently, the energy-based (EB) model [10], [11] has emerged as a promising alternative. This hysteresis model is based on energy considerations and is intrinsically vectorial [12]. The forward EB model takes the magnetic field \(\mathbf{H}\) as input and returns the history-dependent magnetic flux density \(\mathbf{B}\), while its inverse formulation [13], [14], computes \(\mathbf{H}\) from a given \(\mathbf{B}\). The forward model can be evaluated either by solving a minimization problem [15] or a nonlinear vector function [11]. Besides these two formulations, a simplified explicit variant, akin to a vector play model [10], [11], also exists and provides a computationally efficient update at the cost of a slightly reduced accuracy [16]. The inverse formulation, by contrast, requires solving a minimization problem [11], [15], [17], whose robustness has been recently significantly improved [13], [14].

Eddy currents further add to the complexity of modeling ferromagnetic behavior. They form closed current loops whose distribution is governed by the spatial and temporal variations of the magnetic flux density. In laminated cores, eddy currents are largely confined to the plane of each lamination, but three-dimensional (3D) effects arise near edges and geometric discontinuities and cannot be fully captured by the 2D formulations commonly employed in electrical machine design. However, direct 3D modeling of laminated cores, in which the fields are resolved within each individual lamination, is several orders of magnitude more computationally expensive. Although feasible on high-performance computing platforms, such simulations remain impractical for engineering design and optimization. The periodicity of the laminated structure can, however, be exploited through homogenization techniques, yielding 2D formulations that accurately capture the bulk eddy-current behavior, although their accuracy deteriorates near edges and geometric discontinuities. In multiscale methods like the heterogeneous multiscale method (HMM) [18], [19] or the multilevel finite element method (FE\(^2\)) [20], the effective material response unavailable at the macroscopic scale is computed during the simulation by solving a mesoscopic-scale problem on a representative volume element and homogenizing its solution. In the present context, the mesoscopic-scale model can be formulated as a magneto-quasistatic FE simulation solved on the cross-section of a lamination and incorporating the EB hysteresis model. However, this approach remains too computationally demanding for routine machine design for three main reasons: i) The mesoscopic-scale model must be queried in each macroscopic finite element, at each time step and at each nonlinear iteration. ii) The mesoscopic-scale model must additionally provide the corresponding Jacobian, \(\partial\mathbf{H}/\partial\mathbf{B}\) or \(\partial\mathbf{B}/\partial\mathbf{H}\), to construct the macroscopic tangent stiffness matrix required for a Newton-Raphson scheme. If obtained by finite differences, the mesoscopic-scale problem must be solved for the nominal state, and for additional perturbed states, further increasing the number of mesoscopic-scale model queries. iii) Each mesoscopic-scale query itself requires evaluating the EB model at every element of the mesoscopic mesh.

Advances in machine learning have led to a rapid expansion of neural network (NN)-based methods dedicated to magnetic modeling. Among other applications, NNs are used to predict core losses [21], to predict hysteresis model parameters [22], [23], or to predict \(\mathbf{B}\) from \(\mathbf{H}\) and inversely [24]–[27]. This enthusiasm is further evidenced by the emergence of dedicated international competitions in the power electronics community, such as the MagNet challenges, which have been specifically designed to stimulate innovation in the data-driven modeling of core losses [28] and transient magnetic behavior [29]. Yet, despite an intense research activity, and proven effectiveness in related fields such as heterogeneous elasto-plastic materials [30], integrating an NN-based surrogate model to account for vector hysteresis and eddy currents in 2D FE simulations of realistic electrical machines remains an open problem.

In a previous contribution [31], the relation required for magnetic vector-potential formulations was approximated by a parametric homogenized constitutive law \(\tilde{\mathbf{H}}(\mathbf{B},\dot{\mathbf{B}},p_k)\). The parameters \(p_k\) of this law were identified element-wise through a three-step procedure: i) A 2D macroscopic-scale anhysteretic FE simulation was performed to compute characteristic periodic \(\mathbf{H}(t)\) sequences. ii) Corresponding \(\mathbf{B}(t)\) responses were obtained by feeding the resulting \(\mathbf{H}(t)\) sequences to the aforementioned mesoscopic-scale model using the explicit variant of the forward EB model, and homogenizing its solution. iii) The parameters \(p_k\) were then identified by fitting the homogenized constitutive law \(\tilde{\mathbf{H}}(\mathbf{B}, \dot{\mathbf{B}}, p_k)\) to the computed \((\mathbf{B}(t), \mathbf{H}(t))\) pairs. This identification was carried out through NN training, which minimizes the fitting error over all sequences simultaneously and yields a smooth spatial distribution of the parameters \(p_k\) across the macroscopic domain. The resulting macroscopic 2D FE simulations efficiently captured the combined effects of hysteresis and eddy currents through the parametric constitutive law \(\tilde{\mathbf{H}}(\mathbf{B}, \dot{\mathbf{B}}, p_k)\). However, this approach is restricted to periodic steady-state operating conditions. Furthermore, the fixed analytical structure of the parametric law inherently limits the achievable accuracy, and the resulting macroscopic FE simulations may converge slowly.

In the present work, rather than assuming a fixed parametric structure, an effective constitutive relation \(\tilde{\mathbf{H}}(\mathbf{B})\) is directly learned from precomputed data generated by the aforementioned mesoscopic-scale model. A recurrent neural network (RNN) is used to naturally handle the history-dependent nature of hysteresis, lifting the periodicity restriction of the previous approach, while the Jacobian \(\partial\tilde{\mathbf{H}}/\partial\mathbf{B}\) is obtained exactly through automatic differentiation [32]. The improved approach is evaluated by modeling two real-world applications. It demonstrates excellent accuracy and robust convergence properties, while maintaining computational times comparable to those of conventional anhysteretic simulations. In addition, the proposed NN-based model provides an estimation of its error, which constitutes valuable information for assessing the reliability of the computed solution.

The paper is organized as follows. Section 2 first presents the mesoscopic-scale model used to generate the training datasets and to assess the NN-based model accuracy. Section 3 then details the model architecture, training dataset and training procedure, while Section 4 presents the conducted 2D FE simulations. Finally, Section 5 provides conclusions and perspectives for future work.

2 Lamination Model↩︎

This section presents the detailed lamination model used as reference throughout this paper. It provides a homogenized description of the fields at the mesoscopic scale, i.e., at the scale of the lamination thickness. The model combines the explicit variant of the EB hysteresis model [10] with a one-dimensional (1D) magnetodynamic FE formulation accounting for eddy currents. The EB model is first introduced, followed by the associated 1D FE simulation.

2.1 Quasi-Static Mesoscopic Model: The Energy-Based Hysteresis Model↩︎

The EB hysteresis model is grounded in thermodynamic principles; its theoretical derivation is presented in detail in [11], [15]. Without expanding on this derivation, the model captures the distributed nature of pinning strengths through an ensemble of \(N\) elementary cells. Each cell \(n\) decomposes the mesoscopic magnetic field \(\mathbf{h}\) into a reversible component \(\mathbf{h}_\text{rev}^n\) and an irreversible component \(\mathbf{h}_\text{irr}^n\), whose magnitude is constrained within a threshold \(\kappa^n\). As long as \(|\mathbf{h}_\text{irr}^n|\le \kappa^n\), the cell behaves reversibly. Once the threshold is exceeded, the cell transitions into irreversible behavior. Among the update laws proposed to describe this transition, we adopt the explicit variant, illustrated on Fig. 1 and expressed as \[\mathbf{h}_\text{rev}^n = \begin{cases} \mathbf{h}_{\text{rev},\text{p}}^n & \text{if } |\mathbf{h}-\mathbf{h}_{\text{rev},\text{p}}^n| \le \kappa^n\text{,}\\ \mathbf{h} - \kappa^n \frac{\mathbf{h}-\mathbf{h}_{\text{rev},\text{p}}^n}{|\mathbf{h}-\mathbf{h}_{\text{rev},\text{p}}^n|} & \text{otherwise,} \end{cases} \label{eq:update}\tag{1}\] where \(\mathbf{h}_{\text{rev},\text{p}}^n\) denotes the previous reversible component \(\mathbf{h}_\text{rev}^n\). From the weighted sum of the reversible fields \(\mathbf{h}_\text{rev} = \sum_{n=1}^{N} w^n \mathbf{h}_\text{rev}^n\), with the \(w^n\) parameters summing up to 1, the mesoscopic magnetic flux density is expressed as \[\mathbf{b} = \mu_0 \mathbf{h} + \left[ \mathscr{L}\left(\frac{|\mathbf{h}_\text{rev}|}{h_a}\right)\frac{J_a}{|\mathbf{h}_\text{rev}|} + \mathscr{L}\left(\frac{|\mathbf{h}_\text{rev}|}{h_b}\right)\frac{J_b}{|\mathbf{h}_\text{rev}|} \right] \mathbf{h}_\text{rev}, \label{eq:b}\tag{2}\] where \(\mu_0\) is the magnetic permeability of vacuum, \(\mathscr{L}\,(x)=\coth\,(x)-1/x\) denotes the Langevin function, and \(h_a\), \(h_b\), \(J_a\) and \(J_b\) are parameters of the EB model. These parameters, along with the \(\kappa^n\) and \(w^n\) parameters, are determined from measurements, and have to be updated whenever the material changes. This work employs \(N=11\) elementary cells for the modeling of M235-35A non-oriented electrical steel sheets, with the identified parameter values reported in 6.

2.2 Dynamic Mesoscopic Model: Laminations with Eddy Current and Hysteresis↩︎

The aforementioned EB model solely accounts for quasi-static hysteresis. To also account for dynamic effects, it is included into the FE simulation of a stack of ferromagnetic laminations. Considering that all the laminations see the same imposed macroscopic magnetic field \(\mathbf{H}\), and neglecting edge effects, simulating a single lamination is sufficient to describe the whole stack. Furthermore, the mesoscopic magnetic field \(\mathbf{h}\) and the mesoscopic magnetic flux density \(\mathbf{b}\) inside a lamination can reasonably be assumed parallel to the lamination plane and to vary exclusively within the lamination depth. Eventually, the symmetry with respect to the median plane of a lamination is exploited to simulate solely half of a lamination thickness \(d\) (see Fig. 2).

Let \(t\) denote time, \(\rho\) the electrical resistivity, and \((x,y,z)\) a Cartesian coordinate system, with \(x\) and \(y\) lying in the plane of the laminations and \(z\) in the stacking direction. Furthermore, let \(z=0\) correspond to the lamination surface and \(z=d/2\) to its center. The 1D magneto-quasistatic formulation then reads: Find the mesoscopic magnetic field \(\mathbf{h}(z,t)=(h_x(z,t), h_y(z,t))\) such that \[\begin{cases} \partial_t\Big( b_x\left(\mathbf{h}\right) , h_x' \Big)_{\left[0,d/2\right]} + \Big( \rho\,\partial_zh_x , \partial_zh_x' \Big)_{\left[0,d/2\right]}&=0\\ \partial_t\Big( b_y\left(\mathbf{h}\right) , h_y' \Big)_{\left[0,d/2\right]} + \Big(\rho\,\partial_zh_y , \partial_zh_y' \Big)_{\left[0,d/2\right]}&=0 \end{cases} \label{eq:lam}\tag{3}\] holds for all test functions \(h_x'\), \(h_y'\), where the notation \((\cdot,\cdot)_D\) denotes the inner product of its arguments over bounded domain \(D\). The macroscopic excitation \(\mathbf{H}(t)\) is imposed at the lamination surface, i.e. \(\mathbf{h}(0,t) = \mathbf{H}(t)\), and the zero-current condition at the mid-plane leads to the Neumann boundary condition \(\partial_z\mathbf{h}(d/2,t)=0\). The macroscopic magnetic flux density \(\mathbf{B}(t)\) is retrieved by averaging the mesoscopic magnetic flux density \(\mathbf{b}(\mathbf{h})=(b_x(\mathbf{h}),b_y(\mathbf{h}))\) over the half lamination thickness. This formulation accounts for the eddy currents and the associated skin effect in the thickness direction via the resistivity term \(\rho\,\partial_z\mathbf{h}\), and accounts for hysteresis via the \(\mathbf{b}(\mathbf{h})\) term, solution of the EB model.

In practice, the model is implemented as a self-contained C language code, solving the linear system, whose matrix is sparse, using an efficient LU sparse algorithm. The code is applied to a stack of M235-35A non-oriented electrical steel laminations, whose domain of analysis \(z \in \left[0,d/2\right]\) is discretized with 51 equidistant nodes. Given that the maximum differential permeability is \(\mu_\text{max}=\lim_{h\rightarrow0^+}(\mu_0+\mathscr{L}(h/h_a)\,J_a/h+\mathscr{L}(h/h_b)\,J_b/h)\approx25.5\) mH/m, obtained from (2 ), the resulting mesh size \(d/100=3.5~\mu\)m ensures that the skin depth \(\delta = \sqrt{\rho/(\pi f \mu_\mathrm{max})}\) is resolved by at least five nodes up to approximately 44 kHz. Thus, skin-effect phenomena are accurately captured over the frequency range considered in this work (see Section 3.2). The temporal discretization is chosen sufficiently fine to avoid nonlinear iterations, and the time derivatives are discretized with the implicit Euler scheme. Finally, the cells of the EB model are initialized assuming an initially demagnetized material subjected to a slowly increasing excitation field until reaching \(\mathbf{H}(t=0)\), such that no dynamic effects are initially present and \(\mathbf{B}(\mathbf{H}(t=0))\) lies on the virgin magnetization curve.

3 Architecture and Training of the NN-Based Surrogate Model↩︎

The choice of the NN architecture is of predominant importance. Since hysteresis inherently depends on the material history, the NN model must operate on time-sequence data. Some authors process the entire sequence at once using, for example, neural operators [33]. In FE simulations, however, the input sequence length evolves continuously during time integration, which makes whole-sequence architectures poorly suited for direct integration into incremental solvers. Another approach consists in using a fixed-length sliding window that captures only the most recent history [34]–[36]. Nevertheless, such a fixed-length window may omit parts of the history that significantly influence the response, thereby reducing accuracy. Alternatively, attention-based architectures [37] have shown high accuracy on sequence-to-sequence tasks [38], but the computational overhead becomes significant for long sequences. In this contribution, an RNN architecture is preferred. Indeed, RNNs are well suited to handle hysteresis in FE simulations: (i) they maintain and update state variables similarly as the EB model maintains and updates the reversible fields \(\mathbf{h}_\text{rev}^n\); (ii) they are fast to infer a single time step, making them efficient within incremental solvers; and (iii) they achieve high accuracy [39]–[41].

This section discusses the implementation of the NN-based model that will be integrated into the FE simulations of Section 4, and details its training process. At first, Section 3.1 presents the general model architecture and Section 3.2 details the datasets used for training. Section 3.3 later demonstrates that the model cannot be trained naively to minimize the mean squared error, since saturated fields will dominate the training to the detriment of lower fields where the magnetic behavior is more intricate. Section 3.3 therefore introduces a scaled loss, which is further modified in Section 3.4 to assess the model confidence. In Section 3.5, the model is upgraded to ensure correct prediction even in extreme saturation, which lies out of the training range. Section 3.6 later describes how the training sequences, initially defined on constant time steps, are resampled to include variable time steps, thereby enabling adaptive time-stepping strategies in the FE simulations. Finally, Section 3.7 discusses the hyperparameter tuning, while Section 3.8 assesses the model performance, prior to any FE coupling.

3.1 Baseline Architecture↩︎

To process temporal sequences, RNNs ingest a single time step at a time. At each time step \(k\), the network updates a hidden state \(S^k\) based on the current input \(X^k\) and the previous hidden state \(S^{k-1}\). In its simplest form, this update is given by \(S^{k}=\phi(W_x\,X^k+W_s\,S^{k-1}+b)\), where \(\phi(\cdot)\) is a nonlinear activation function, \(W_x\) and \(W_s\) are weight matrices, \(b\) is a bias term, and the output \(Y^k\) is computed directly from the hidden state \(S^k\). This basic formulation however suffers from the vanishing- and exploding-gradient problems, which impair its ability to learn long-range dependencies. To address this issue, advanced recurrent units have been developed, the most prominent being the Long Short-Term Memory (LSTM) [42] and the Gated Recurrent Unit (GRU) [43]. The GRU is adopted here, as it achieves accuracy comparable to LSTMs with fewer parameters and a simpler gating structure, making it both easier to train and more efficient to deploy.

Without going into the details of the GRU cell, Fig. 3 presents the basic model architecture, which is further enhanced in Sections  3.4 and 3.5. At time step \(k\), the model takes as input the \(x\)- and \(y\)-components of the magnetic flux density \(\mathbf{B}^k\), together with the time increment \(\Delta t^k\). The latter may vary between successive iterations to enable temporal refinement during FE simulations. Both \(\mathbf{B}^k\) and \(\Delta t^k\) are scaled, concatenated, and passed through an encoder block implemented as a multilayer perceptron, which progressively increases the feature dimensionality to match the GRU hidden-state size. The GRU cell is followed by a decoder block, also implemented as a multilayer perceptron, which maps the hidden state back to the \(x\)- and \(y\)-components of \(\tilde{\mathbf{H}}^k\), that is the predicted magnetic field, intended to approximate the target magnetic field \(\mathbf{H}^k\). Rather than initializing the hidden state to zero, it is initialized by a trainable vector, plus the output of a shallow network taking \(\mathbf{B}^0\) and \(\Delta t^0\) as inputs. This offloads the initial magnetization response, which is history-independent and lies on the virgin curve (see Section 2.2), to a dedicated subnetwork. The GRU cell therefore starts with a meaningful hidden state, and can devote its full capacity to learning history-dependent evolution, reducing the number of parameters that would otherwise be required.

3.2 Generation of the Training Datasets↩︎

In order to propose an NN-based surrogate model readily usable in a wide range of applications, the training dataset should be populated with a broad variety of field profiles. Since using data from simulations or measurements on existing devices may bias the model toward specific applications, the training \(\mathbf{H}\) sequences are generated artificially to span the diversity of waveforms encountered in electrical machines. The corresponding \(\mathbf{B}\) sequences are later obtained by feeding them to the lamination model of Section 2. As an initial strategy, the training \(\mathbf{H}\) sequences can be expressed as a sum of harmonics in which each component is defined by two orthogonal contributions with random amplitudes and phases:

\[\mathbf{H}(t)=\sum_{h\in\mathcal{H}}\frac{\alpha_h H^\star}{\sqrt{h}} \left[ \cos(\beta_h)\sin(2\pi fht+\varphi_h^x)\,\hat{\mathbf{x}} +\sin(\beta_h)\sin(2\pi fht+\varphi_h^y)\,\hat{\mathbf{y}} \right]. \label{eq:dataBasic}\tag{4}\]

In this expression, the set \(\mathcal{H}\) is generated by first randomly selecting the maximum harmonic \(h_\text{max}\in[1,50]\), later randomly selecting the number of harmonics \(|\mathcal{H}|\in[1,h_\text{max}]\), and finally randomly selecting the remaining harmonic indices without replacement in \(\{1,\ldots,h_\text{max}-1\}\). The characteristic amplitude \(H^\star\) is drawn from a log-uniform distribution over \([10~\mathrm{A/m}, 100~\mathrm{kA/m}]\), while the fundamental frequency \(f\) is uniformly sampled between 0.1 Hz and 1 kHz, with harmonic \(h\) having frequency \(h\,f\). The selected ranges for \(H^\star\) and \(h_\text{max}\) are justified in Sections 3.5 and 3.6, respectively, where the associated anhysteretic extension and resampling strategy are introduced. The coefficients \(\alpha_h \in [0,1]\) are also uniformly sampled and, together with the decay factor \(1/\sqrt{h}\), modulate the relative amplitude of each harmonic, thereby controlling the harmonic distribution and enabling the generation of both smooth and highly distorted waveforms. The parameters \(\beta_h \in [0,2\pi]\) in turn govern how each harmonic contribution is shared between the orthogonal directions \(\hat{\mathbf{x}}\) and \(\hat{\mathbf{y}}\), thus defining its orientation in the plane. Eventually, the phase terms \(\varphi_h^x\) and \(\varphi_h^y \in [0,2\pi]\) independently shift the temporal evolution of each directional contribution, while their difference determines the trajectory of each harmonic component in the plane: identical phases lead to unidirectional oscillations, whereas a non-zero phase difference produces a rotating component, resulting in circular or elliptical trajectories.

Up to now, (4 ) describes a broad class of excitations, encompassing alternating, rotating, and distorted magnetic field trajectories. However, it poorly accounts for several features commonly encountered in electrical machines: i) purely unidirectional fields, ii) DC-biased fields, iii) pulse-like fields, iv) and ramp-up fields. To improve the representation of these excitation types, several independent modifications are introduced. Each modification is activated with a probability of 50 %, avoiding both under-representation of the corresponding waveform characteristics and excessive bias toward them.

3.2.0.1 i) Unidirectional Fields

Without specific constraints on the parameters \(\beta_h\), \(\varphi_h^x\), and \(\varphi_h^y\), (4 ) almost never generates purely unidirectional fields. Such fields arise only when all \(\beta_h\) are identical and when \(\varphi_h^x\) and \(\varphi_h^y\) are either equal or differ by an integer multiple of \(\pi\) for all \(h \in \mathcal{H}\). However, unidirectional fields are commonly encountered in electrical machines. To ensure their adequate representation in the training dataset, all \(\beta_h\) are assigned the same value with a probability of 50 %. Independently, the phase terms are generated such that \(\varphi_h^x = \varphi_h^y\) for all \(h \in \mathcal{H}\), again with a probability of 50 %. Since these two conditions are imposed independently, approximately 25 % of the generated sequences are purely unidirectional.

3.2.0.2 ii) DC-biased Fields

Among other environments, DC-biased fields are encountered in the rotor of synchronous machines, in the stator of DC machines and in the presence of permanent magnets. With a probability of 50 %, (4 ) is therefore augmented by a DC component \[\mathbf{H}_\text{dc}=H_\text{dc}\,\Big(\cos\left(\beta_\text{dc}\right)\,\hat{\mathbf{x}}+\sin\left(\beta_\text{dc}\right)\,\hat{\mathbf{y}}\Big)\] where both the amplitude \(H_\text{dc}\in[0,H^\star]\) and direction \(\beta_\text{dc}\in[0,2\pi]\) are uniformly sampled.

3.2.0.3 iii) Pulse-Like Fields

In rotating electrical machines, the slotted geometry of the stator and rotor often causes the air-gap reluctance to vary periodically. As a result, the field near the air-gap takes a pulse-like shape: high when rotor and stator teeth are facing each other, and significantly reduced when a tooth faces a slot. In principle, (4 ) can build such pulse-like fields, but it is unlikely. To ensure adequate coverage of this field shape in the dataset, a pulse-like modulation is applied with a probability of 50 %. When applied, the \(x\) and \(y\) components of both (4 ) and \(\mathbf{H}_\text{dc}\) are scaled by \(\alpha_p^x + (1-\alpha_p^x)\,p(t)\) and \(\alpha_p^y + (1-\alpha_p^y)\,p(t)\), respectively, where \(\alpha_p^x\) and \(\alpha_p^y \in [0,1]\) are uniformly sampled, and \[p(t)=\sin(2\pi fh_pt+\varphi_p)^{e_p}\] is a pulse-like function. The number of pulses per period \(h_p\in[1,20]\) and the phase \(\varphi_p\in[0,2\pi]\) are both drawn from uniform distributions. The upper bound on \(h_p\) was chosen to ensure the resulting pulses are sufficiently resolved, while still allowing variability within each pulse owing to the harmonic content of (4 ). On the other hand, the exponent \(e_p\), sampled as a random integer between 3 and 20, controls the sharpness and symmetry of the pulses. When \(e_p\) is odd, \(p(t)\) preserves the sign of the sine term, resulting in alternating positive and negative pulses, as can be observed in the stator teeth. Conversely, when \(e_p\) is even, the function remains non-negative, producing pulses of constant sign, as is typical in the rotor near the air gap, a behavior that will be illustrated in Fig. 17 at the end of Section 4.2.

3.2.0.4 iv) Ramp-up Transients

The magnetic excitation can progressively build up from an initially demagnetized state. To account for such excitations, the generated field is multiplied, with a probability of 50 %, by the ramp-up function \[r(t)=1-\exp(-\gamma\pi ft)\] where \(\gamma\) is drawn from a log-uniform distribution over \([1, 10]\). This parameter controls the rise speed of the excitation, ranging from smooth progressive growth to rapid establishment of the underlying waveform.

From this procedure, \(100\,000\) \(\mathbf{H}\) sequences are generated and fed to the lamination model to form the training dataset. As will be shown in Section 3.8, this dataset size provides sufficient accuracy. Smaller datasets, on the other hand, lead to a degradation in model performance. While the lamination model is evaluated with a fine temporal discretization, as discussed in Section 2.2, the resulting sequences are downsampled to 500 uniform time steps spanning half of the fundamental period \(T=1/f\). Although the second half is not strictly redundant — owing to phase differences, pulse modulation, and ramp-up transients — it contains largely similar information. Furthermore, restricting the sequence length speeds up the RNN training while reducing the risk of vanishing or exploding gradients. Validation and test datasets are generated in the same manner, each comprising \(10\,000\) sequences. Overall, the dataset generation process requires approximately four days on a single CPU core of an HPE ProLiant XL225n compute node from the Lucia cluster. The workload is, however, embarrassingly parallel and can therefore be completed in a much shorter wall-clock time when distributed across multiple CPU cores.

3.3 Permeability-Based Error Scaling↩︎

A common loss for regression tasks is the Mean Squared Error (MSE). Restraining its definition to the time step \(k\) of a single sequence, it is given by

\[\mathcal{L}_\text{MSE} = \left|\mathbf{H}^k-\tilde{\mathbf{H}}^k\right|^2.\]

Poor results are however obtained using this loss as is. Indeed, in the saturated regime, the magnetic response can vary over tens of kA/m or more, so even a small relative error produces a large absolute contribution to the loss. Conversely, low-field regions (where the material behavior is most sensitive and non-linear) contribute comparatively little to the total error, leading to poor accuracy in this regime. Common strategies to improve accuracy across a wide range of output values include: (i) using a relative loss, e.g. \(|\mathbf{H}^k-\tilde{\mathbf{H}}^k|^2/(|\mathbf{H}^k|^2+\epsilonup)\), with an \(\epsilonup\) added in the denominator term to control the trade-off between accuracy at low and high fields; (ii) defining the loss on the time derivative, although it may suffer the same issues as the MSE, and may cause the predicted trajectories to drift progressively over time; and (iii) applying a logarithmic scaling, which improves overall accuracy but does not fully resolve the imbalance between low- and high-field regions. To overcome these limitations, this paper introduces a loss that conjointly scales the predicted field \(\tilde{\mathbf{H}}^k\) and the target field \(\mathbf{H}^k\) by the anhysteretic magnetic permeability \(\mu_\text{anh}(|\mathbf{H}^k|)\). This scaling essentially expresses the error in terms of a magnetic flux density, naturally reducing the dominance of high saturated fields. The resulting loss, referred to as the Mean Squared Scaled Error (MSSE), is specifically designed for applications involving magnetic saturation and writes

\[\mathcal{L}_\text{MSSE} = \left(\mu_\text{anh}\left(|\mathbf{H}^k|\right)\left|\mathbf{H}^k-\tilde{\mathbf{H}}^k\right|\right)^2, \label{eq:MSSE}\tag{5}\]

with

\[\mu_\text{anh}\left(|\mathbf{H}^k|\right) = \mu_0 + \mathscr{L}\left(\frac{|\mathbf{H}^k|}{h_a}\right)\frac{J_a}{|\mathbf{H}^k|} + \mathscr{L}\left(\frac{|\mathbf{H}^k|}{h_b}\right)\frac{J_b}{|\mathbf{H}^k|} \label{eq:mu},\tag{6}\]

which directly stems from (2 ). Fig. 4 demonstrates the suitability of (5 ), by comparing the predictions made on two simple sequences after having trained the model of Fig. 3 on the aforementioned losses.

3.4 Learning the Prediction Confidence↩︎

Beyond predicting the material response itself, it is valuable for the NN to also indicate how much confidence it places in each of its own predictions. Nix and Weigend [44] proposed in 1994 a loss function that enables an NN to learn both the mean and standard deviation of a normal probability distribution, originally to capture the uncertainty inherent to data collection. Although this source of uncertainty is absent here – since the lamination model used to generate the training dataset is deterministic (see Sections 2 and 3.2) – a similar loss can be used to make the NN predict an expected error. Evaluating the NN previously-trained with the MSSE loss on the whole validation dataset shows that the error \(\mu_\text{anh}(|\mathbf{H}^k|)\,|\mathbf{H}^k-\tilde{\mathbf{H}}^k|\) is exponentially distributed. Following the steps presented in [45], the loss is therefore derived by taking the negative log-likelihood of an exponential distribution with scale parameter \(\varepsilon\). This yields the exponential negative log-likelihood (ENLL) loss, adopted as the final training loss and expressed here for a single time step \(k\): \[\mathcal{L}_\text{ENLL} =\log\left(\varepsilon^k\right) +\frac{\mu_\text{anh}(|\mathbf{H}^k|)\,|\mathbf{H}^k-\tilde{\mathbf{H}}^k|}{\varepsilon^k}.\] The two terms in the ENLL loss play complementary roles: the logarithmic term penalizes overly large values of \(\varepsilon^k\), whereas the second term penalizes values that underestimate the prediction error. Their balance encourages \(\varepsilon^k\) to match the actual prediction error. Consequently, the scale parameter \(\varepsilon^k\) acts as a predicted error, expressed in tesla, which varies from sequence to sequence and from time step to time step. At evaluation, \(\varepsilon^k\) is expected to grow when the model is given sequences it poorly handled during training. In practice, \(\varepsilon^k\) is also expected to grow on sequences lying outside of the training domain, but this behavior remains a tendency rather than a certainty, as the NN response to unseen inputs cannot itself be guaranteed.

To ensure stable training, \(\varepsilon^k\) must remain strictly positive. In practice, a minimum value of \(\varepsilon_\text{min}=1\) mT is enforced using a softplus-like activation function.

\[\varepsilon^k = \frac{1}{\varkappa}\ln\left[1+\exp\bigg(\varkappa\left(\varsigma o^k+1-\varepsilon_\text{min}\right)\bigg)\right] + \varepsilon_\text{min},\]

where \(o^k\) denotes the raw output of the decoder, \(\varsigma=0.1\) controls the slope of the activation function, and \(\varkappa=5\) its sharpness.

3.5 Robustness at High Fields↩︎

The model is trained on \(\mathbf{H}\) sequences with maximum amplitudes of roughly \(100\) kA/m (see Section 3.2), which, via (6 ), correspond to anhysteretic amplitudes \(|\mathbf{B}|=\mu_\text{anh}(|\mathbf{H}|)\,|\mathbf{H}|\) of approximately \(2.05\) T. The model is therefore expected to provide reliable predictions within this range. Predictions for \(\mathbf{B}\) amplitudes exceeding \(2.05\) T, however, lie outside the training domain and their accuracy cannot be guaranteed. While exceeding \(2.05\) T is uncommon in practice, excursions to higher field values may occur during the Newton–Raphson iterations of the 2D macroscopic-scale FE simulations (see Section 4). Therefore, the model must remain well-behaved even when \(|\mathbf{B}|\) exceeds 2.05 T.

Beyond \(100\) kA/m, ferromagnetic materials are fully saturated. In this regime, the material response follows a simple anhysteretic law, such as (6 ). A practical strategy is thus to detect when the model operates outside its training domain, and to replace it with an anhysteretic law. However, this approach may introduce issues near the switching point, as it may create a discontinuity. The proposed solution relies on two points: i) having the neural network predict only the deviation from an anhysteretic curve; and ii) clipping the amplitude of the input \(\mathbf{B}^k\), ensuring that the model never receives data outside its training domain. The resulting enhanced model architecture is illustrated in Fig. 5.

3.5.0.1 i) Residual formulation

The predicted field \(\tilde{\mathbf{H}}^k\) is now written as the sum of an anhysteretic contribution \(\mathbf{H}_\text{anh}^k\) and a hysteretic deviation \(\tilde{\mathbf{H}}_\text{dev}^k\), output of the decoder. The term \(\mathbf{H}_\text{anh}^k\) is computed as \(\nu_\text{anh}(|\mathbf{B}^k|)\,\mathbf{B}^k\), where \(\nu_\text{anh}(|\mathbf{B}^k|)=1/\mu_\text{anh}(|\mathbf{B}^k|)\) expresses a magnetic reluctivity, and is implemented as a piecewise cubic Hermite interpolating polynomial (PCHIP). By inverting (6 ), it is defined on a uniform partition of \(|\mathbf{B}|^2\), each interval having a width \(\Delta B^2\). The index of the appropriate spline segment is therefore directly given by \(\lfloor (B_x^2+B_y^2)/\Delta B^2\rfloor\). If the PCHIP was defined on \(\Delta B\) intervals, an additional square root would be required. Alternatively, if it was defined on \(\Delta H\) intervals, finding the corresponding spline index would require a lookup operation. For \(|\mathbf{B}^k|\ge3\) T, the PCHIP is replaced by the analytical expression \([(|\mathbf{B}^k|-B_s)/\mu_0+H_s]/|\mathbf{B}^k|\), where the parameters \(B_s\) and \(H_s\) are chosen to guarantee \(\nu_\text{anh}(|\mathbf{B}^k|)\) is continuous at \(3\) T.

3.5.0.2 ii) Input Clipping

When the material is fully saturated, \(\mathbf{H}^k(\mathbf{B}^k)\) follows the anhysteretic law, so that \(\mathbf{H}_\text{anh}^k(\mathbf{B}^k)\) alone captures the behavior of \(\tilde{\mathbf{H}}^k(\mathbf{B}^k)\) and \(\tilde{\mathbf{H}}_\text{dev}^k(\mathbf{B}^k)\) remains constant. In this situation, the magnitude of the input \(\mathbf{B}^k\) therefore has no impact on the prediction and can be clipped without impacting the model accuracy. Equivalently, the NN can be given an approximate magnetic polarization \(\tilde{\mathbf{J}}^k=[1-\mu_0\,\nu_\text{anh}(|\mathbf{B}^k|)]\,\mathbf{B}^k\) as input, since it is naturally bounded by \(J_\text{max}=J_a+J_b\approx1.95\) T (see Section 2.1). With a training dataset populated with \(\mathbf{B}\) sequences whose maximum amplitude reaches \(B_\text{max}\approx 2.05\) T, the NN is trained on magnetic polarizations up to \([1-\mu_0\,\nu_\text{anh}(B_\text{max})]\,B_\text{max}\approx 1.92\) T, which is slightly smaller, yet sufficiently close to \(J_\text{max}=1.95\) T. Note that extending the training dataset beyond \(100\) kA/m is therefore unnecessary. It would also be counterproductive, as it would increase the cost of data generation and training, despite the saturated regime being fully described by \(\mathbf{H}_\text{anh}^k\).

In addition to ensuring robustness at high fields, the residual formulation offers further advantages related to the learning problem itself. By constraining the network to predict only the hysteretic deviation \(\tilde{\mathbf{H}}_\text{dev}^k\) from the known anhysteretic baseline \(\mathbf{H}_\text{anh}^k\), the target function is significantly smoother and smaller in magnitude than the full material response. This simplification reduces the complexity of the function the NN must approximate, allowing good accuracy to be achieved with a more compact architecture. Furthermore, since the dominant, well-understood part of the behavior is handled analytically, the NN is implicitly guided toward physically meaningful predictions, which tends to reduce the risk of spurious oscillations.

3.6 Sequence Resampling for Adaptive Time-Stepping↩︎

At each training iteration, a batch of sequences is drawn at random from the \(100\,000\) sequences composing the training dataset. To encourage accurate predictions across different sequence lengths, these sequences are randomly truncated, resulting in sub-sequences of lengths between 10 and 500 uniformly-spaced time steps. These sub-sequences are then resampled to expose the model to non-uniform temporal discretizations, making the trained model compatible with adaptive time-stepping schemes.

The complete resampling algorithm is presented in 7. In brief, the algorithm first selects a random subset of points to retain from the source sub-sequence, while enforcing a maximum gap \(\bar g = 5\) between consecutive retained points. Since the highest considered harmonic (\(50^\text{th}\) harmonic) is initially sampled by 20 points per period (see Section 3.2), enforcing \(\bar{g} = 5\) guarantees that at least 4 points per period are retained after resampling, keeping it sufficiently resolved. The algorithm then inserts up to \(n_\text{ins}\) linearly interpolated points according to a skewed random distribution, allowing the interpolated points to be either smoothly spread across the source sub-sequence or concentrated within localized regions. In the present work, \(n_\text{ins}\) is chosen equal to the length of the source sub-sequence, allowing substantial variation in the time steps while keeping the resampled sequence length comparable to that of the source sub-sequence. This limits training time and mitigates vanishing and exploding gradient issues. A minimum normalized spacing \(\delta_\text{min} = 0.05\) is also enforced between consecutive output points within each segment, corresponding to 5 % of the source-sequence time step \(\Delta t\). This minimum spacing prevents the insertion of nearly coincident points that could cause numerical instabilities during training. With this spacing, up to \(1/\delta_\text{min}-1=19\) points can be inserted within a segment. However, the algorithm limits the number of inserted points per segment to \(\bar m = 10\), in order to leave sufficient room for the remaining points to be distributed with varying, non-degenerate time steps. Moreover, during training, the loss is evaluated only at the retained source points to prevent biasing the model toward the interpolation scheme. Note that Neural Ordinary Differential Equations (Neural ODEs) could inherently handle variable time steps, as Neural ODEs learn the system dynamics and integrates it over the requested time interval, rather than directly predicting the system output at a given time. In the context of constitutive surrogate models for FE simulations, this approach may however be inappropriate, since it may require multiple NN evaluations to perform a single nonlinear iteration, thereby increasing the computational cost.

3.7 Hyperparameter Tuning and Training↩︎

Fig. 6 compares the training, validation, and test metrics obtained for different architectural variants, and justifies the hereunder selected architecture. Preliminary experiments revealed that using a single GRU cell is sufficient to achieve high quality predictions. Slightly better results are observed by stacking two GRU cells but, since it doubles the evaluation time, a single GRU cell is preferred. This GRU cell uses a hidden size of 300, chosen as a compromise between accuracy and computational cost. The encoder and decoder use respectively two and three layers, as deeper architectures did not improve performance significantly, whereas shallower configurations were found less accurate. The decoder’s final layer uses a linear activation function to produce continuous and unbounded \(\tilde{\mathbf{H}}_\text{dev}^k\) values, whereas the activation functions used in the encoder and in the first two decoder layers were determined through a grid search. To guarantee a smooth Jacobian \(\partial \mathbf{H}/\partial \mathbf{B}\) and thus well-behaved Newton–Raphson iterations, the activation function search was restricted to continuously differentiable candidates: tanh, softplus, and CELU. Although different activation function combinations produced similar final results, some exhibited instabilities during training. The final architecture uses softplus activations in both encoder layers, and CELU and softplus activations in the first and second decoder layers, respectively.

The Adam optimizer was employed throughout the optimization procedure, while the learning rate and batch size were progressively adjusted, resulting in a final learning rate of \(2\cdot10^{-4}\) for a batch size of 128. Overall, the NN-based model contains approximately 670 000 parameters, corresponding to a storage size of 2.8 MB. It was trained two days on a single NVIDIA A100 GPU from the Lucia cluster using Pytorch [46]. Although training was performed in single precision (float32), inference is carried out in double precision (float64), as the round-off errors inherent to single-precision arithmetic were observed to prevent the FE simulations of Section 4 from reaching sufficiently small residuals. Importantly, no additional training or fine-tuning was performed; all subsequent results reported in this work were obtained using this single trained model.

3.8 Standalone Performance↩︎

To assess the model performance, the error \(\mu_\text{anh}(|\mathbf{H}|)\,|\mathbf{H}-\tilde{\mathbf{H}}|\) is evaluated on the test dataset described in Section 3.2, yielding a low average value of 8.78 mT. This error is defined as the norm of the difference between the target \(\mathbf{H}\) and predicted \(\tilde{\mathbf{H}}\) magnetic field vectors, and expressed in T through the scaling factor \(\mu_\text{anh}(|\mathbf{H}|)\), allowing a direct comparison with the predicted error \(\varepsilon\). Assuming the error follows an exponential distribution, the expected fractions of predictions whose observed error lies within the \(1\varepsilon\), \(2\varepsilon\), and \(3\varepsilon\) predicted error intervals are 63.2 %, 86.5 % and 95 %, respectively. In practice, the observed fractions on the test dataset are 63.5 %, 93.6 %, and 99.2 %. The coverage curve shown on the left of Fig. 7 extends this comparison to all confidence levels by plotting the observed proportion of predictions lying within the corresponding predicted error interval. A perfectly calibrated model would yield a coverage curve lying exactly on the diagonal, where the observed coverage systematically equals the nominal confidence level. Curves below the diagonal indicate overconfidence, whereas curves above the diagonal correspond to conservative error estimates. As shown in Fig. 7, the model is slightly overconfident for small confidence levels, where the predicted error intervals are somewhat too narrow, and slightly conservative at higher confidence levels. The right panel of Fig. 7 also depicts how the mean observed error and mean predicted error vary with the magnetic flux density magnitude \(|\mathbf{B}|\). In particular, both decrease as saturation is reached, reflecting the fact that the magnetic response becomes essentially anhysteretic, so that \(\mathbf{H}_\text{anh}\) alone provides an accurate prediction. In the strongly saturated regime, the observed error drops below the imposed lower bound \(\varepsilon_\text{min}\), making the predicted error conservative.

Fig. 8 compares three target \(\mathbf{H}\) sequences from the test dataset, representative of low-, medium-, and high-field regimes, with the corresponding predicted \(\tilde{\mathbf{H}}\) sequences. It also reports the observed error \(|\mathbf{H}-\tilde{\mathbf{H}}|\) and the scaled predicted error \(\varepsilon/\mu_\text{anh}(|\mathbf{H}|)\). The figure demonstrates close agreement between the target and predicted sequences across the full range of field amplitudes, including regions where small variations in \(\mathbf{B}\) produce significantly amplified responses in \(\mathbf{H}\). Moreover, \(\varepsilon/\mu_\text{anh}(|\mathbf{H}|)\) provides a good approximation of the error \(|\mathbf{H}-\tilde{\mathbf{H}}|\).

Similarly, Fig. 9 compares the \(\tilde{\mathbf{B}}\) field, obtained by inverting the NN-based model via Newton-Raphson for a given \(\mathbf{H}\) field, against the reference lamination solution \(\mathbf{B}\). The Jacobian \(\partial \tilde{\mathbf{B}} / \partial \tilde{\mathbf{H}}\) required by the scheme is obtained via automatic differentiation of the NN-based model, yielding exact derivatives at low computational cost. This formulation also naturally enables adaptive time-stepping strategies. To illustrate this capability, the imposed \(\mathbf{H}\) field is designed to switch from strongly saturated regions to low-field regions, and an adaptive time-stepping strategy, controlled by the predicted error \(\varepsilon\), is employed: the time step \(\Delta t\) is doubled when \(\varepsilon < 25\) mT and halved when \(\varepsilon > 40\) mT, with minimum and maximum time steps of \(1.25\) \(\mu\)s and \(40\) \(\mu\)s, respectively. When fully-saturated, the anhysteretic contribution \(\mathbf{H}_\text{anh}\) alone governs the model response, resulting in low predicted error \(\varepsilon\) and allowing the adaptive scheme to increase the time step. Conversely, at lower field levels, the irreversible contribution \(\tilde{\mathbf{H}}_\text{irr}\) becomes significant and the predicted error \(\varepsilon\) increases, prompting the scheme to reduce the time step. This adaptive strategy proves effective: the predicted \(\tilde{\mathbf{B}}\) is in close agreement with the reference lamination solution \(\mathbf{B}\). Fig. 9 also demonstrates the ability of the model to extend out of its training range in the saturation regime, as \(|\mathbf{B}|\) reaches 3.2 T.

4 Macroscopic-Scale Simulation↩︎

The magnetic vector potential \(\mathbf{A}\) formulation is used for 2D FE simulations of electrical machines because it reduces to a single scalar unknown (the out of plane component of the vector potential) discretized with nodal Lagrange elements (rather than edge elements), makes \(\mathrm{div\,}\mathbf{B}=0\) exact with flux lines as its contours, lets prescribed winding currents enter trivially as a source while giving flux linkage and torque directly, and handles the dominant non-conducting regions without any topological cuts [47].

Let \(\Omega\) denote the computational domain, \(\Omega_{\text{core}} \subset \Omega\) the laminated ferromagnetic core region, \(\Omega_s \subset \Omega\) the stranded conductor region and \(\mathbf{j}_s\) the imposed current density, it seeks \(\mathbf{A}\) such that

\[\Big( \mu_0^{-1}\mathrm{curl\,}\mathbf{A} , \mathrm{curl\,}\mathbf{A}' \Big) _{\Omega \backslash \Omega_{\text{core}}} +\;\Big( \mu^{-1}\mathrm{curl\,}\mathbf{A} , \mathrm{curl\,}\mathbf{A}' \Big) _{\Omega_{\text{core}}} -\;\Big( \mathbf{j}_s , \mathbf{A}' \Big)_{\Omega_s} = 0 \label{eq:nrformulation}\tag{7}\]

holds for all test function \(\mathbf{A}'\) in the appropriate function space. The coupling with the NN-based model is introduced by replacing the constitutive relation \(\mu^{-1}\mathrm{curl\,}\mathbf{A} = \mu^{-1}\mathbf{B}\) in \(\Omega_{\text{core}}\) by the NN-based prediction \(\tilde{\mathbf{H}}(\mathbf{B})\). Owing to the nonlinearity of \(\tilde{\mathbf{H}}(\mathbf{B})\), the resulting problem is solved using a Newton-Raphson scheme. At iteration \(i\), the constitutive law is linearized around \(\mathbf{B}_{i-1}\) as

\[\tilde{\mathbf{H}}(\mathbf{B}_i)\approx\tilde{\mathbf{H}}(\mathbf{B}_{i-1})+\frac{\partial\tilde{\mathbf{H}}}{\partial\mathbf{B}}(\mathbf{B}_{i-1})\left(\mathbf{B}_i-\mathbf{B}_{i-1}\right),\]

where \(\tilde{\mathbf{H}}(\mathbf{B}_{i-1})\) is evaluated in a forward pass of the NN, and the Jacobian \(\partial\tilde{\mathbf{H}}/\partial\mathbf{B}(\mathbf{B}_{i-1})\) is computed by reverse-mode automatic differentiation. Note that the GRU hidden state is kept fixed during the Newton-Raphson iterations and updated only upon convergence. Consequently, once convergence has been reached, the NN-based model must be queried one additional time to update its hidden state. Alternatively, two hidden states can be stored per element in \(\Omega_\text{core}\): the input hidden state \(S^{k-1}\), which remains fixed during the Newton-Raphson iterations, and the output hidden state \(S^k_i\), returned by the GRU at each iteration, which becomes the new input hidden state \(S^k\) once convergence is reached. The latter approach is adopted in this work, requiring the storage of \(2 \cdot 300 = 600\) double-precision numbers (4.8 kB) per element in \(\Omega_\text{core}\).

The proposed methodology is evaluated on two representative case studies, both making use of the same NN-based model trained for stacks of M235-35A steel laminations. Subsection 4.1 first considers a three-phase transformer, and subsection 4.2 then considers a synchronous reluctance motor. In both applications, the models are discretized using the free and open-source software Gmsh [48], while the magnetic vector potential field \(\mathbf{A}\) is discretized with first order elements using the free and open-source FE solver GetDP [49]. All simulations reported in this work also perform element-wise integration using a single Gauss quadrature point, and consider relative and absolute residual tolerances of \(10^{-7}\) and \(10^{-10}\), respectively. The convergence behavior of the Newton-Raphson scheme is discussed in in subsection 4.3, while subsection 4.4 analyses the computational cost.

4.1 Three Phase Transformer↩︎

For each phase \(p \in \{A,B,C\}\) of the three-phase transformer, the current density in the primary and secondary windings is imposed as \(j_p(t)=10^5\;(1-\exp(-2\pi ft))\;\sin(2\pi ft-\phi_p)\), where \(f=50\) Hz is the frequency, \(\phi_A=0^\circ\), \(\phi_B=120^\circ\) and \(\phi_C=240^\circ\) are phase lags, and the factor \((1-\exp(-2\pi ft))\) introduces a gradual startup. A homogeneous Dirichlet boundary condition is also imposed on the boundary of \(\Omega\). Fig. 10 depicts the mesh of the analyzed three-phase transformer, together with the curves of the imposed current density in each phase and the corresponding \(B_y(t)\) and \(H_y(t)\) curves at points \(a\), \(b\), and \(c\) located at the center of each leg. These curves result either from an anhysteretic simulation (dotted lines) or from a simulation using the NN-based model (solid lines), both running over 800 constant time steps for a single period. The impact of hysteresis is visible on the \(B_y\) curves which, compared to the anhysteretic case, are delayed due to remanent magnetization.

To evaluate the effect of time-step refinement on simulation accuracy, two additional simulations are performed over a single period \(T=1/f\) using the same trained NN-based model. The first considers half less (400) constant time steps, and the second uses an adaptive time step \(\Delta t^k\). At each time step \(k\), the converged solution is accepted if the maximum predicted error \(\max_{\Omega}\varepsilon^k\) is lower than \(100\) mT, it is otherwise rejected and the time step is halved. In case \(\max_{\Omega}\varepsilon^k\) is lower than \(70\) mT, the time step is doubled back, with a maximum time step of \(T/100=200~\mu\)s, and a minimum time step of \(T/1600=12.5~\mu\)s, resulting in a simulation of 210 accepted time steps. Note that this adaptive strategy is guided by the predicted error \(\varepsilon^k\) rather than by an estimate of the local truncation error, and should therefore be regarded as a heuristic approach. Nevertheless, classical error-controlled time-stepping schemes remain applicable.

To measure the error introduced by the NN-based model, the \(\mathbf{H}\) sequences (one per element in \(\Omega_\text{core}\)) obtained during the simulation with 800 time steps are injected into the lamination model. The resulting sequences, noted \(\bar{\mathbf{B}}\), serve as target and can be compared to the \(\mathbf{B}\) sequences obtained directly during the three simulations using the NN-based model. They are used to compute the mismatch \(|\mathbf{B} - \bar{\mathbf{B}}|\), which is a direct measure of the simulation accuracy. If \(|\mathbf{B} - \bar{\mathbf{B}}|\) was zero at each Gauss point and at each time step, it would mean the NN-based model is a perfect substitute to the lamination model. The same results could therefore be obtained using the lamination model as material law, using sufficiently fine temporal refinement and at a computational cost and complexity largely superior.

The first two columns of Fig. 11 are used to assess the simulation accuracy. The first column displays the spatial distribution of the maximum-in-time mismatch \(\max_t |\mathbf{B} - \bar{\mathbf{B}}|\), while the second displays the spatial distribution of the time-averaged mismatch \(\langle |\mathbf{B} - \bar{\mathbf{B}}| \rangle_t\), both computed element-wise over the ferromagnetic core region \(\Omega_{\text{core}}\), and for the three aforementioned NN-based simulations. For all three simulations, the maximum-in-time mismatch remains below 66 mT throughout \(\Omega_{\text{core}}\), while the mean mismatch \(\langle |\mathbf{B} - \bar{\mathbf{B}}| \rangle_{t,\Omega}\), averaged over time and space, is about 6 mT. When compared locally with the norm of \(\mathbf{B}\), these mismatches correspond to relative values \(\langle |\mathbf{B} - \bar{\mathbf{B}}| \rangle_t / \langle |\mathbf{B}| \rangle_t\) that remain below 5.4 % throughout \(\Omega_{\text{core}}\), with a spatial average \(\langle \langle |\mathbf{B} - \bar{\mathbf{B}}| \rangle_t / \langle |\mathbf{B}| \rangle_t \rangle_{\Omega}\) of about 0.7 %.

Without access to the lamination model, the mismatch \(|\mathbf{B} - \bar{\mathbf{B}}|\) cannot be computed, and one should rely exclusively on the predicted error \(\varepsilon\). The third column of Fig. 11 shows the spatial distribution of its time average \(\langle\varepsilon\rangle_t\). It is observed that, except in the vicinity of point \(c\), \(\langle\varepsilon\rangle_t\) remains larger than the time-averaged mismatch \(\langle|\mathbf{B}-\bar{\mathbf{B}}|\rangle_t\). Furthermore, locations where the time-averaged mismatch is larger correspond to locations where the time-averaged predicted error is larger, showing the ability of the NN to detect challenging regions.

The time-averaged mismatch \(\langle|\mathbf{B}-\bar{\mathbf{B}}|\rangle_t\) and time-averaged predicted error \(\langle\varepsilon\rangle_t\) however provide only aggregate measures over time and do not indicate whether the predicted error consistently bounds the actual mismatch, or if the mismatch frequently exceeds the predicted error. The last column of Fig. 11 reports the percentage of time the mismatch is larger than the predicted error, hereafter referred as the exceedance \(\langle\mathbf{1}(|\mathbf{B}-\bar{\mathbf{B}}|>\varepsilon)\rangle_t\). The mean exceedance \(\langle\mathbf{1}(|\mathbf{B}-\bar{\mathbf{B}}|>\varepsilon)\rangle_{t,\Omega}\) is obtained by further averaging over space. It reaches 16.4 %, 13.6 % and 7.8 % for respectively the 800-time-step, 400-time-step and variable-time-step simulations, showing the predicted error indeed remains larger than the mismatch most of the time. The mean exceedance decreases to 4.9 %, 4 % and 2.2 % when considering twice the predicted error, and to 2 %, 1.9 % and 0.9 % when considering three times the predicted error.

While the previous analysis focuses on global and time-averaged quantities, it is also instructive to examine the behavior at specific points. Fig. 12 shows the time evolution of \(\mathbf{B}\), \(\varepsilon\), and \(|\mathbf{B}-\bar{\mathbf{B}}|/\varepsilon\), together with the corresponding \(\mathbf{B}(\mathbf{H})\) curves, at points \(c\), \(d\), and \(e\), identified in Fig. 10 and 11. Point \(c\) stands in the middle of the region where the time-averaged mismatch is greater than the time-averaged predicted error, while point \(d\) exhibits the worst maximum-in-time mismatch and point \(e\) showcases the largest time-averaged predicted error. Primarily focusing on the \(\mathbf{B}\) and \(\mathbf{B}(\mathbf{H})\) curves, it is observed that the NN curves closely match the target curves. The \(B_y\) curve of point \(c\) also shows that the higher time-average mismatch does not originate from a discrepancy in shape, but rather results from a slight amplitude offset between the target and prediction curves. This curve explains why the time-averaged mismatch is slightly larger than the time-averaged predicted error in the vicinity of point \(c\). Despite this deviation, the predicted \(B_y(H_y)\) loop still accurately reproduces the reference curve. By comparing the \(\varepsilon\) and \(\mathbf{B}\) curves, it appears that \(\varepsilon\) is minimum when \(|\mathbf{B}|\) is sufficiently large to reach saturation, in which case the model prediction mainly relies on the anhysteretic contribution \(\mathbf{H}_\text{anh}\). Fig. 12 also depicts the evolution of the time step for the adaptive-time-step simulation, and the \(B_y\) curves at the center of each of the transformer leg. The time step is decreased whenever a transformer leg switches from nearly-zero field to saturation or inversely, reflecting a rise in the predicted error during these transitions. Eventually, the \(|\mathbf{B}-\bar{\mathbf{B}}|/\varepsilon\) curves reveal that the previously reported exceedance translates locally into ratios that stay close to unity for most of the time horizon.

4.2 Synchronous Reluctance Motor↩︎

The considered synchronous reluctance motor is supplied with three-phase sinusoidal currents of \(217.31~A_{\mathrm{RMS}}\), delivering a power of \(220\) kW. It has two pole pairs and operates at a frequency of \(65\) Hz, resulting in a nominal speed of \(\dot{\,\theta_n}=1950\) rpm. The flux barriers in the rotor result from a shape optimization process aimed at minimizing torque ripple [50]. Fig. 13 presents the considered mesh — which relies on a moving band technique to handle the relative motion between stator and rotor — together with a map of the time-averaged magnetic-flux-density amplitude \(\langle|\mathbf{B}|\rangle_t\). It also introduces two coordinate systems: a fixed \(xy\) frame and a rotating \(x'y'\) frame attached to the rotor. When evaluating the NN-based model, the rotor sequences must be expressed according to the \(x'y'\) frame, in order to eliminate artificial time variations due to rotation.

As with the transformer, the target \(\bar{\mathbf{B}}\) sequences (one per element in \(\Omega_{\text{core}}\)) are obtained by first running an NN-based simulation, and subsequently feeding the resulting \(\mathbf{H}\) sequences to the lamination model. This simulation uses 4800 constant time steps to model a single period, a deliberately large number to ensure a sufficiently fine temporal resolution and accurately capture rapid variations due to rotor–stator tooth alignment. Simulations with a coarser temporal resolution may lack precision in the predicted \(\mathbf{H}\) sequences and, although the impact on subsequent NN predictions remains limited, it could drive the EB hysteresis model cells into a spurious state from which they cannot automatically recover, ultimately compromising the accuracy of the target \(\bar{\mathbf{B}}\) sequences.

At first, these \(\bar{\mathbf{B}}\) sequences are compared against the \(\mathbf{B}\) sequences obtained from an NN-based simulation with 1200 time steps for a single period. Unlike the simulations of the transformer, the three-phase currents are imposed directly as sinusoidal waveforms at their nominal amplitude, without any ramp-up phase. As a result, the first time step is more challenging to converge, since the solution must evolve from an initially completely demagnetized state along the initial magnetization curve. To facilitate convergence, a relaxation factor of 0.8 is applied for this initial time step, and no strict limit is imposed on its number of iterations.

Fig. 14 presents the spatial distributions of the

maximum-in-time mismatch \(\max_t|\mathbf{B}-\bar{\mathbf{B}}|\), time-averaged mismatch \(\langle|\mathbf{B}-\bar{\mathbf{B}}|\rangle_t\), time-averaged predicted error \(\langle\varepsilon\rangle_t\), and exceedance \(\langle\mathbf{1}(|\mathbf{B}-\bar{\mathbf{B}}|>\varepsilon)\rangle_t\). The peak mismatch \(\max_{t,\Omega} |\mathbf{B} - \bar{\mathbf{B}}|\) reaches 73 mT at point \(c\), while the time-averaged mismatch remains below 27.5 mT throughout \(\Omega_{\text{core}}\), reaching its maximum at point \(e\). The time-averaged predicted error attains its maximum at the rotor air-gap interface (point \(f\)), where the field undergoes the most abrupt variations, ranging from strong saturation when aligned with a stator tooth, to nearly zero otherwise. The predicted error also remains strictly above 1 mT, as discussed in Section 3.4, and, except for a region near point \(a\), consistently bounds the time-averaged mismatch. The mean relative mismatch is here computed as \(\langle \langle |\mathbf{B} - \bar{\mathbf{B}}| \rangle_t / (\langle |\mathbf{B}| \rangle_t + \epsilon) \rangle_{\Omega}\), where \(\epsilon = 1\) mT is introduced to avoid near-zero divisions in regions where \(\langle|\mathbf{B}|\rangle_t\) vanishes (see Fig. 13). The resulting mean relative mismatch is 0.5 %.

In Fig. 15 and 16, the \(\mathbf{B}\), \(\varepsilon\), \(|\mathbf{B}-\bar{\mathbf{B}}|/\varepsilon\), and \(\mathbf{B}(\mathbf{H})\) curves are displayed at points \(a\), \(b\) and \(c\) in the stator, and at points \(d\), \(e\) and \(f\) in the rotor. Similarly to point \(c\) in the transformer, the \(\mathbf{B}\) and \(\bar{\mathbf{B}}\) curves at point \(a\) exhibit a slight amplitude shift, which explains the localized larger error near point \(a\). For all the other curves, \(\mathbf{B}\) and \(\bar{\mathbf{B}}\) are in close agreement.

To assess the sensitivity of the method to reduced temporal resolution, three additional NN-based simulations are performed with 600, 300, and 150 nominal time steps per period. Note that, in practice, an adaptive time-stepping strategy is employed: except for the first time step, the time step is halved whenever convergence is not achieved within 20 iterations, and doubled when convergence is reached in fewer than 10 iterations. Table 1 reports the peak mismatch \(\max_{t,\Omega} |\mathbf{B} - \bar{\mathbf{B}}|\), mean mismatch \(\langle |\mathbf{B} - \bar{\mathbf{B}}| \rangle_{t,\Omega}\), and mean predicted error \(\langle\varepsilon\rangle_{t,\Omega}\) at the stator and at the rotor for the aforementioned NN-based simulations. It is first worth noting that the simulations maintain good accuracy on average, as the mean stator and rotor mismatches increase only modestly, from 6.1 to 6.9 mT and 1.8 to 2.6 mT, respectively, when reducing the number of nominal time steps from 1200 to 150. Moreover, the mean and peak mismatches remain nearly unchanged between 1200 and 300 nominal time steps, suggesting that a resolution of 300 time steps is already sufficient. This observation is further supported by the behavior of the mean predicted error, which increases only slightly from 1200 to 300 time steps, but rises more markedly from 300 to 150 time steps.

As the temporal resolution is reduced, the peak error increases more markedly at the rotor than at the stator, which can be explained by the coarser temporal discretization of the rotor sequences. Indeed, with 48 stator teeth and two pole pairs, the \(24^\text{th}\) harmonic dominates at the rotor and is sampled by only about 6 time steps in the simulation with 150 nominal time steps per period. Nevertheless, Fig. 17, which shows the \(\mathbf{B}\) and \(\mathbf{H}\) fields at point \(f\) in the rotor for the 1200, 600, 300, and 150 nominal-time-step simulations, demonstrates the ability of the NN-based model to handle such poorly resolved temporal variations. This robustness is particularly valuable during the design phase of electrical machines, where coarse time discretizations are often preferred to reduce computational costs. In this regime, the model strongly benefits from the anhysteretic contribution \(\mathbf{H}_\text{anh}\), which provides a physically consistent base to the NN predictions, regardless of the temporal discretization. Fig. 18 further compares these NN-based simulations by depicting the \(\mathbf{B}(\mathbf{H})\) curves at point \(b\) in the stator. Owing to the contribution of \(\mathbf{H}_\text{anh}\), the NN-based model provides perfect predictions in the strongly saturated regime, even for coarse temporal discretizations. At lower field levels, although some noise is observed, even the 150-nominal-time-step simulation remains capable of capturing the minor hysteresis loops.

4.3 Convergence of the Newton-Raphson scheme↩︎

Fig. 19 presents the convergence behavior of the Newton-Raphson scheme for the synchronous reluctance motor simulations performed with 1200, 600, 300, and 150 nominal time steps, using either the NN-based surrogate model or the anhysteretic law. Except for a limited number of instances in the 150-time-step case, where the NN-based simulation reduces its time step after reaching the maximum of 20 Newton iterations, both models exhibit nearly identical convergence behavior. This close superposition indicates that the NN-based model preserves the numerical behavior of the anhysteretic model. The initial time step also requires more Newton iterations to converge since, as previously mentioned, the three-phase currents are imposed directly as sinusoidal waveforms at their nominal amplitude. Consequently, the solution must evolve from a completely demagnetized state and stabilize along the first magnetization curve.

4.4 Computational Cost↩︎

All simulations were carried out on a MacBook Pro laptop equipped with an Apple M1 Max chip, using GetDP for the FE computations [49] and Python to evaluate either the anhysteretic law or the NN-based surrogate model. In practice, the transformer simulations run in 12.7, 6.8, and 4.3 minutes for the 800-, 400-, and adaptive-time-step cases respectively, while the synchronous reluctance motor simulations run in 55.4, 34.2, 20.9 and 14.8 minutes for 1200, 600, 300, and 150 nominal time steps respectively. Tables 2 and 3 provide a detailed breakdown of these timings for the transformer and synchronous reluctance motor simulations, decomposed into three contributions: evaluation of the material law in Python, assembly and solution of the system in GetDP, and input/output (I/O) operations corresponding to data exchange between Python and GetDP. This exchange is currently implemented via text files written to and read from disk, which, while functional, accounts for approximately 20 % of the total simulation time. In future work, this file-based communication will be replaced by shared-memory communication to remove this overhead.

For both applications, the time required to evaluate the NN-based model in Python is approximately that needed to assemble and solve the system in GetDP. Compared to the anhysteretic simulations, for which the time spent in Python is negligible, and assuming that the I/O overhead would be removed through improved communication, a simulation using the NN-based model is therefore limited to roughly twice the cost of an anhysteretic simulation. Beyond these absolute timings, the computational cost scales linearly with problem size. The number of elements in \(\Omega\) increases by a factor of 2.8 from the transformer to the synchronous reluctance motor (10500 to 29600), and the GetDP assembly and solve time per system solve grows by a consistent factor of 2.8 (0.09 to 0.25 s). The same holds for the Python evaluation, where a 3.7-fold increase in \(\Omega_\text{core}\) elements (4600 to 17000) yields a 3.7-fold increase in evaluation time (0.08 to 0.28 s).

5 Conclusions and Perspectives↩︎

This paper has discussed the 2D FE modeling of electrical machines with laminated ferromagnetic cores. The results show that homogenization of the laminated ferromagnetic core behavior via a surrogate NN-based model offers key implementation and numerical advantages over direct computational homogenization: i) The automatic differentiation capabilities of NNs allow an exact evaluation of the NN-based ferromagnetic model Jacobian, eliminating the need for finite differences. Furthermore, since the proposed NN is smooth by construction, owing to its differentiable activation functions, the resulting Jacobian matrices are smooth, yielding Newton-Raphson convergence comparable to that obtained with conventional anhysteretic constitutive laws. In contrast, Jacobians obtained using direct homogenization would inherit the numerical noise of finite differences and could be affected by the non-smoothness of the EB hysteresis model. ii) The NN allows coupling an \(\mathbf{H}\)-conform mesoscopic-scale model (to use the fast vector-play EB approximation) with an efficient \(\mathbf{B}\)-conform macroscopic machine model. iii) The NN-based model inherits the accuracy of the fine temporal discretization used during dataset generation and largely preserves this accuracy at inference, even when coarser time steps are used at the macroscopic scale. Furthermore, after training on variable time steps, the NN-based model seamlessly handles variable time increments, so that macroscopic-scale adaptive time-stepping strategies are naturally supported. Achieving the same with a direct homogenization approach would require maintaining a fine temporal discretization at the mesoscopic scale throughout the simulation, thereby incurring the full computational cost of fine mesoscopic-scale solves [51] iv) From the perspective of the FE solver, the proposed constitutive model behaves as a single function call. The magnetic history required to describe hysteresis is entirely encoded in the NN hidden state, eliminating the need to explicitly manage the internal variables of an underlying mesoscopic-scale model. v) Once trained for a given lamination (material grade and sheet thickness), the proposed surrogate enables computationally efficient FE simulations. Although extending the approach to a new lamination requires generating a new dataset (approximately 4 days on a single CPU core, but embarrassingly parallel) and retraining the network (approximately 2 days), the model can subsequently be used in a wide range of FE simulations, with computational times roughly twice those of conventional anhysteretic simulations. Moreover, the predicted error returned by the NN-based model provides an immediate accuracy measure of the FE simulations.

Several directions are left for future work. The present study considers sinusoidal current excitations, whereas motors are often operated using pulse-width modulation (PWM); assessing the performance of the surrogate model under such conditions is a natural next step. The NN architecture could also be extended to predict ferromagnetic core losses and to account for temperature effects, following an approach similar to [41], where NN-based models are used for AC loss prediction in superconducting conductors. On a broader scale, the method currently assumes isotropic magnetic properties; extending it to anisotropic laminations, such as grain-oriented electrical steels, would require adapting both the dataset generation and the network architecture. A 3D extension would also be needed to accurately capture eddy-current paths near edges and geometric discontinuities.

6 Parameters of the Energy-Based Model↩︎

7 Resampling algorithm↩︎

The resampling algorithm is designed to produce sequences with irregular time steps of varying density. It first thins the source sub-sequence by retaining only a random subset of points, then re-enriches it by inserting linearly interpolated points, whose distribution can range from spatially uniform to highly localized.

Acknowledgments↩︎

Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant n^o2.5020.11 and by the Walloon Region. The present research also benefited from computational resources made available on Lucia, the Tier-1 supercomputer of the Walloon Region, infrastructure funded by the Walloon Region under the grant agreement n^o1910247.

Funding sources↩︎

This work was supported in part by the Walloon Region of Belgium under grant 8139 M&SSCoT and 9001 ELECTRON. The work of Louis Denis was supported by the F.R.S-FNRS.

Declaration of generative AI and AI-assisted technologies in the manuscript preparation process↩︎

During the preparation of this work, the authors used ChatGPT in order to revise portions of the manuscript and improve its language and readability. After using this service, the authors reviewed and edited the content as needed and take full responsibility for the content of the published article.

References↩︎