1 Introduction↩︎

Real ferromagnetic materials exhibit complex hysteresis effects, requiring memory-dependent constitutive mapping between the magnetic field \(\mathbf{H}\) and the flux density \(\mathbf{B}\) for solving low-frequency Maxwell’s equations. Operator-based phenomenological approaches [1] are widely used, expressed as a vector-valued, rate-independent operator \(\mathbf{H} = \boldsymbol{\textit{V}}[\mathbf{B}]\). The generalized Prandtl–Ishlinskii model [2] effectively characterizes hysteresis nonlinearities and its stop-type formulation defines the mapping \(\boldsymbol{\textit{V}}\) through a weighted superposition of generalized stop operators. Its rheological representation [3], [4] combines a nonlinear elastic spring in parallel with multiple generalized vector-stop hysterons, and a thermodynamic interpretation of such nonlinear vector-stop elements for anisotropic materials is discussed in [5]. Thermodynamically formulated models provide physically meaningful constructions and allows decomposition of dissipative and reactive power [3].

In the thermodynamically formulated generalized Prandtl–Ishlinskii stop-type model [5], each hysteron is split into an elastic (reversible) part and a plastic (irreversible) part, where the plastic part constitutes the internal memory state. The hysteron output is a nonlinear function of the elastic part, yielding an implicit update equation for each hysteron and can only be resolved iteratively through the Newton method. This local, nonlinear, iterative update incurs considerable computational effort, particularly in finite element (FE) simulations.

Although \(\boldsymbol{\textit{V}}\) is local and parallelizable across integration points, the per-point hysteretic updates remain expensive. In this work, we propose a simplification of the generalized Prandtl–Ishlinskii stop operator within the thermodynamic framework, while retaining the essential hysteresis properties. This simplification linearly weights the stop-operator outputs, and for isotropic cases, this yields an analytical return-point mapping for the local hysteron updates, eliminating the need for iterative Newton procedures. Integrated into an FE framework, the simplified constitutive mapping enables faster evaluation of \(\boldsymbol{\textit{V}}[\mathbf{B}]\) at each integration point, reducing the overall computational cost.

The paper is organized as follows. In Sec. 2.1, we introduce the operator simplification and its rheological interpretation. The anhysteretic and hysteretic formulations are presented in Sec. 2.2 and Sec. 2.3, respectively. Thermodynamic consistency of the simplification is verified in Sec. 2.4, and parameter identification is briefly discussed in Sec. 2.5. In Sec. 3, we present a 2D FE simulation comparing computational performance and losses, and finally conclude in Sec. 4.

2 Simplified model and its formulation↩︎

2.1 Operator simplification↩︎

The generalized Prandtl–Ishlinskii stop-type operator defines the mapping between the output field intensity \(\mathbf{H}\) and input flux density \(\mathbf{B}\) in a discrete formulation [2], [5] as \[\mathbf{H} = \boldsymbol{\textit{V}}[\mathbf{B}] = w_1\,g(\mathbf{B}) + \sum_{i=2}^{N} w_i\,g(\boldsymbol{\textit{S}}_i[\mathbf{B}]), \label{eq:GeneralizedPrandtl}\tag{1}\] where \(g\) is a nonlinear memoryless mapping, \(\boldsymbol{\textit{S}}_i\) are the memory-dependent operators of \(i\)-th stop hysteron, \(w_1\) and \(w_i\) are the weight coefficients, and \(N\) denotes the total number of rheological elements (Fig. 1 (a)). The first term represents the anhysteretic response, while the summation accounts for the hysteretic contribution through a weighted superposition of \(g\) applied to the operators \(\boldsymbol{\textit{S}}_i\).

Within the thermodynamic framework, the memoryless mapping is identified as \(g(\cdot)=\nabla f(\cdot)\), where \(f\) is the Helmholtz free energy potential. Hence, each hysteron contribution takes the nonlinear form \(\nabla f\) applied to the stop operators \(\boldsymbol{\textit{S}}_i\). We introduce a simplification of 1 , retaining the anhysteretic term \(w_1\nabla f(\mathbf{B})\) while simplifying the hysteretic contribution by treating \(\nabla f\) as an identity mapping. The stop hysterons are thus weighted directly through their outputs \(\boldsymbol{\textit{S}}_i\) as follows \[\mathbf{H} = \boldsymbol{\textit{V}}[\mathbf{B}] = w_1 \nabla f(\mathbf{B}) + \sum_{i=2}^{N} w^\mathrm{new}_{i}\,\boldsymbol{\textit{S}}_i[\mathbf{B}]. \label{eq:SimplifiedPrandtlThermodynamical}\tag{2}\] Here, \(w^\mathrm{new}_{i}\) are refitted hysteron weights. In the rheological network, this is equivalent to replacing the nonlinear spring in each hysteron with a linear spring of unit stiffness (\(1\,\)N/m) as shown in Fig. 1 (b).

2.2 Anhysteretic characteristic↩︎

The spring-only branch is modeled by the isotropic free energy potential \(f\). It is expressed as a linear combination of \(m\) piecewise basis potentials \(p_{s_i}\in\mathbb{P}_2\) that depend on flux density \(\mathbf{B}\), with weights \(v_i\) as follows \[f(\mathbf{B}) = \sum_{i=1}^{m} v_i\,p_{s_i}. \label{eq:freeEnergy}\tag{3}\] Here, \(s_i\) denote the dead-zone thresholds [6]. The potential gradient \(\nabla f\) and the Hessian \(\nabla^2 f\) are further obtained as \[\nabla f(\mathbf{B}) = \sum_{i=1}^{m} v_i \nabla p_{s_i}\,, \qquad \nabla^2 f(\mathbf{B}) = \sum_{i=1}^{m} v_i \nabla^2 p_{s_i}. \label{eq:potentialGradHessian}\tag{4}\] As shown in Fig. 2, a zero dead-zone threshold \(s_i\) yields a linear basis potential gradient

while a nonzero threshold introduces a piecewise linear characteristic. The anhysteretic response \(\nabla f\) is thus constructed as a weighted superposition of such piecewise linear basis functions. To ensure smoothness of the gradient and Hessian for numerical robustness, a regularization is applied using the Dirac delta hat function \(\delta^{\kappa\varepsilon}\) (Fig. 3 (a)), where \(\kappa\) and \(\varepsilon\) are the regularization parameters. The coordinate shift \(z = \|\mathbf{B}\|-s_i\) is introduced to relocate the \(i\)-th threshold to the origin. Integrating \(\delta^{\kappa\varepsilon}\) twice yields \(C^2\) continuous basis potential gradients, as illustrated in Fig. 3 (b).

2.3 Hysteretic characteristic↩︎

The hysteretic response is produced by the stop hysteron (Fig. 4). For this section, a representative hysteron is considered and the hysteron index \(i\) in 1 , 2 is omitted for readability.

In the vector stop model, the flux density decomposes as \(\mathbf{B} = \mathbf{B}^e+\,\mathbf{B}^p\), where \(\mathbf{B}^e\), \(\mathbf{B}^p\) are the elastic and plastic parts, respectively, while \(\mathbf{H} = \mathbf{H}^e = \mathbf{H}^p\). The plastic element in both models (Fig. 4) is described by a strongly convex quadratic dissipation potential \(k\), dependent on the plastic field \(\mathbf{H}^p\) and scalar threshold \(r(\mathbf{B})\) as \[k(\mathbf{H}^p,r(\mathbf{B})) = \frac{1}{2} (\mathbf{H}^p)^\top \mathbf{H}^p - \frac{1}{2} r(\mathbf{B})^2 \leq 0. \label{eq:diss95pot}\tag{5}\] The threshold \(r\) depends on the flux density to deactivate hysterons under large \(\mathbf{B}\), where the response is nearly reversible. The constraint \(k\leq0\) defines a strongly convex admissible set \[K(r(\mathbf{B})) = \left\{\mathbf{H}^p\,|\, k(\mathbf{H}^p,r(\mathbf{B})) \leq 0\right\}, \label{eq:admissibleset}\tag{6}\] where \(k < 0\) inside the elastic domain and \(k=0\) on the boundary \(\partial K\) (Fig. 5). For isotropic cases, \(\partial K\) forms a circle in 2D or a sphere in 3D with radius determined by threshold \(r\). The maximum dissipation principle from plasticity theory [7] determines a unique plastic flow direction by finding the optimal hysteron field response \(\mathbf{H}^*\) that maximizes the dissipative power \(\mathbf{H}^\top\dot{\mathbf{B}}^p\) subject to \(k \leq 0\), i.e. \[\mathbf{H}^* = \arg\max_{\mathbf{H} \in K} \left(\mathbf{H}^\top\dot{\mathbf{B}}^p\right). \label{eq:MaximumDissipationPrinciple}\tag{7}\]

Solving 7 yields the evolution equation \[\dot{\mathbf{B}}^p = \lambda\nabla k, \label{eq:evolutionEquation}\tag{8}\] with Lagrange parameter \(\lambda\) and Karush-Kuhn-Tucker (KKT) conditions: \(k \leq 0\), \(\lambda \geq 0\) and \(\lambda k = 0\). For the generalized hysteron, the nonlinear elastic spring gives \(\mathbf{H} = \nabla f(\mathbf{B}^e)\) and 8 reduces to \(\dot{\mathbf{B}}^p = \lambda\mathbf{H} = \lambda\nabla f(\mathbf{B} - \mathbf{B}^p)\). Applying backward-Euler time discretization yields \[\mathbf{B}_{n+1}^p - \mathbf{B}_{n}^p = \eta_{n+1} \nabla f(\mathbf{B}_{n+1} - \mathbf{B}_{n+1}^p), \label{eq:GeneralizedDiscreteEvolutionEquation}\tag{9}\] where \(\eta\) is the discrete-time Lagrange parameter and \(n+1\) denotes the current time step.

For the simplified hysteron, a quadratic free energy potential with unit stiffness (\(c=1\) in Fig. 4 (b)) gives the hysteron response \(\mathbf{H} = \mathbf{B}^e\). As a result, 8 reduces to \(\dot{\mathbf{B}}^p = \lambda \mathbf{B}^e\). Further time discretization yields the update equation \[\mathbf{B}_{n+1}^p - \mathbf{B}_{n}^p = \eta_{n+1} (\mathbf{B}_{n+1} - \mathbf{B}_{n+1}^p). \label{eq:SimplifiedDiscreteEvolutionEquation}\tag{10}\] Both the updates 9 and 10 are solved using an elastic predictor-plastic corrector return point mapping scheme [8]. First, a trial state is computed based on previous plastic state \(\mathbf{B}^p_n\), assuming no plastic evolution: \[\begin{align} \mathbf{H}^\mathrm{trial}_{n+1} &= \nabla f(\mathbf{B}_{n+1}-\mathbf{B}^p_n), \\ \mathbf{H}^\mathrm{trial}_{n+1} &= \mathbf{B}_{n+1}-\mathbf{B}^p_n, \label{eq:HtrialForSimplified} \end{align}\tag{11}\] for the generalized and simplified models, respectively. If \(\|\mathbf{H}^\mathrm{trial}_{n+1}\| > r(\mathbf{B}_{n+1})\), a plastic correction projects the trial state back onto \(\partial K\) (Fig. 5). In the generalized model, this requires iterative Newton solves due to the nonlinear free energy potential. In the simplified model, the linear elastic relation yields the fully analytical return-point mapping \[\mathbf{H}_{n+1} = r(\mathbf{B}_{n+1}) \frac{\mathbf{H}^\mathrm{trial}_{n+1}}{\|\mathbf{H}^\mathrm{trial}_{n+1}\|}\,,\quad \eta_{n+1} = \frac{\|\mathbf{H}^\mathrm{trial}_{n+1}\|}{r(\mathbf{B}_{n+1})} -1\,, \label{eq:returnPointMappingSimplified}\tag{12}\] enabling direct analytical hysteron updates without iteration.

2.4 Thermodynamic consistency↩︎

The energy balance for the \(i\)-th hysteron decomposes the input power into stored and dissipated contributions as \[\mathbf{H}_i^\top\dot{\mathbf{B}}_i = \mathbf{H}_i^\top\dot{\mathbf{B}}_i^e + \mathbf{H}_i^\top\dot{\mathbf{B}}_i^p,\] where \(\dot{\mathbf{B}}_i^e\) and \(\dot{\mathbf{B}}_i^p\) are the elastic and plastic flux rates, respectively. The hysteresis loss density \(p[\mathbf{B}]\) is obtained as a weighted sum of dissipative contributions from all hysterons in Fig. 1 (b) as, \[p[\mathbf{B}] = \mathbf{H}^\top\dot{\mathbf{B}}^p = \sum_{i=2}^{N} w^\mathrm{new}_{i}\,\mathbf{H}_i^\top\dot{\mathbf{B}}_i^p. \label{eq:dissipativeTerm}\tag{13}\] Substituting \(\mathbf{H}_i = \mathbf{B}^e_i\) from the simplified model and the evolution equation \(\dot{\mathbf{B}}_i^p = \lambda \mathbf{B}_i^e\) into 13 yields \[\mathbf{H}^\top\dot{\mathbf{B}}^p = \sum_{i=2}^{N} w^\mathrm{new}_{i}\,{\mathbf{B}_i^e}^\top \lambda \mathbf{B}_i^e, \label{eq:finalDissipativeExpression}\tag{14}\] which involves two key quantities: \(\lambda \geq 0\) from KKT conditions and \({\mathbf{B}_i^e}^\top\mathbf{B}_i^e \geq 0\). The non-negativity of weights \(w^\mathrm{new}_{i}\) is enforced during parameter identification, ensuring overall non-negative dissipation and thus thermodynamic consistency of the simplified formulation.

2.5 Parameter identification↩︎

Parameter identification proceeds in two stages. The first stage identifies the anhysteretic parameters through a nested nonlinear program for the dead-zone thresholds \(s_i\) and an inner quadratic program for the coefficients \(v_i\), both defined in 3 . To ensure strongly convex potential \(f\) and monotonous gradient, certain constraints are imposed on \(v_i\) that can be found in [6]. Deep saturation is further captured by extrapolating the potential gradient to the inverse of the permeability of free space. In the second stage, the hysteron weights \(w^\mathrm{new}_{i}\) are obtained by minimizing the error between the overall model response and the measured field value \(\mathbf{H}_\mathrm{meas}\), \[\mathbf{e} = \sum_{i=1}^{N}w^\mathrm{new}_{i}\,\mathbf{H}_i - \mathbf{H}_\mathrm{meas}, \label{eq:errorHysteronParameters}\tag{15}\] where \(\mathbf{H}_i\) is the field response from \(i\)-th branch in Fig. 1 (b).

3 Numerical results↩︎

3.1 Comparison of \(B\)–\(H\) curves and hysteresis loss density↩︎

The material response is first assessed at a single point, driven by a 2D oscillating flux of increasing amplitude (Fig. 6) in \(B_x\) and \(B_y\) without field coupling. Both models produce nearly identical \(B_x\)–\(H_x\) responses, capturing similar minor loop behavior across the excitation range (Fig. 7).

The hysteresis loss density comparison (Fig. 8) yields an overall relative difference of approximately \(4.71\%\), confirming comparable accuracy of the simplified model.

3.2 FE Simulation of 2D electric machine↩︎

Both models are then applied to a 2D Permanent Magnet Synchronous Machine (PMSM) FE model (Fig. 9), discretized with 10,000 elements. The stator and rotor are assigned an isotropic hysteresis material (M330-35A). The stator coils are excited with 3-phase sinusoidal currents and a static simulation was performed for 90 rotor positions.

Two simulation cases are considered: case (i) evaluates hysteresis and excess losses without classical eddy-current coupling, and case (ii) additionally includes eddy-current losses at the lamination scale. The generalized model is evaluated using two approaches. The online approach computes potential, gradient, and Hessian values during simulation whereas the offline approach retrieves precomputed values via interpolation, resulting in a negligible maximum interpolation error (\(\approx\) \(0.02\%\)). The simplified model is evaluated with the online approach only. In case (i), the simplified model achieves \(88\%\) and \(61\%\) runtime reductions compared to the generalized model using online and offline approaches, respectively. The maximum relative difference (relative to generalized model with offline approach) in hysteresis loss is \(\approx\) \(6.8\%\) (Fig. 10),

while the maximum excess loss difference remains negligible at \(\approx\) \(0.12\%\), since the same excess loss parameters are used for both models. Furthermore, the number of Newton iterations required to solve the nonlinear system of equations was higher at a few rotor positions, reaching a maximum of \(8\) compared to \(4\) for the generalized model. However, each iteration is computational cheaper due to the simplified material response evaluation, and convergence is achieved with a consistent decrease in the residual. In case (ii), the maximum hysteresis loss difference reduces to approximately \(0.3\%\) (Fig. 11), as eddy current coupling dominates the field solution and drives both models toward a similar field distribution, resulting in close agreement in the loss values. The maximum relative difference for excess and eddy current losses is \(\approx\) \(0.16\%\) and \(\approx\) \(0.18\%\), respectively.

4 Conclusion and Outlook↩︎

In this work, a thermodynamically consistent simplification of the generalized Prandtl–Ishlinskii stop operator is introduced by retaining the nonlinear anhysteretic response and simplifying the hysteretic part. For isotropic case, this enables a fully analytical return-point mapping, replacing iterative Newton updates with a closed-form solution. Applied to a 2D PMSM FE simulation, significant runtime reductions of \(88\%\) and \(61\%\) are achieved compared to the generalized model with online and offline approaches, respectively, with comparable loss accuracy. In future work, we will extend the simplification to anisotropic cases and investigate similar analytical updates.

References↩︎