# Magnetoelastic Vector Hysteresis Modeling for Electromagnetic Devices: A Combination of a Multiscale Model With the Energy-Based Hysteresis Framework

K. Roppert, M. Kaltenbacher, L. Domenig, L. Daniel

► To cite this version:

K. Roppert, M. Kaltenbacher, L. Domenig, L. Daniel. Magnetoelastic Vector Hysteresis Modeling for Electromagnetic Devices: A Combination of a Multiscale Model With the Energy-Based Hysteresis Framework. IEEE Transactions on Magnetics, 2025, 61 (8), pp.1-14. 10.1109/TMAG.2025.3584819. hal-05288689

HAL Id: hal-05288689

https://hal.science/hal-05288689v1

Submitted on 29 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Magnetoelastic Vector Hysteresis Modeling for Electromagnetic Devices: A Combination of a Multiscale Model With the Energy-Based Hysteresis Framework

K. Roppert<fnref n="1" />,<fnref n="2" />✉, M. Kaltenbacher<fnref n="1" />, L. Domenig<fnref n="1" />, and L. Daniel<fnref n="1" />,<fnref n="3" />,<fnref n="4" />

In this work, the simplified multiscale model (SMSM) is incorporated into the energy-based (EB) quasi-static vector hysteresis model to represent the anhysteretic part of the material behavior. This approach enables the inclusion of effects such as mechanical stress, magnetostriction, material anisotropy, and crystallographic texture. By integrating the anhysteretic model into the EB framework, it becomes possible to account for dissipative effects (in our case, domain wall pinning) while utilizing detailed material information. To solve the EB model in conjunction with the SMSM, two approaches are pursued: a numerical optimization of a free energy functional and an explicit approximate variant, known as the vector play model (VPM). Both methods are compared in terms of computational performance, and the differences in results are demonstrated through the simulation of the cross section of an electric machine. Furthermore, the local as well as global behavior is investigated. It is shown that in an electrical machine configuration, the VPM provides a very satisfactory approximation to both local and global responses, together with a reduced computation time compared to the EB model. In the provided application case, it is shown that a shrink-fitting operation can lead to a 30% increase in the overall hysteresis losses.

**Index Terms**— Energy-based (EB) hysteresis, finite element (FE) modeling, multiscale model (MSM), numerical optimization, vector hysteresis, vector play model (VPM).

## I. INTRODUCTION

THE accurate simulation of modern electrical devices, such as electric machines and transformers, requires advanced material models that capture a wide range of effects beyond simple ferromagnetic saturation. This includes the impact of mechanical stress, magnetostriction, material anisotropy, the crystallographic texture, and grain morphology of materials. Standard models (material models used in commercial software) often fail to adequately account for these complex interactions, necessitating the use of more sophisticated approaches. Accurate modeling of these effects is critical in various applications, such as electrical machines and transformers, where material behavior directly influences device performance.

For instance, residual stresses in stator cores of electric machines can significantly affect magnetic properties, as demonstrated in [1]. Similarly, the mechanical effects caused by cutting electric steel sheets, commonly used in e-drives and transformers, influence the electromagnetic performance of both non-oriented [2] and grain-oriented steels [3]. Additionally, core losses can change when

shrink-fitting techniques are employed in assembly, as shown in [4], [5], and [6]. These practical examples highlight the importance of incorporating a broad range of effects into material models to ensure accurate simulation results.

This article focuses on the integration of the simplified multiscale model (SMSM) into the energy-based (EB) vector hysteresis model. The SMSM is a versatile tool capable of accurately modeling a wide variety of material behaviors, including the magnetoelastic effect, magnetostriction, and the effects of texture. On its own, the SMSM provides a comprehensive representation of reversible (anhysteretic) material behavior, including the response of materials to both mechanical and magnetic excitations (see [7], [8] for further details). On the other hand, the EB vector hysteresis (EB) model excels at capturing dissipative effects (hysteresis) and can model both isotropic and anisotropic dissipative behavior [9]. The EB model is highly flexible, allowing for increasing its precision (ability to be fit to real measurements) through the inclusion of additional material parameters called pinning forces, which can be efficiently fit using a certain probability density [10] and there exist adaptions such that rotational losses are depicted correctly [11], [12], [13]. This model considers quasi-static domain wall pinning as the main source of dissipation, neglecting, following the naming convention of Bertotti, excess losses resulting from microscopic eddy currents on the grain and sub-grain lengthscale. Furthermore, no rate dependence of the hysteretic description is used,

Received 1 November 2024; revised 5 March 2025, 13 May 2025, and 25 June 2025; accepted 26 June 2025. Date of publication 1 July 2025; date of current version 28 July 2025. Corresponding author: K. Roppert (e-mail: klaus.roppert@tugraz.at).

Color versions of one or more figures in this article are available at https://doi.org/10.1109/TMAG.2025.3584819.

Digital Object Identifier 10.1109/TMAG.2025.3584819

<footnote n="1">¹Institute of Fundamentals and Theory in Electrical Engineering, 8010 Graz, Austria</footnote>

<footnote n="2">²Institute for Accelerator Science and Electromagnetic Fields (TEMF), Technische Universität Darmstadt, 64289 Darmstadt, Germany</footnote>

<footnote n="3">³CentraleSupélec, CNRS, Laboratoire de Génie Électrique et Électronique de Paris, Université Paris-Saclay, 91192 Gif-sur-Yvette, France</footnote>

<footnote n="4">⁴CNRS, Laboratoire de Génie Électrique et Électronique de Paris, Sorbonne Université, 75252 Paris, France</footnote>
making the material model(s) used in this contribution purely rate independent. Importantly, the EB model has an efficient solution scheme, making it well-suited for device simulations. Furthermore, a simplified variant of the EB model exists, the approximate vector play model (VPM) [14], which is an explicit and faster analytical approximation, providing computational efficiency while maintaining acceptable accuracy under certain conditions, which is also shown in this work.

## A. Modeling Aspect and Main Contribution

The effects discussed above, such as magnetoelasticity, magnetostriction, and material anisotropy, are modeled using the SMSM [7], [8], a reduced version of the full multiscale model (MSM) [15]. The SMSM is integrated into the EB vector hysteresis model to provide a physics-based approach for modeling anhysteretic (reversible) behavior, while the EB model captures dissipative (irreversible) effects. This decoupling of reversible and dissipative behavior allows for the flexible and accurate representation of material behavior under various conditions. This approach replaces the classical atanh or double Langevin function, typically used to model the anhysteretic behavior in EB models (see [9], [16], [17]) and opens the field to a more physics-based approach to model effects influencing the anhysteretic behavior, using the SMSM. Similar work has been carried out in [6] where the SMSM was used as the anhysteretic function in the Jiles Atherton hysteresis model, or in [18] where the SMSM together with the explicit VPM were used.

In this work, we follow in a very similar direction by including the SMSM in the (fully optimized) EB vector hysteresis model as well as in the VPM, whereas the latter one can be found in [18]. However, the primary focus of this article is the efficient numerical inclusion of the SMSM into the EB framework, and finally incorporating both approaches (SMSM + VPM and SMSM + EB) into a finite element (FE) formulation, extending the previous work of [19] and [20]. Throughout our examples, we assume an ideal polycrystal, that is, no preferred magnetization direction, which resembles an ideally isotropic non-oriented electric steel with negligible crystallographic anisotropy. This choice stems from typical industrial-grade non-oriented steels, where any residual anisotropy is small relative to, for instance, grain-oriented steel with Goss texture. However, crystallographic anisotropy can easily be introduced via the anisotropic free energy functional, introduced later on. The only source of anisotropy considered in this work results from mechanical loading. By doing so, a comparison between the very inexpensive and fast explicit VPM approach and the solution of the fully optimized EB model can be obtained not only from the material model point of view but rather from a more practical, application-oriented simulation side.

The inclusion of the SMSM into the EB framework builds on the idea of reformulating energy minimization as the search for the angle of the irreversible field [9]. We derive and implement a suitable minimization procedure using Newton-Raphson and quasi-Newton methods, which is compared against the explicit VPM approach.

**Scope and Experimental Validation:** The present work is intentionally methodological. We demonstrate how the sim-

plified magnetoelastic MSM can be embedded within an EB vector hysteresis framework and explore the impact of shrink-fitting stress numerically. A rigorous device-level validation requiring the separation of hysteresis, eddy current, and excess losses, control of the mechanical stress state in rotating parts, and accurate voltage or torque measurements lies outside the scope of this article. Material level verifications of the SMSM + VPM have already been reported in [18].

## II. MODEL DESCRIPTION

In this section, the two primary models used in this work, the SMSM and the EB vector hysteresis model (together with its VPM approximation), are briefly introduced. The goal is to establish consistent and precise nomenclature for both models to improve clarity and reproducibility. It is important to note that both models already exist, and the focus of this article is not on further developing these models but rather on developing a way to couple them for practical use.

### A. Simplified MSM

The SMSM, as described in [8], provides a computationally efficient alternative to the full MSM [15], [21]. To avoid the complexity of the full MSM, the SMSM models the material as a fictitious single crystal, simplifying the treatment of local variations in magnetic field strength $H$ and mechanical stress tensor $\sigma$ at the grain and domain levels. This allows for a domain description based on magnetic domain families, where each domain family shares the same orientation $\alpha$.

In a domain family with orientation $\alpha$, the domain magnetization $M_\alpha$ and the domain magnetostriction strain tensor $\epsilon_\alpha^\mu$ are defined as

$$M_\alpha(\alpha) = M_s\alpha, \quad \epsilon_\alpha^\mu(\alpha) = \frac{3}{2}\lambda_s\left(\alpha \otimes \alpha - \frac{1}{3}I\right) \quad (1)$$

where $M_s$ denotes the saturation magnetization and $\lambda_s$ the maximum magnetostriction strain. Each domain is assigned a certain domain volume fraction

$$\hat{f}_{\alpha}(\mathbf{H}, \boldsymbol{\sigma}) = \frac{\exp(-A_s G_{\alpha}(\mathbf{H}, \boldsymbol{\sigma}, \boldsymbol{\alpha}))}{\int_{\mathrm{dir}(S^3)} (\exp(-A_s G_{\alpha}(\mathbf{H}, \boldsymbol{\sigma}, \boldsymbol{\alpha}))) d\boldsymbol{\alpha}} \quad (2)$$

where $A_s = \text{const.}$ is a constant related to the permeability in the low-field region and $G_\alpha$ represents the total free energy of domain family $\alpha$ (comprising anisotropy, magnetic, and magnetoelastic<fnref n="1" />) contributions and the integral $\int_{\mathrm{dir}(S^3)} (\cdot) d\boldsymbol{\alpha}$ is evaluated over all directions of the unit sphere $\mathrm{dir}(S^3)$.

The total free energy $G_\alpha$ can be given, repeating [8], as the sum of magnetic-$G_\alpha^{\mathrm{mag}}$, elastic-$G_\alpha^{\mathrm{el}}$, and anisotropic-$G_\alpha^{\mathrm{an}}$ free energy, whereas the different contributions can be defined as

$$\begin{aligned} G_{\alpha}(\mathbf{H}, \boldsymbol{\sigma}, \boldsymbol{\alpha}) &= G_{\alpha}^{\mathrm{mag}}(\mathbf{H}, \boldsymbol{\alpha}) + G_{\alpha}^{\mathrm{el}}(\boldsymbol{\sigma}, \boldsymbol{\alpha}) + G_{\alpha}^{\mathrm{an}}(\boldsymbol{\alpha}) \\ &= -\mu_0 \mathbf{H} \cdot \mathbf{M}_{\alpha} - \boldsymbol{\sigma} : \boldsymbol{\epsilon}_{\alpha}^{\mu} + G_{\alpha}^{\mathrm{an}}(\boldsymbol{\alpha}). \end{aligned} \quad (3)$$

The anisotropic free energy $G_\alpha^{\mathrm{an}}$ depends on the type of crystal (uniaxial, biaxial, and cubic anisotropies) that is used. However, if one considers a perfectly arbitrary and equally

<footnote n="1">¹Note that there are different versions of the magnetoelastic energy and magnetostrictive strain tensor, including different effects, like the non-monotonicity of the permeability over the applied stress state; see [18].</footnote>
distributed orientation of (enough) grains, there is no macroscopic anisotropy, resulting in a vanishing anisotropic free energy. In the present work (mechanically unloaded), material is considered isotropic; therefore, $G_{\alpha\alpha}^{\text{an}}$ is set to zero.

The macroscopic magnetization $\hat{M}$ and magnetostrictive strain tensor $\hat{\epsilon}^{\mu}$ are computed by weighting each domain family with its volume fraction and integrating over all directions on the sphere

$$ 
\begin{aligned}
\hat{M}(\mathbf{H}, \boldsymbol{\sigma}) &= \int_{\text{dir}(S^3)} \mathbf{M}_{\alpha}(\boldsymbol{\alpha}) \hat{f}_{\alpha}(\mathbf{H}, \boldsymbol{\sigma}, \boldsymbol{\alpha}) \, d\boldsymbol{\alpha}, \\
\hat{\epsilon}^{\mu}(\mathbf{H}, \boldsymbol{\sigma}) &= \int_{\text{dir}(S^3)} \boldsymbol{\epsilon}_{\alpha}^{\mu}(\boldsymbol{\alpha}) \hat{f}_{\alpha}(\mathbf{H}, \boldsymbol{\sigma}, \boldsymbol{\alpha}) \, d\boldsymbol{\alpha}.
\end{aligned}
\quad (4) $$

To avoid the complexity of the continuous Boltzmann distribution and integration, the continuous directions on the sphere are discretized using a tri-meshed sphere with uniformly distributed nodes (also known as an icosphere or icosahedron), as detailed in [8]. In the present work, two different versions were tested, one with a meshed sphere using 2562 directions and a 2d simplification with discretizing a circle using 360 directions. Even though the spherical version includes significantly more directions, compared to the circle version, the performance difference was negligible due to the vectorization used in the implementation of the SMSM.

The discrete domain volume fraction becomes (now it is a discrete volume fraction associated with direction $\boldsymbol{\alpha}_i$ for $i \in \mathcal{N}_S$)

$$ \hat{f}_{\alpha}(\mathbf{H}, \boldsymbol{\sigma}) \approx f_{\alpha}(\mathbf{H}, \boldsymbol{\sigma}) = \frac{\exp(-A_s G_{\alpha}(\mathbf{H}, \boldsymbol{\sigma}, \boldsymbol{\alpha}_i))}{\sum_{i \in \mathcal{N}_S} (\exp(-A_s G_{\alpha}(\mathbf{H}, \boldsymbol{\sigma}, \boldsymbol{\alpha}_i)))} \quad (5) $$

and the macroscopic quantities are approximated as

$$ 
\begin{aligned}
\hat{M}(\mathbf{H}, \boldsymbol{\sigma}) &\approx M(\mathbf{H}, \boldsymbol{\sigma}) = \sum_{i \in \mathcal{N}_S} \mathbf{M}_{\alpha}(\boldsymbol{\alpha}_i) f_{\alpha}(\mathbf{H}, \boldsymbol{\sigma}, \boldsymbol{\alpha}_i) \\
\hat{\epsilon}^{\mu}(\mathbf{H}, \boldsymbol{\sigma}) &\approx \epsilon^{\mu}(\mathbf{H}, \boldsymbol{\sigma}) = \sum_{i \in \mathcal{N}_S} \boldsymbol{\epsilon}_{\alpha}^{\mu}(\boldsymbol{\alpha}_i) f_{\alpha}(\mathbf{H}, \boldsymbol{\sigma}, \boldsymbol{\alpha}_i).
\end{aligned}
\quad (6) $$

For clarity in the algorithmic description, the SMSM is represented as the operator SMSM:$(\mathbb{R}^3, \mathbb{R}^{3 \times 3}) \rightarrow (\mathbb{R}^3, \mathbb{R}^{3 \times 3})$

$$ (\mathbf{M}, \boldsymbol{\epsilon}^{\mu}) = \text{SMSM}(\mathbf{H}, \boldsymbol{\sigma}) \quad (7) $$

where the necessary material parameters for the SMSM were omitted in this definition for a more concise notation. However, they can be stated as follows.

1) $A_s$, which is composed of the initial susceptibility $\chi_0$ of the stress-free anhysteresis curve and the saturation magnetization $M_s$, computed by $A_s = 3\chi_0/(\mu_0 M_s^2)$ (from [15]). Both $\chi_0$ and $M_s$ can either be found in the literature or can easily be extracted from measurements.

2) Magnetostriction parameter $\lambda_s$ can be found in literature or extracted from measurements via strain gauges.

3) Anisotropy constants, depending on whether a smired anisotropy from an underlying crystallographic texture or other influencing factors is considered or not (see [8] for further details).

<figure><img src="image_2.png" /><figcaption>Fig. 1. Example of an anhysteretic $M-H$ curve (saturation curve) and the corresponding free energies ($F$ denoting the Helmholtz free energy and $G$ denoting the Gibbs free energy).</figcaption></figure>

## B. EB Vector Hysteresis Model

Without delving into the thermodynamics-based derivation of the EB vector hysteresis model or exploring the various formulation approaches, for example, the variational versus the differential/inclusion approaches, as discussed in [16] and [22], we adopt the variational approach with some conceptual modifications compared to [16]. These modifications, detailed in [12], frame the problem as an incremental free energy minimization, similar to elastoplastic formulations [23], [24], [25]

$$ \arg \min_{\mathbf{H}_{\text{rev}} \in \mathcal{A}} \left\{ -\Delta G_{\mathbf{H}_{\text{rev},p}}^{\mathbf{H}_{\text{rev}}} (\mathbf{H}_{\text{rev}}) - \mu_0 \mathbf{M}_{\text{p}} \cdot (\mathbf{H}_{\text{rev}} - \mathbf{H}_{\text{rev},p}) \right\} \quad (8) $$

where the subscript $[-]_p$ denotes the state variables in the previous step, as hysteretic behavior requires information from the past. The term $\mathbf{H}_{\text{rev}}$ refers to the reversible part of the total magnetic field strength $\mathbf{H}$, while $\mathbf{M}$ represents the magnetization vector. The term $\Delta G_{\mathbf{H}_{\text{rev},p}}^{\mathbf{H}_{\text{rev}}}$ describes the change in the (reversible/anhysteretic) Gibbs free energy between the current and previous steps, that is, between the fields $\mathbf{H}_{\text{rev},p}$ and $\mathbf{H}_{\text{rev}}$. The set $\mathcal{A}$ denotes a sphere (in the case of isotropic dissipation) with radius $\kappa$, often referred to as the pinning force. If the dissipation behavior is strongly anisotropic, which is the case for, for example, oriented electric steel with Goss texture, this pinning force can be considered tensorial, with the additional cost of a more elaborate identification in several directions. The optimization problem ensures that the vector difference $\mathbf{H} - \mathbf{H}_{\text{rev}}$ remains in the sphere, meaning that the irreversible field can never be larger than $\kappa$. The expression for the functional $\mathcal{L}(\mathbf{H}_{\text{rev}})$ will be used for a more concise notation in the subsequent optimization steps.

Regarding the description of the Gibbs free energy [same quantity as in the SMSM, e.g., in (2)], this can be stated as a functional of the reversible field

$$ -\Delta G_{\mathbf{H}_{\text{rev},p}}^{\mathbf{H}_{\text{rev}}}(\mathbf{H}_{\text{rev}}) = \mu_0 \int_{\mathbf{H}_{\text{rev},p}}^{\mathbf{H}_{\text{rev}}} \mathbf{M}(\tilde{\mathbf{H}}_{\text{rev}}) \, d\tilde{\mathbf{H}}_{\text{rev}} \quad (9) $$

as illustrated in Fig. 1.

The shape of the saturation curve suggests that typical approximation functions such as arctan, coth, and tanh are commonly used for isotropic behavior in hysteresis models [26]. Other approaches involve interpolating or fitting $M-H$ (or $B-H$) curves at various excitation angles [27], [28], or using magnetic energy and co-energy expressions [29], [30]. However, such approximations can impose assumptions
on the anhysteretic curve, obscuring the relationship between the model parameters and the material's physical properties. To overcome this, the SMSM is used as the description of the anhysteretic behavior, as detailed in Section III. To simplify the initial development, we start by imagining an analytic approximation of the anhysteretic curve (such as arctan) before introducing the SMSM in Section III. For clarity, we distinguish two key cases in (8).

1) Case 1: The excitation field **H** is such that the irreversible field remains within the sphere of radius $κ$. In this scenario, no change in magnetization occurs, so no optimization is needed. The reversible field **H**<sub>rev</sub> remains unchanged, and the irreversible field is calculated as **H**<sub>irr</sub> = **H** - **H**<sub>rev</sub>.

2) Case 2: The excitation field **H** causes the irreversible field **H**<sub>irr</sub> = **H** - **H**<sub>rev</sub> to reach the sphere's boundary. In this case, the reversible field vector must be adjusted such that **H**<sub>irr</sub> stays within the sphere. This introduces the need for solving the optimization problem in (8). The constrained problem can be reformulated into an unconstrained one by recognizing that the length of the irreversible field vector is $κ$, resulting in **H**<sub>irr</sub> = $κ**e**<sub>irr</sub>$, where $**e**<sub>irr</sub>$ is the direction vector of the irreversible field. Thus, instead of solving for the reversible field, we solve for the direction of the irreversible field. This idea, first proposed in [9], leads to the following unconstrained minimization problem

$$ \arg \min_{e_{\text{irr}}} \left\{ -\Delta G_{H_{\text{rev,p}}}^{H_{\text{rev}}}(e_{\text{irr}}) - \mu_0 \mathbf{M}_p \cdot (\mathbf{H} - \kappa e_{\text{irr}} - \mathbf{H}_{\text{rev,p}}) \right\}. \quad (10) $$

The direction vector of the irreversible field can easily be parametrized for the 2-D and 3-D cases as

$$ \begin{aligned} e_{\text{irr},2\text{-D}} &= (\cos \varphi, \sin \varphi)^{\top} \\ e_{\text{irr},3\text{-D}} &= (\cos \varphi \sin \vartheta, \sin \varphi \sin \vartheta, \cos \vartheta)^{\top} \end{aligned} \quad (11) $$

with $\varphi \in [0, 2\pi]$ and $\vartheta \in [0, \pi]$. For solving this system, a Newton-Raphson algorithm can be employed since the $\Delta G_{H_{\text{rev,p}}}^{H_{\text{rev}}}$ term just includes the analytically integrated anhysteretic curve, which can easily be derived with respect to the irreversible direction angles, as described for 2-D in [9] and for 3-D in [12].

### C. Vector Play Model

Although the VPM is not new and has been previously described in the literature (with [31] introducing the vector stop model and [32] extending it to include the vector stop and vector play *hysteron*), it is briefly introduced here to ensure consistent and unified nomenclature throughout this article. The VPM serves as an explicit, simplified approximation of the fully optimized EB model, providing a fast update rule for the reversible field $**H**_{rev}$ based on the new excitation **H** and the previous reversible field $**H**_{rev,p}$

$$ H_{\text{rev}}^{\text{VPM}} = H - \kappa \frac{H - H_{\text{rev,p}}}{|H - H_{\text{rev,p}}|}. \quad (12) $$

This update rule is derived using the following simplifying assumptions.

1) **Approximation 1:** The direction of the magnetization change $((M - M_{\text{prev}})/(|M - M_{\text{prev}}|))$ is identified as the direction of the irreversible field $e_{irr} = ((H - H_{\text{rev,p}})/(|H - H_{\text{rev,p}}|))$.

2) **Approximation 2:** Replace the current reversible field $H_{\text{rev}}$ by the previous reversible field $H_{\text{rev,p}}$ to obtain $H_{\text{rev}}^{\text{VPM}}$.

While these approximations mean that the VPM does not fully solve the energy minimization problem of the EB model, it still provides a fast and computationally efficient approach. Furthermore, due to the relatively small differences between the VPM and the full EB model for typical non-oriented electrical steel materials, the VPM often serves as an excellent starting point for solving the full optimization problem in (10).

Approximation 1 constrains only the change of the magnetization $\dot{M}$, the magnetization vector **M** itself remains free to orient according to the minimization of the reversible thermodynamic potential. Consequently, the characteristic rate-independent phase shift between **H** and **M** and hence the quasi-static hysteresis loop area, is fully preserved. In the uniaxial limit, where no directional change occurs, the EB reduces exactly to the VPM, so that both models coincide.

### III. OPTIMIZATION PROCEDURE OF SMSM + EBM

Introducing the SMSM from Section II-A as a physics-based alternative to the anhysteretic approaches described in Section II-B involves replacing the direct relationship between magnetization and the reversible field $M(H_{\text{rev}})$ with the SMSM operator SMSM($H_{\text{rev}}$, $\sigma$) from (7). This change results in different independent parameters for the Gibbs free energy as a function of the reversible field (or rewritten in terms of the direction vector of the irreversible field), as well as the applied mechanical stress state

$$ \arg \min_{e_{\text{irr}}} \left\{ -\Delta G_{H_{\text{rev,p}}}^{H_{\text{rev}}}(e_{\text{irr}})(\boldsymbol{e}_{\text{irr}}, \boldsymbol{\sigma}) - \mu_0 \mathbf{M}_p \cdot (\mathbf{H} - \kappa \boldsymbol{e}_{\text{irr}} - \mathbf{H}_{\text{rev,p}}) \right\}. \quad (13) $$

This formulation complicates the integration in (9), particularly when both the magnetic field (through $e_{irr}$) and the mechanical stress state vary simultaneously. It also complicates the differentiation needed for gradient-based optimization methods. This scenario represents the most challenging case in solving the optimization problem and is referred to as *Case MagMech*.

A significant simplification arises when the mechanical stress state is considered fixed while the magnetic field varies. This situation, referred to as *Case PureMag*, greatly reduces the complexity of the numerical procedure.

#### Important Notes

1) **Convexity:** It is neither assumed nor necessary to consider the convexity of the functional $\mathcal{L}(e_{irr}, \sigma)$, meaning the solution obtained is guaranteed to be a local minimum, but not necessarily the global one. However, the numerical examples in Section V will demonstrate that the VPM provides such a good starting value that it always lies within the basin of the global minimum (for the tested material parameters, representative of typical electrical steel).
2) Reversible Mechanical Process: For the loading case with varying mechanical stress and a constant magnetic field, no dissipation is assumed to occur. This process is considered fully reversible and is captured entirely by the SMSM without the need for an additional dissipative mechanism. Readers interested in a dissipative process for mechanical loading similar to the magnetic case may refer to [18], but this is beyond the scope of this work.

3) Numerical Methods: The proposed numerical algorithms in Sections III-A and III-B are gradient-based, utilizing a quasi-Newton method, Broyden-Fletcher-Goldfarb-Shanno (BFGS) for Case MagMech and Newton-Raphson for Case PureMag. While fixed-point methods could also be used, they typically converge more slowly compared to gradient-based methods.

## A. Numerical Optimization in Case MagMech

For Case MagMech, the BFGS algorithm, a quasi-Newton method, is employed, but other algorithms, such as Davidon-Fletcher-Powell (DFP), could also be used. However, special attention must be given to an effective line search handling. Before presenting the BFGS algorithm, we introduce some key remarks and findings regarding the derivatives of the functional $\mathcal{L}(e_{\text{irr}}, \sigma)$.

1) Jacobian of $\mathcal{L}(e_{\text{irr}}, \sigma)$: Let us start with the derivative of $\mathcal{L}$ with respect to the direction of the irreversible field. In this case, the parameterization using one angle $\varphi$ resulting in a 2-D setup is used due to a more concise notation, although the extension to 3-D with an additional angle, like in Section II-B is straightforward

$$ 
\begin{aligned}
\frac{\partial \mathcal{L}(e_{\text{irr}}, \sigma)}{\partial \varphi} &= -\left( \frac{\partial \Delta G_{H_{\text{rev,p}}}^{H}(\mathbf{H}_{\text{rev}}, \boldsymbol{\sigma})}{\underbrace{\partial \mathbf{H}_{\text{rev}}}_{-\mathbf{M}}} \cdot \frac{\partial \mathbf{H}_{\text{rev}}}{\partial \varphi} \right) - \mu_0 \mathbf{M}_{\text{p}} \cdot \left( -\kappa \frac{\partial e_{\text{irr}}}{\partial \varphi} \right) \\
&= \mu_0 (\mathbf{M}_{\text{p}} - \mathbf{M}) \cdot \kappa \frac{\partial e_{\text{irr}}}{\partial \varphi} \qquad (14)
\end{aligned}
$$

where the derivative of the direction vector with respect to $\varphi$ is given by

$$ \frac{\partial e_{\text{irr}}}{\partial \varphi} = (-\sin \varphi, \cos \varphi)^{\top}. \qquad (15) $$

The derivative with respect to the mechanical stress state $\sigma$ is

$$ \frac{\partial \mathcal{L}(e_{\text{irr}}, \sigma)}{\partial \sigma} = - \frac{\partial \Delta G_{H_{\text{rev,p}}}^{H}(\mathbf{H}_{\text{rev}}, \sigma)}{\partial \sigma} = \epsilon^{\mu}. \qquad (16) $$

Thus, the Jacobian can now be stated as

$$ \mathcal{J} = \left[ \mu_0 (\mathbf{M}_{\text{p}} - \mathbf{M}) \cdot \kappa \frac{\partial e_{\text{irr}}}{\partial \varphi} \right] \qquad (17) $$

which requires evaluating the SMSM model at the current magnetization and magnetostrictive strain for the current reversible field (corresponds to the applied excitation field) and the current applied mechanical stress state ($\mathbf{M}, \epsilon^{\mu}$) = SMSM($\mathbf{H}_{\text{rev}}, \sigma$).

2) Hessian of $\mathcal{L}(e_{\text{irr}}, \sigma)$: For the Hessian $\mathcal{H}$, the derivatives of the Jacobian (17) are needed. However, obtaining the full analytical derivatives of the SMSM model is challenging for Case MagMech, so we employ the BFGS algorithm, which approximates the Hessian. The algorithm proceeds as follows.

1) Initial Guess: Obtain an initial guess for $e_{\text{irr},0}$ by evaluating one step of the VPM (12) and one step of the SMSM($\mathbf{H}_{\text{rev}}, \sigma$) with the applied $\sigma$ and $\mathbf{H}_{\text{rev}} = \mathbf{H} - \kappa e_{\text{irr}}$ to obtain $\epsilon_0^\mu$.

2) Set Initial Values: Set $\mathbf{x}_0 = (e_{\text{irr}}, \epsilon^{\mu})^{\top}$ and initialize the inverse of the Hessian to be the identity matrix $\mathcal{H}_0 = \mathbf{I}$

3) **Iterate**: For each iteration $k \to k + 1$, perform the following.

a) Obtain new search direction $\mathbf{p}_k = -\mathcal{H}_k \mathbf{J}_k$.

b) Perform a linesearch to obtain the linesearch parameter $\alpha$ (discussed later).

c) Update the solution $\mathbf{x}_{k+1} = \mathbf{x}_k + \underbrace{\alpha \mathbf{p}_k}_{s_k, \dots, \text{step}}$

d) $\mathbf{y}_k = \mathcal{J}_{k+1} - \mathcal{J}_k$ with $\mathcal{J}_k$ from (17).

e) Update the inverse Hessian using the BFGS update formula

$$ \mathcal{H}_{k+1} = \mathcal{H}_k + \frac{(\mathbf{s}_k^{\top} \mathbf{y}_k + \mathbf{y}_k^{\top} \mathcal{H}_k \mathbf{y}_k)(\mathbf{s}_k^{\top} \mathbf{s}_k)}{(\mathbf{s}_k^{\top} \mathbf{y}_k)^2} - \frac{(\mathcal{H}_k \mathbf{y}_k \mathbf{s}_k^{\top} + \mathbf{s}_k \mathbf{y}_k^{\top} \mathcal{H}_k)}{\mathbf{s}_k^{\top} \mathbf{y}_k} \qquad (18) $$

A crucial part of this process is the line search, which ideally aims at minimizing the functional of the next step $\mathcal{L}_{k+1} = \mathcal{L}(\mathbf{x}_k + \alpha \mathbf{p}_k)$ with $\alpha$ as the independent variable. Although there are many different methods to solve this problem, a robust and rather fast method is following [33], where the new functional is derived with respect to $\alpha_k$, which gives, through the application of the first order optimality condition:

$$ \frac{\partial \mathcal{L}_{k+1}}{\partial \alpha_k} = \left( \frac{\partial \mathcal{L}_{k+1}}{\partial \mathbf{x}_k} \right) \frac{\partial \mathbf{x}_k}{\partial \alpha_k} = \mathcal{J}_{k+1}^{\top} \mathcal{H}_k \mathcal{J}_k = G = 0 \qquad (19) $$

which can either be solved using Brent's method or by evaluating $G$ for a few $\alpha$ parameters between zero and one and choosing the one with the smallest resulting $G$, which is a quick approximation to the optimal linesearch parameter.

## B. Numerical Optimization in Case PureMag

In Case PureMag, where the mechanical stress state $\sigma$ is constant, an analytical derivative of the magnetization with respect to the reversible field exists, as detailed in [8, Sec. 3.2]. Consequently, it is not necessary to approximate the Hessian via quasi-Newton methods. The second derivative of $\mathcal{L}(e_{\text{irr}})$ (since stress is not a variable in this case) is given, in the 2-D case with $\varphi$ as parameterization, as

$$ \frac{\partial^2 \mathcal{L}(e_{\text{irr}}, \sigma)}{\partial \varphi^2} = \mu_0 (\mathbf{M}_{\text{p}} - \mathbf{M}) \cdot \kappa \frac{\partial^2 e_{\text{irr}}}{\partial \varphi^2} + \frac{\partial e_{\text{irr}}}{\partial \varphi} \underbrace{\frac{\partial \mathbf{M}}{\partial \mathbf{H}_{\text{rev}}}}_{\text{from [8]}} \frac{\partial e_{\text{irr}}}{\partial \varphi} \kappa^2. \qquad (20) $$
Furthermore, the Jacobian consists now only of the derivative with respect to $\varphi$, which is given in (14).

The Newton-Raphson algorithm can now be given as follows.

1) **Initial Guess:** Obtain an initial guess for $e_{irr,0}$ by evaluating one step of the VPM (12) and one step of the SMSM($H_{rev}$, $\sigma$) with the applied $\sigma$ and $H_{rev} = H - \kappa e_{irr}$, which yields the initial magnetostrictive strain $\epsilon_0^\mu$.

2) **Iterate:** For each iteration $k \to k+1$, perform the following.

a) Compute the new search direction

$$p_k = \left( \left. \frac{\partial^2 \mathcal{L}(e_{irr}, \boldsymbol{\sigma})}{\partial \varphi^2} \right|_{\varphi_k} \right)^{-1} \left. \frac{\partial \mathcal{L}(e_{irr}, \boldsymbol{\sigma})}{\partial \varphi} \right|_{\varphi_k} . \quad (21)$$

b) Perform a linesearch to determine the linesearch parameter $\alpha$ (similar to the method used in Case MagMech).

c) Update the solution for the next iteration $\varphi_{k+1} = \varphi_k + \alpha p_k$.

## IV. COMPOSITE MODE

The dissipation models described in Sections II and III cannot fully capture the real behavior of magnetic materials because they rely on a single pinning force $\kappa$. This implies that once the applied magnetic field exceeds the pinning force, the magnetization changes suddenly. However, in reality, magnetization changes occur on a much finer scale, leading to smoother transitions observed in actual hysteresis loops.

A more accurate representation of the real material behavior is achieved by introducing multiple pinning forces with varying magnitudes. This is a common approach in EB models, where the magnetization is calculated for several pinning forces. The overall magnetization is obtained through a weighted superposition of these individual contributions

$$M = \sum_{i=0}^{N} \omega^{(i)} M^{(i)} \quad (22)$$

where $N$ is the number of pinning forces and $\omega^{(i)}$ the respective weight of the $i$th force. These weights typically follow a particular distribution, as described in [10].

It is important to emphasize that these parameters pertain only to the dissipative part of the model and do not affect the anhysteretic behavior described by the SMSM, which governs the overall shape of the hysteresis loop. To avoid confusion, the parameters required for the dissipative modeling are listed below, compared to the SMSM parameters outlined in Section II-A.

1) Set of pinning forces $\underline{\kappa} := \kappa_1, \kappa_2, \dots, \kappa_N$.

2) Set of weights $\underline{\omega} := \omega_1, \omega_2, \dots, \omega_N$.

By introducing a range of pinning forces and corresponding weights, the model captures the gradual transitions in magnetization that occur in real materials, leading to a more accurate depiction of the hysteresis loop. These additional parameters make the model flexible enough to match experimental data, enabling the identification of pinning forces and weight distributions from real measurements. However, both sets of

<table><caption>TABLE I
ANHYSTERETIC MATERIAL PARAMETERS FOR THE SMSM</caption><thead><tr><th>M<sub>s</sub><br>(A/m)</th><th>λ<sub>s</sub><br>(ppm)</th><th>A<sub>s</sub><br>(m<sup>3</sup>/J)</th></tr></thead><tbody><tr><td>1.41 · 10<sup>6</sup></td><td>9.4</td><td>1.1 · 10<sup>-2</sup></td></tr></tbody></table>

parameters ($\underline{\kappa}$ and $\underline{\omega}$) are assumed to be independent of the mechanical stress state. The influence of mechanical stress is solely captured by the SMSM, representing the anhysteretic character, and not by the dissipative hysteresis part. There are approaches to include the mechanical stress also in the dissipative part, resulting in a stress-induced anisotropic dissipation potential, carried out (see [18]), but in this work, the dissipation is considered isotropic and stress independent.

Finally, the computational effort grows in principle with the number of pinning forces because each force-weight combination calls for a separate optimization. However, because the magnetization contributions for each pinning force are mutually independent, one can exploit both parallelization and vectorization to mitigate the overall cost. As a result, the runtime scales sub-linearly with the number of pinning forces in practical implementations.

## V. NUMERICAL TESTS OF THE SMSM + EB AND SMSM + VPM

In this section, we compare the fully optimized EB model using the SMSM as the anhysteretic model (denoted as SMSM + EB) with the explicit, approximative VPM also using the SMSM (denoted as SMSM + VPM). The objective is to evaluate the differences between the approximated VPM, as described in [18], and the fully optimized EB model for typical non-oriented electrical steel. The material parameters for the SMSM are listed in Table I and are based on measurements carried out in [34].

The dissipation parameters for the EB and VPM models, listed in Table II, are taken from [12]. In [12], the pinning forces were split into x- and y-directions, but for simplicity, here we use an average scalar pinning force, resulting in isotropic dissipation. These parameters were obtained for a non-oriented electric steel sheet, based on the assumed pinning force and weight distribution described by [10]. As mentioned earlier, in this work, we vary only the magnetic field while keeping the mechanical stress state fixed, allowing us to use the analytical derivative in Case PureMag from Section III-B.

### A. Example I: Rotating Increasing Amplitude

In this example, we apply a circular excitation with increasing amplitude $H = 50000t/(50\pi)(\cos t, \sin t)^\top$, designed to test both models across a wide range of possible excitations, as shown in Fig. 2. This excitation represents a combination of a rotating and a uniaxial component, which helps to highlight the differences between the models. For purely uniaxial or purely rotating excitations, no significant differences occur between the EB and VPM models (independent of the anhysteretic curve). However, this combined excitation brings out subtle differences, although they remain minor for the chosen material parameters of typical electrical steel.
<table><caption>TABLE II</caption><caption>κ AND ω VALUES</caption><thead><tr><th>κ<sup>(i)</sup><br>(A/m)</th><th>ω<sup>(i)</sup></th></tr></thead><tbody><tr><td>1e-12</td><td>0.10742</td></tr><tr><td>417.7675</td><td>0.070312</td></tr><tr><td>737.83495</td><td>0.1123</td></tr><tr><td>964.2539</td><td>0.15723</td></tr><tr><td>1134.6232</td><td>0.10547</td></tr><tr><td>1330.7816</td><td>0.086914</td></tr><tr><td>1562.82545</td><td>0.092773</td></tr><tr><td>1852.8662</td><td>0.058594</td></tr><tr><td>2192.9904</td><td>0.057617</td></tr><tr><td>2585.8123</td><td>0.037109</td></tr><tr><td>3053.7581</td><td>0.03125</td></tr><tr><td>3566.2705</td><td>0.024414</td></tr><tr><td>4070.19735</td><td>0.021484</td></tr><tr><td>4663.0516</td><td>0.012695</td></tr><tr><td>5280.0725</td><td>0.0097656</td></tr><tr><td>6992.98935</td><td>0.0058594</td></tr><tr><td>8765.8392</td><td>0.0019531</td></tr><tr><td>11056.891</td><td>0.0039062</td></tr><tr><td>14530.1853</td><td>0.0019531</td></tr><tr><td>27446.22335</td><td>0.00097656</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Excitation H-field used for testing the differences between the SMSM + EB and SMSM + VPM.</figcaption></figure>

1) *Stress-Free*: In this case, the mechanical stress tensor $\sigma$ is set to zero. Using Voigt notation, the stress tensor is represented as $\sigma_V = (\sigma_{xx}, \sigma_{yy}, \sigma_{xy}) = (0, 0, 0)$, where positive values indicate tensile stress components. As seen in Fig. 3, at most 6 iterations are required to reduce the relative incremental error $||\varphi_{k+1} - \varphi_k||_2/||\varphi_k||_2$ to the specified threshold of $1 \cdot 10^{-10}$. With the BFGS algorithm from Section III-A, the number of iterations remains approximately the same. In Figs. 3-5, the differences between the SMSM + VPM and the fully optimized SMSM + EB solutions are small (qualitatively speaking), given the material parameters and excitation.

In the following, we increase the anisotropy by introducing a non-zero (but constant) stress state for further evaluation.

2) *With Mechanical Stress*: In this case, the stress tensor is set to be $\sigma_V = (100, 0, 0)$ MPa, which drastically increases the anisotropy of the material model, clearly visible in Figs. 6-8. Although the differences between the approximated SMSM + VPM and the SMSM + EB become more evident, it is still a qualitative comparison and based on the results it seems as if the region where both models differ the most is the case when the applied field, and the magnitude of the magnetization, are small and when they suddenly change directions, which can be seen in Fig. 7. These direction changes are not only

<figure><img src="image_4.png" /><figcaption>Fig. 3. Overview over magnetization, excitation, and number of iterations for $\sigma_V = (0, 0, 0)$ MPa.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Magnetization loci for $\sigma_V = (0, 0, 0)$ MPa.</figcaption></figure>

caused by sudden changes of the excitation field but also when a switching occurs, that is, when one of the pinning forces is overcome, there is a sudden change in the magnetization, which leads to slightly unphysical behavior of the SMSM + VPM (can clearly be seen in Fig. 7 at $M_x \approx 1 \cdot 10^6$ A/m and $M_y \approx 0.1 \cdot 10^6$ A/m, so that the magnetization suddenly grows overproportionally compared to the homogeneously growing excitation field). This behavior can also be observed in Fig. 4, although much less visible because there is less inherent anisotropy due to the zero-stress tensor. The Appendix includes three additional tests with progressively increasing and then decreasing excitation amplitudes, further illustrating this behavior (see Appendix A). Notably, reducing the quasi-timestep, defined as the increment of change in the magnetic field strength input, does not reduce the discrepancies between
<figure><img src="image_6.png" /><figcaption>Fig. 5. Magnetostriction strain over magnetization for $\sigma_V$ = (0, 0, 0) MPa.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Overview over magnetization, excitation, and number of iterations for $\sigma_V$ = (100, 0, 0) MPa.</figcaption></figure>

the two models. This is because the VPM faces limitations in capturing the abrupt directional shifts that arise when a specific threshold in pinning force is surpassed, as shown in Appendix B. Furthermore, increasing the number of pinning forces to better resolve these abrupt directional changes does not result in the convergence of the SMSM + VPM model toward the SMSM + EB model either, even with 1000 pinning forces, as carried out in Appendix C.

Up to this point, both models (SMSM + EB and SMSM + VPM) have only been evaluated and compared on the local material scale, simulating what could be seen as controlled

laboratory conditions. However, the more practically relevant comparisons should occur at the device level, where these models are used in electromagnetic simulations for engineering applications. For example, we can utilize both models as material inputs within a field simulation, such as an FE method (FEM) simulation of an electric machine. This is where the true value of the models will become apparent, as the differences between the two may influence the accuracy of performance predictions for real-world devices. In the next section, we extend the model comparison from the material level to a device-scale simulation using FEM, assessing how these models perform under more complex and practically relevant conditions.

<figure><img src="image_8.png" /><figcaption>Fig. 7. Magnetization loci for $\sigma_V$ = (100, 0, 0) MPa.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Magnetostriction strain over magnetization for $\sigma_V$ = (100, 0, 0) MPa.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 9. Geometry of the simplified PMSM with $V_{1,2,3}$ as the three stator windings and $E_{1,2,3,4,5}$ marking the location of evaluation points where the magnetization loci are evaluated.</figcaption></figure>

<table><caption>TABLE III
$\kappa$ AND $\omega$ VALUES</caption><thead><tr><th>$\kappa^{(i)}$ (A/m)</th><th>$\omega^{(i)}$</th></tr></thead><tbody><tr><td>1e-12</td><td>0.25</td></tr><tr><td>46.7</td><td>0.25</td></tr><tr><td>93.3</td><td>0.25</td></tr><tr><td>140.0</td><td>0.25</td></tr></tbody></table>

## VI. USING SMSM + EB AND SMSM + VPM IN AN FE SIMULATION

To quantify the differences between the two vector hysteresis approaches, both models (SMSM + EB and SMSM + VPM) were implemented as the material model in an FE formulation using the magnetic scalar potential approach described in [19] and [20]. These simulations were carried out in the open-source FE software openCFS [35]. The simulated device is based on the cross section of a six-pole, 36-slot permanent magnet synchronous machine (PMSM), using a 2-D triangular mesh, as shown in Fig. 9. Inside the permanent magnet region, a remanent flux density of $B_r = 0.6$ T was prescribed, and, additionally, the rotor is kept fixed, without rotation. The stator and the rotor are composed of electrical steel, with the anhysteretic parameters for the SMSM provided in Table I. Due to the nature of the device considered (laminated steel sheets), no electric conductivity was assigned to any parts of the computational domain (no macroscopic eddy currents).

Regarding the hysteretic parameters (set of $\kappa$ pinning parameters), only 4 pinning forces and weights were used due to a faster numerical simulation, given in Table III, which are uniformly distributed, depicting a purely isotropic dissipative behavior. This is in contrast to the previously used set of $\kappa$ and $\omega$ values in Table II (the weights and pinning forces do not follow the distribution described in [10]). In this case, the weights chosen were equally distributed, and the $\kappa$ parameters were obtained by fitting to uniaxial hysteresis

loop measurements of a non-oriented electric steel. These dissipative parameters were used for both subsequent SMSM + VPM as well as SMSM + EB simulations. The air gap and stator windings are modeled as air with vacuum permeability $\mu_0$, and a three-phase excitation is applied by prescribing the current density in the stator windings with a frequency of $f = 50$ Hz. This includes a start-up phase where the current amplitude is linearly increased. The current densities in the three phases are defined as

$$J_i(t) = \begin{cases} t J_i \sin(2\pi f t + \phi_i), & \text{for } t \le 0.005 \\ J_i \sin(2\pi f t + \phi_i), & \text{if } t > 0.005 \end{cases}$$

where $\phi_i \in \{0^\circ, 120^\circ, 240^\circ\}$ and $J_i = 3.468025 \cdot 10^6$ A/m$^2$. The simulation proceeds with a pseudo-time step of $\Delta t = 0.001$ seconds, and it is run for 300 time steps. Note that since eddy currents and rate-dependent effects are excluded, the excitation frequency does not influence the results.

The goal of this numerical experiment is to quantify the impact of choosing between the SMSM + EB and SMSM + VPM models at the device level, extending the material-scale comparisons presented in Section V. The comparison includes both local magnetization loci at different locations in the stator and rotor, as well as global hysteresis losses by calculating the area of the $B_x-H_x$ and $B_y-H_y$ loops (as described in [36]) over one cycle (one 360° rotation of the stator field), once the steady-state is reached. The energy loss is computed element-wise and then summed over all stator and rotor elements.

Additionally, to assess the significance of hysteresis models (SMSM + EB and SMSM + VPM) versus purely anhysteretic models, the simulations are also performed with and without hysteresis, as well as with and without a prescribed mechanical stress state. For the stress state, a mechanical simulation was conducted to model a press fit of the rotor and stator sheet onto a rod or casing by prescribing a radial displacement at the rotor's inner and stator's outer edge, resulting in the stress distribution shown in Fig. 10, with a maximum von Mises stress of 67 MPa.

To ensure a fair comparison of computational performance between the SMSM + VPM and SMSM + EB material models in the FEM simulation, all non-linear simulations were iterated until the $L_2$ norm of the residual reached $1 \cdot 10^{-8}$. The non-linear iteration scheme used is detailed in [19] and [20], as previously mentioned.

Without the prescribed stress state from Fig. 10, the number of global iterations required by the FEM model was approximately the same for both SMSM + VPM and SMSM + EB. However, due to the additional iterations needed for the underlying material optimization in the SMSM + EB model, the wall clock time was roughly four times higher compared to the SMSM + VPM model (see Table IV). To provide an unbiased performance comparison, all FEM simulations were carried out without parallelization on a single CPU core and a single thread of an apple m3 pro chip. For comparison, also the purely anhysteretic (without hysteresis) wall clock time was added in Table IV. When including the prescribed stress state, the SMSM + VPM model's iteration count remained
<figure><img src="image_11.png" /><figcaption>Fig. 10. Mechanical stress (in Pa) after prescribing a purely radial displacement at the rotor's inner and stator's outer edge. (a) von Mises stress. (b) $\sigma_{rr}$. (c) $\sigma_{\phi\phi}$. (d) $\sigma_{\tau\phi}$.</figcaption></figure>

<table><caption>TABLE IV
COMPUTATION TIMES FOR 300 TIMESTEPS WITH THE VARIOUS MODELS UNDER NO STRESS AND WITH STRESS CONDITIONS</caption><thead><tr><th></th><th>SMSM<br>(s)</th><th>SMSM+VPM<br>(s)</th><th>SMSM+EB<br>(s)</th></tr></thead><tbody><tr><td>No Stress</td><td>5745</td><td>7820</td><td>30473</td></tr><tr><td>With Stress</td><td>5862</td><td>7889</td><td>51859</td></tr></tbody></table>

largely unaffected, whereas the SMSM + EB model showed a significant increase in the number of iterations, leading to a wall clock time approximately six to seven times higher than that of the SMSM + VPM simulation.

## A. Local and Global Comparisons

Figs. 11–15 show the local magnetization loci comparisons between SMSM + VPM and SMSM + EB, both with and without mechanical stress, at evaluation points $E_1 - E_5$ as marked in Fig. 9. Similar to the earlier local-scale comparisons, the differences between the SMSM + VPM and SMSM + EB models are small. Qualitatively, they seem negligible under these conditions. Only at evaluation point $E_4$, a more significant difference between both models can be observed in the case with non-zero stress (see Fig. 14), since the magnetization values are far from saturation and in this regime, the influences of mechanical stress are the highest.

However, the impact of considering hysteresis is more pronounced when comparing the hysteresis models with the purely anhysteretic SMSM, as seen in the first row of Figs. 11–15. The inclusion of dissipation significantly alters the local magnetization behavior, highlighting the importance

<figure><img src="image_12.png" /><figcaption>Fig. 11. Evaluation point E1: Comparison of both vector hysteresis models and the pure SMSM with and without mechanical stress. (a) With mech. stress. (b) Zero mech. stress. (c) SMSM + VPM comparison. (d) SMSM + EB comparison.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Evaluation point E2: Comparison of both vector hysteresis models and the pure SMSM with and without mechanical stress. (a) With mech. stress. (b) Zero mech. stress. (c) SMSM + VPM comparison. (d) SMSM + EB comparison.</figcaption></figure>

of considering vector hysteresis in the macroscopic simulation, at least from a local perspective.
<figure><img src="image_14.png" /><figcaption>Fig. 13. Evaluation point E3: Comparison of both vector hysteresis models and the pure SMSM with and without mechanical stress. (a) With mech. stress. (b) Zero mech. stress. (c) SMSM + VPM comparison. (d) SMSM + EB comparison.</figcaption></figure>

Furthermore, the influence of mechanical stress is evident even though the differences in the local magnetization loci seem negligible at evaluation points $E_1$ and $E_2$ in Figs. 11 and 12, respectively. The reason for this is the small stress value and the already high saturation of the material, partly caused by the permanent magnets, which introduce a constant offset (bias). In these regions, mechanical stress has less influence on the magnetization curve than in less saturated regions. This is demonstrated by plotting the magnetization loci at evaluation points $E_3$ — $E_5$, which are located in the rotor sheet closer to the shaft, right next to a stator coil an at the outer stator sheet, respectively (see Figs. 13–15). In these regions, the difference between simulations with and without the inclusion of mechanical stress becomes more evident.

Notably, these regions also seem to cause a significant change in the global level, which will be discussed next.

Until now, only local comparisons between the different models with and without dissipation and with and without stress have been discussed. To compare them on a more global scale, the local element-wise hysteresis losses (energy loss per cycle—one 360° rotation of the stator field), as previously described, were evaluated, scaled by the element volume, and summed over all elements of the stator and rotor regions. The influence of mechanical stress on the local element-wise hysteresis loss densities can also be seen in Figs. 16 and 17, where the loss density is visibly higher in the region with large compressive stress components (stator sheet press fit region),

<figure><img src="image_15.png" /><figcaption>Fig. 14. Evaluation point E4: Comparison of both vector hysteresis models and the pure SMSM with and without mechanical stress. (a) With mech. stress. (b) Zero mech. stress. (c) SMSM + VPM comparison. (d) SMSM + EB comparison.</figcaption></figure>

<table><caption>TABLE V
ENERGY LOSS PER CYCLE (ONE 360° ROTATION OF THE STATOR FIELD)</caption><thead><tr><th></th><th>SMSM+VPM<br>no stress<br>(J)</th><th>SMSM+EB<br>no stress<br>(J)</th><th>SMSM+VPM<br>with stress<br>(J)</th><th>SMSM+EB<br>with stress<br>(J)</th></tr></thead><tbody><tr><td>Energy loss</td><td>1.388</td><td>1.383</td><td>1.822</td><td>1.810</td></tr></tbody></table>

while being not as pronounced in regions with tensile stress, for example, in the rotor sheet press fit region.

The results of the integrated hysteresis loss densities over the rotor and stator sheet region are presented in Table V. The relative difference between the SMSM + VPM and SMSM + EB models is 0.38% in the zero-stress case and 0.65% in the non-zero-stress case. Although these differences appear negligible for this particular use case, they reflect the same trend observed in the local comparisons—specifically, the SMSM + VPM shows slight deviations from the SMSM + EB model in scenarios with increased anisotropy, leading to more rapid changes of the magnetization vector. However, in this case, the differences are insignificant.

Regarding the impact of mechanical stress, the SMSM + VPM predicts an increase in energy loss of 31.25%, whereas the SMSM + EB predicts a 30.90% increase when the stress state from Fig. 10 is applied.

This behavior, although still of an academically simplified nature, definitely justifies performing further, more realistic
<figure><img src="image_16.png" /><figcaption>Fig. 15. Evaluation point E5: Comparison of both vector hysteresis models and the pure SMSM with and without mechanical stress. (a) With mech. stress. (b) Zero mech. stress. (c) SMSM + VPM comparison. (d) SMSM + EB comparison.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. Local hysteresis loss density evaluated with SMSM + EB. (a) Without mechanical stress. (b) With mechanical stress.</figcaption></figure>

simulations, for example, with a rotating rotor or an additional stress state caused by the rotation of the rotor.

## VII. CONCLUSION

In this work, the SMSM was successfully integrated into the fully optimized EB vector hysteresis framework, providing a physics-based anhysteretic model capable of capturing important material effects such as magnetostriction and the magnetoelastic effect. Two efficient numerical schemes were developed: one leveraging the Newton-Raphson method for the constant-stress case and the BFGS algorithm for cases with varying mechanical stress. These schemes enable the

<figure><img src="image_18.png" /><figcaption>Fig. 17. Difference between local hysteresis loss density without stress minus with stress. Positive values indicate higher losses with stress.</figcaption></figure>

effective solution of the coupled SMSM + EB model, making it practical for device-level simulations.

The applicability of the SMSM + EB model for device simulations was demonstrated through its implementation in an FE formulation for a simplified PMSM. Despite the VPM being a simplification of the full EB model, results with and without mechanical stress showed a close match to the fully optimized EB model, but significant differences in their computational performance became obvious. This suggests that, for certain applications, the SMSM + VPM may be a more practical choice without substantial loss of accuracy.

Further investigations are necessary to confirm the sufficiency of the SMSM + VPM in more general cases. While our current use case suggests that using the SMSM + VPM in a macroscopic device-level FEM simulation provides comparable results to the full SMSM + EB model, this conclusion is based on specific assumptions such as an initially isotropic material and no rotation of the rotor. These factors need to be explored in future works to fully validate the use of the SMSM + VPM over the full SMSM + EB material model in a wider range of scenarios.

**Future Work:** Several promising extensions of this work are identified as follows.

1) *Incorporating Texture Information:* The next logical step is to integrate crystallographic texture data, for example, from electron backscatter diffraction (EBSD) measurements, using the simplified texture MSM [8]. This will allow for even more accurate modeling of material anisotropy at the grain level. A first step in this direction would be to incorporate macroscopic anisotropy in the SMSM by including the anisotropic energy in the global energy minimization of the SMSM model.

2) *More Realistic PMSM Simulations:* Performing an extensive set of FE simulations on a PMSM with included rotation of the rotor would more accurately reflect the real magnetic field distribution in an operating machine, showing the practical relevance of including mechanical stress in electromagnetic device simulations.

3) *Machine-Level Experiments and Validations:* Although this article has demonstrated feasibility and consistency
<figure><img src="image_19.png" /><figcaption>Fig. 18. Increasing and decreasing excitation.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 19. Magnetization loci for $\sigma_V$ = (100, 0, 0) MPa and increasing-decreasing excitation.</figcaption></figure>

at the material and simplified machine cross section level, a possible next step is to conduct full-scale measurements on an actual rotating electric machine. This will allow detailed comparisons of modeled vs. measured quantities such as iron losses, local flux densities, mechanical deformations, and induced voltages at various loads.

4) **Accounting for Additional Mechanical Stresses:** Including mechanical stresses due to the rotor's rotation would introduce more complex stress distributions and allow for the evaluation of these effects on the magnetization and hysteresis behavior, further improving the estimation of the model's applicability in real-world conditions.

# APPENDIX A
## INCREASING-DECREASING EXCITATION

Using the same material parameters as in Section V, an increasing and then decreasing magnetic excitation is prescribed, as displayed in Fig. 18, highlighting the region of difference between SMSM + VPM and SMSM + EB, which occurs specifically in regions when a pinning force is overcome, displayed in Fig. 19. In these regions, the permeability changes drastically, caused by fast changes of the magnetization, where the approximation of the VPM becomes obvious. A particularly illustrative behavior occurs at $M \approx [1.0 \cdot 10^6, -0.3 \cdot 10^6]$ A/m (and also at $M \approx [-0.75 \cdot 10^6, 0.2 \cdot 10^6]$ A/m), where switching events can be observed. The SMSM + EB optimization prevents the magnetization from abruptly changing. The fact that the deviations are indeed caused by these switching events is evident in Fig. 19, as the trajectories begin to realign after the final switching event,

<figure><img src="image_21.png" /><figcaption>Fig. 20. Magnetization loci for $\sigma_V$ = (100, 0, 0) MPa with 10× the number of datapoints compared to Fig. 19. However, for plotting purposes and comparability, the dataset was downsampled to the same timesteps as in the original version.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 21. Magnetization loci for $\sigma_V$ = (100, 0, 0) MPa with 1000 pinning forces.</figcaption></figure>

which occurs at approximately $M \approx [0.8 \cdot 10^6, -0.5 \cdot 10^6]$ A/m, following a period of increased deviation.

# APPENDIX B
## INCREASING TEMPORAL RESOLUTION

As already mentioned in the main text body, increasing the temporal resolution (smaller increments in the change of the magnetic field strength input) cannot lessen the differences in those regions, as displayed in Fig. 20.

# APPENDIX C
## INCREASING NUMBER OF PINNING FORCES

Increasing the number of pinning forces for a fixed time step width decreases the differences slightly, as displayed in
Fig. 21. However, even when increasing the number of pinning forces up to 1000, the differences remain visible, and further increasing the number does not decrease the difference.

REFERENCES

[1] Y. Kai, Y. Tsuchida, T. Todaka, and M. Enokizono, "Evaluation of local residual stress distribution of stator core in a rotating machine," *Electr. Eng. Jpn.*, vol. 181, no. 3, pp. 1-8, Jul. 2012.

[2] P. Baudouin, A. Belhadj, F. Breaban, A. Deffontaine, and Y. Houbaert, "Effects of laser and mechanical cutting modes on the magnetic properties of low and medium Si content nonoriented electrical steels," *IEEE Trans. Magn.*, vol. 38, no. 5, pp. 3213-3215, Sep. 2002.

[3] R. Sundaria, A. Hemeida, A. Arkkio, A. Daem, P. Sergeant, and A. Belahcen, "Effect of different cutting techniques on magnetic properties of grain oriented steel sheets and axial flux machines," in *Proc. 45th Annu. Conf. IEEE Ind. Electron. Soc.*, vol. 1, Oct. 2019, pp. 1022-1027.

[4] N. Takahashi and D. Miyagi, "Effect of stress on iron loss of motor core," in *Proc. IEEE Int. Electric Mach. Drives Conf. (IEMDC)*, May 2011, pp. 469-474.

[5] K. Yamazaki, H. Mukaiyama, and L. Daniel, "Effects of multi-axial mechanical stress on loss characteristics of electrical steel sheets and interior permanent magnet machines," *IEEE Trans. Magn.*, vol. 54, no. 3, pp. 1-4, Mar. 2018, doi: 10.1109/TMAG.2017.2757531.

[6] L. Bernard and L. Daniel, "Effect of stress on magnetic hysteresis losses in a switched reluctance motor: Application to stator and rotor shrink fitting," *IEEE Trans. Magn.*, vol. 51, no. 9, pp. 1-13, Sep. 2015.

[7] L. Daniel, O. Hubert, and M. Rekik, "A simplified 3-D constitutive law for magnetomechanical behavior," *IEEE Trans. Magn.*, vol. 51, no. 3, pp. 1-4, Mar. 2015.

[8] L. Bernard, B. J. Mailhé, N. Sadowski, N. J. Batistela, and L. Daniel, "Multiscale approaches for magneto-elasticity in device simulation," *J. Magn. Magn. Mater.*, vol. 487, Oct. 2019, Art. no. 165241.

[9] L. Prigozhin, V. Sokolovsky, J. W. Barrett, and S. E. Zirka, "On the energy-based variational model for vector magnetic hysteresis," *IEEE Trans. Magn.*, vol. 52, no. 12, pp. 1-11, Dec. 2016.

[10] R. Scorretti and F. Sixdenier, "An analytical formula to identify the parameters of the energy-based hysteresis model," *J. Magn. Magn. Mater.*, vol. 548, Apr. 2022, Art. no. 168748.

[11] A. Sauseng, L. Domenig, K. Roppert, and M. Kaltenbacher, "Adaptions of the energy-based hysteresis model for correct rotational losses," in *Proc. IEEE 20th Biennial Conf. Electromagn. Field Comput. (CEFC)*, Oct. 2022, pp. 01-02.

[12] A. Sauseng, M. Kaltenbacher, and K. Roppert, "Revisiting the dry friction-like magnetic vector hysteresis model," *J. Magn. Magn. Mater.*, vol. 604, Aug. 2024, Art. no. 172285.

[13] L. Domenig, K. Roppert, A. Gschwentner, A. Sauseng, and M. Kaltenbacher, "Identification of parameters to correctly adapt energy-based hysteresis models regarding rotational losses," *IEEE Trans. Magn.*, vol. 60, no. 3, pp. 1-4, Mar. 2024.

[14] F. Henrotte, A. Nicolet, and K. Hameyer, "An energy-based vector hysteresis model for ferromagnetic materials," *COMPEL, Int. J. Comput. Math. Electr. Electron. Eng.*, vol. 25, no. 1, pp. 71-80, Jan. 2006.

[15] L. Daniel, O. Hubert, N. Buiron, and R. Billardon, "Reversible magneto-elastic behavior: A multiscale approach," *J. Mech. Phys. Solids*, vol. 56, no. 3, pp. 1018-1042, Mar. 2008.

[16] V. Francois-Lavet, F. Henrotte, L. Stainier, L. Noels, and C. Geuzaine, "An energy-based variational model of ferromagnetic hysteresis for finite element computations," *J. Comput. Appl. Math.*, vol. 246, pp. 243-250, Jul. 2013.

[17] S. Steentjes, F. Henrotte, C. Geuzaine, and K. Hameyer, "A dynamical energy-based hysteresis model for iron loss calculation in laminated cores," *Int. J. Numer. Modelling: Electron. Netw., Devices Fields*, vol. 27, no. 3, pp. 433-443, Oct. 2013.

[18] L. G. da Silva, L. Bernard, M. Domenjoud, and L. Daniel, "A magneto-elastic vector-play model including piezomagnetic behavior," *J. Magn. Magn. Mater.*, vol. 609, Nov. 2024, Art. no. 172439.

[19] L. D. Domenig, K. Roppert, and M. Kaltenbacher, "Incorporation of a 3-D energy-based vector hysteresis model into the finite element method using a reduced scalar potential formulation," *IEEE Trans. Magn.*, vol. 60, no. 6, pp. 1-8, Jun. 2024.

[20] H. Egger, F. Engertsberger, L. Domenig, K. Roppert, and M. Kaltenbacher, "On nonlinear magnetic field solvers using local quasi-Newton updates," *Comput. Math. Appl.*, vol. 183, pp. 20-31, Apr. 2025.

[21] L. Daniel and N. Galopin, "A constitutive law for magnetostrictive materials and its application to Terfenol-D single and polycrystals," *Eur. Phys. J. Appl. Phys.*, vol. 42, no. 2, pp. 153-159, May 2008.

[22] M. Kaltenbacher, K. Roppert, L. D. Domenig, and H. Egger, "Comparison of energy based hysteresis models," in *Proc. 23rd Int. Conf. Comput. Electromagn. Fields (COMPUMAG)*, Jan. 2022, pp. 1-4.

[23] M. Canadija and J. Mosler, "On the thermomechanical coupling in finite strain plasticity theory with non-linear kinematic hardening by means of incremental energy minimization," *Int. J. Solids Struct.*, vol. 48, nos. 7-8, pp. 1120-1129, Apr. 2011.

[24] K. M. Lee, S. Y. Park, M. Y. Huh, J. S. Kim, and O. Engler, "Effect of texture and grain size on magnetic flux density and core loss in non-oriented electrical steel containing 3.15% Si," *J. Magn. Magn. Mater.*, vol. 354, pp. 324-332, Mar. 2014.

[25] A. Mielke, "Existence of minimizers in incremental elasto-plasticity with finite strains," *SIAM J. Math. Anal.*, vol. 36, no. 2, pp. 384-404, Jan. 2004.

[26] G. Mörée and M. Leijon, "Review of hysteresis models for magnetic materials," *Energies*, vol. 16, no. 9, p. 3908, May 2023.

[27] M. Enokizono and N. Soda, "Magnetic field analysis by finite element method using effective anisotropic field," *IEEE Trans. Magn.*, vol. 31, no. 3, pp. 1793-1796, May 1995.

[28] F. Martin, D. Singh, A. Belahcen, P. Rasilo, A. Haavisto, and A. Arkkio, "Analytical model for magnetic anisotropy of non-oriented steel sheets," *COMPEL-Int. J. Comput. Math. Electr. Electron. Eng.*, vol. 34, no. 5, pp. 1475-1488, Sep. 2015.

[29] P. P. Silvester and R. P. Gupta, "Effective computational models for anisotropic soft B-H curves," *IEEE Trans. Magn.*, vol. 27, no. 5, pp. 3804-3807, Sep. 1991.

[30] T. Pera, F. Ossart, and T. Waeckerle, "Field computation in non linear anisotropic sheets using the coenergy model," *IEEE Trans. Magn.*, vol. 29, no. 6, pp. 2425-2427, Nov. 1993.

[31] T. Matsuo and M. Shimasaki, "Isotropic vector hysteresis represented by superposition of stop hysteron models," *IEEE Trans. Magn.*, vol. 37, no. 5, pp. 3357-3361, May 2001.

[32] T. Matsuo, "Rotational saturation properties of isotropic vector hysteresis models using vectorized stop and play hysteresons," *IEEE Trans. Magn.*, vol. 44, no. 11, pp. 3185-3188, Nov. 2008.

[33] K. Fujiwara, Y. Okamoto, A. Kameari, and A. Ahagon, "The Newton-raphson method accelerated by using a line search-comparison between energy functional and residual minimization," *IEEE Trans. Magn.*, vol. 41, no. 5, pp. 1724-1727, May 2005.

[34] L. G. da Silva, "Energy-based models for the magneto-elastic behavior of ferromagnetic materials," Ph.D. dissertation, Dept. Laboratoire de Génie Electrique et Electronique de Paris, Universidade Federal de Santa Catarina, Florianópolis, Brazil, 2024.

[35] M. Kaltenbacher, *Numerical Simulation of Mechatronic Sensors and Actuators: Finite Elements for Computational Multiphysics*. Cham, Switzerland: Springer, 2015.

[36] G. Mörée and M. Leijon, "Iron loss models: A review of simplified models of magnetization losses in electrical machines," *J. Magn. Magn. Mater.*, vol. 609, Nov. 2024, Art. no. 172163.