Flexible identification procedure for thermodynamic constitutive models for magnetostrictive materials

We present a novel approach for identifying a multiaxial thermodynamic magneto-mechanical constitutive law by direct bi- or trivariate spline interpolation from available magnetization and magnetostriction data. Reference data are first produced with a multiscale model in the case of a magnetic field and uniaxial and shear stresses. The thermodynamic model fits well to the results of the multiscale model, after which the models are compared under complex multiaxial loadings. A surprisingly good agreement between the two models is found, but some differences in the magnetostrictive behaviour are also pointed out. Finally, the model is fitted to measurement results from an electrical steel sheet. The spline-based constitutive law overcomes several drawbacks of analytical approaches used earlier. The presented models and measurement results are openly available.


Introduction
Modelling tools for coupled magneto-mechanical effects in magnetic materials are needed for analysing losses and vibrations in electrical machines, as well as magnetostrictive devices and actuators. The problem is rather complex due to its three-dimensional and multiaxial nature, since many different combinations and orientations of magnetic fields and mechanical stresses can occur in real applications. In electrical machines, mechanical stresses typically cause adverse effects by increasing the losses [1][2][3]. On the other hand, the same effects allow converting mechanical energy into electrical energy for energy-harvesting purposes [4][5][6]. The coupled magneto-mechanical effects have been studied experimentally [7][8][9][10], while the modelling in the field has evolved from purely uniaxial studies [11] to equivalent stress or strain approaches for multiaxial loadings [12][13][14][15], and further towards genuinely multiaxial models [16][17][18][19].
Thermodynamic approaches provide one possible way of deriving coupled constitutive laws. Couplings between magnetic, electric and thermal fields [20][21][22] as well as elastic and plastic deformations [23][24][25][26] have been considered through thermodynamic frameworks. The applications have included both isotropic [20,22] and anisotropic [26] solids as well as granular materials [27][28][29][30][31]. In [32], an extensive review of earlier magneto-elastic constitutive models is given, and implicit constitutive equations based on 21 scalar invariants are derived for bodies under large deformations. These equations are then simplified for small fields and applied for solving magneto-elastic boundary value problems for slabs and an annulus. However, systematic procedures for identification of such models based on experimental data have not been reported in details so far. It is particularly mentioned in [32] that the given equations are 'too general to be used to correlate with experimental data as there are so many material functions that depend on numerous invariants.' In this paper, we focus on identification of thermodynamic constitutive models based on the theory of invariants.
In our earlier works, we have developed a thermodynamic approach for deriving coupled magneto-mechanical constitutive equations for ferromagnetic materials [33]. The approach has been based on expressing a Helmholtz free energy per unit volume ψ(B, ε) as a function of the magnetic flux density vector B and total small strain tensor ε, i.e. the symmetric part of the displacement gradient. The magnetization and stress have then been derived by the Coleman-Noll procedure from the dissipation inequality, and are given by ∂ψ ∂B T and σ = ∂ψ ∂ε . (1.1) The field strength then becomes where μ 0 is the permeability of free space. By assuming an isotropic material, the dependency of the Helmholtz free energy density on B and ε reduces to the following six invariants I 1 = tr ε, I 2 = tr ε 2 , I 3 = tr ε 3 , (1.3) where e = ε − 1 3 (tr ε)I is the deviatoric strain, I being the second-order identity tensor and tr denotes the trace of a tensor. Equation (1.1) thus becomes So far, we have mainly used analytical expressions for ψ(I 1 , I 2 , I 3 , I 4 , I 5 , I 6 ). The exact form of the function has varied (see different variations in [33][34][35][36][37][38]), but has been more or less similar to (1.5) in which λ and μ are the Lamé parameters yielding linear Hooke's Law (appendix A), and the nonlinear magneto-elastic part is described by the last three polynomials with α i , β i and γ i as fitting parameters. Invariant I 3 has usually been neglected in order to reproduce linear behaviour under purely mechanical loading. Although these analytical expressions are relatively easy to implement, they have sometimes failed to properly fit to measured magnetization and magnetostriction curves. More complicated analytical expressions for the magneto-elastic part might improve the fitting, but finding a suitable expression is rather time-consuming and does not bring any added physical insight into the phenomena. In addition, if the state variables are to be changed, a different analytical expression might need to be used. Analytical expressions for the energy functions were also used in [21] and [22]. In this paper, we propose a new approach for overcoming the disadvantages of analytical energy density expressions. The idea is to express the energy density as a spline of the invariants, identifying the free energy density function as a direct least-squares spline fit to measured magnetization and magnetostriction curves. The main advantage of the proposed approach compared to [33][34][35][36][37][38] is that no a priori assumptions on how to analytically express the free energy density need to be made. This provides more flexibility when fitting the model against measurements, and also allows keeping the model fitting procedure unchanged, even if the input variables need to be changed. Using B as the magnetic input variable is convenient if the constitutive law is to be used in combination with a finite-element (FE) formulation based on the magnetic vector potential, while H is more suitable for magnetic scalar potential formulations. Using ε as the mechanical input variable is convenient in displacement-based FE formulations, while using σ is often straightforward in simple statically determined structures.
In the following, we first describe our measurement set-up and demonstrate the problems related to fitting the analytical energy density expression (1.5) against measurements from an electrical steel sheet. The new spline-based approach overcoming the problems of the analytical approach is then described. The splines are fitted to simulation results produced by the simplified multiscale (SMS) model of [17], and the behaviour of the models under complex multiaxial loadings is compared. Finally, it is demonstrated that the proposed approach provides a significantly better fit against measurements than the analytical expression (1.5).

Methods
(a) Measurement set-up and problems with the analytical models A custom-designed stressing mechanism was built offering the possibility of in-plane uniaxial stressing of a modified single sheet tester (SST) sample. The mechanism was designed to exert both compressive and tensile stress up to 100 MPa for a 25 × 0.5 mm cross-sectioned SST sample. A manual two-way locking screw system in series with a helical compression spring was used to exert the required force. Use of spring was especially considered to better control the stress levels, thus a force resolution of 1 N was achieved with this device. In addition to the stressing mechanism, designing the SST sample is essential as well. In order to achieve a uniform stress at the measurement zone, i.e. centre of the SST sample, proper propagation of force from the clamped zone to the centre of the sample is required. Thus the shape of the SST sample is chosen according to tensile testing standards.
For the magnetic measurements, tunnelling magneto-resistance (TMR) sensors of dimension 6 × 5 × 1.5 mm arranged in a 2 × 2 grid on a printed circuit board of 20 × 20 mm were used to measure the surface magnetic field strength. When placed on the SST sample, the sensing element  The stressing mechanism applies a static force F into the sample. The SST is excited by a voltage waveform U(t) and the voltage waveforms from the field strength sensor (U H (t)), strain gauge (U λ (t)) and search coil (U B (t)) are measured. An input voltage waveform U ref (t) is iteratively searched such that sinusoidal U B (t) is obtained. (Online version in colour.) of the TMR sensor chips was effectively at 0.4-0.5 mm above the sample surface. Similarly, a search coil wound around the SST sample was used to measure the average magnetic flux density. At the low measurement frequency of 6 Hz, the demagnetization field by the eddy currents is small and the surface field strength and average flux density give a sufficient approximation of the static constitutive behaviour. Furthermore, a non-inductive three element strain gauge rosette of 60°delta arrangement (H-series rosette from Micro-Measurements), glued on the surface of the sample (with the insulation coating removed) was used to measure the magnetostriction. The individual element of the rosette has the gauge length and resistance of 3.18 mm and 700 Ω, respectively. Figure 1 shows a schematic diagram of the measurement set-up.
A programmable power source (AMETEK CSW 5550VA) and a data acquisition system (DAQ-NI USB-6251 BNC) with analogue output were used in conjunction with a PC to control the magnitude and the waveform of the supply voltage so as to produce a sinusoidal induction in the SST sample. The feedback control of the supply voltage was programmed using the Matlab/DAQ toolbox. Low-noise/high-gain signal amplifiers were also used to amplify the signal obtained from the search coils. In addition to that, a high sampling rate DAQ system (DEWETRON) controlled by PC/DEWEsoft was used to retrieve the measured signals for the field strength and the flux density. The sample was magnetized using two vertical cores and the excitation coils with a total of 1000 turns. Within the limitation of the power amplifier's operation, the total turns of the excitation coils were sufficient for a wide range of induction amplitudes and frequencies. Figure 2 shows measurement results for magnetization curves and magnetostriction curves under different compressive (−) and tensile (+) stresses in M400-50A electrical steel. The single-valued curves have been obtained by averaging the hysteresis loops in order to identify anhysteretic material models. Although ferromagnetic hysteresis is a major source of losses in electromagnetic devices, anhysteretic material properties are usually sufficient to estimate the magnetic field distribution in devices with air gaps, and hysteresis losses can be calculated a posteriori with a reasonable accuracy [39,40]. Thus, only anhysteretic constitutive laws are considered in this paper.  The problem related to the analytical energy density expressions is demonstrated by setting n α = 5, n β = 1 and n γ = 1 in (1.5) and fitting α i , β i and γ i against the measurement results. In the fitting, the magnetization and magnetostriction curves derived from (1.5) at four different stress values were compared against the measured curves. The curves were normalized by their maximum values to give equal weights to the magnetization and magnetostriction curves. Since the stress is uniaxial in the measurements, the components of ε were first iterated separately for each value of B x such that the desired stress was obtained. This makes the problem nonlinear with respect to the parameters and requires using a nonlinear minimization algorithm. The fitting was attempted using both nonlinear least-squares minimization and a genetic algorithm (Matlab functions lsqcurvefit and ga). Figure 3 shows the best results obtained from the fitting. Although the magnetization curves correspond relatively well to the measurements, the magnetostriction curves do not. This is caused by the limitations in the analytical expression (1.5), which fails to correctly describe the physical behaviour of the material. Derivation of a better expression would require a time-consuming trial and error process and would probably lead to complex analytical formulae. In addition, finding suitable initial values for the nonlinear fitting problem may be difficult. In the following, we propose a new identification approach based on linear least-squares spline fitting which does not require analytically expressing the energy density. From now on, we focus mainly on the magneto-elastic part of the free energy density ψ discussed in §1, and denote this by φ. In addition, ψ and φ are referred to simply as energy densities (instead of Helmholtz energy density), since their physical meaning varies depending on the choice of state variables. We start by deriving the thermodynamic model by choosing the magnetic field strength vector H and the stress tensor σ as the state variables. This choice is preferred, since it will allow straightforward comparison to the SMS model of [17], which uses the same state variables. In this case, the relevant invariants become I 4 = H · H, I 5 = H · sH and I 6 = H · s 2 H, where s = σ − 1 3 (tr σ )I is the deviatoric stress. The magnetic flux density B and magnetostriction tensor λ can then be obtained as It is noted that if the magnetic polarization μ 0 M is preferred as the output variable instead of B, the energy can be changed to φ → φ − 1 2 μ 0 I 4 without affecting the magnetostriction. However, we prefer to write the constitutive models directly in terms of B and H since these are the variables of interest when numerically solving Maxwell equations.
The measurement set-up described in §2a allows characterizing the magneto-mechanical constitutive behaviour in the form of B(H, σ ) and λ(H, σ ), where σ is uniaxial and parallel to B and H (figure 2). Devices for biaxial stressing [13] and shear stressing [10] have also been presented. In the following derivations, it will become apparent that consideration of uniaxial and shear stresses is sufficient to identify the complete energy density. Let us thus assume for

so that the invariants become
Since the material is assumed to be isotropic, only stress-induced magnetic anisotropy is present, and B y can be non-zero only if σ xy = 0. Under uniaxial stress, σ xy = 0, and only two of the invariants are independent. This means that we can identify the energy density only as a bivariate function, which is denoted φ (2) (I 4 , I 5 ). However, if measurement data were available also for σ xy = 0, we could identify trivariate φ (3) (I 4 , I 5 , I 6 ). It is emphasized, that with the given magnetic field in the x-direction, I 6 is independent of I 4 and I 5 only if at least one of the shear stress components σ xy or σ zx is non-zero.
In this paper, we aim for a general identification approach, which does not require making assumptions on the form of functions φ (2) (I 4 , I 5 ) or φ (3) (I 4 , I 5 , I 6 ). Thus, instead of analytical expressions for the functions, we express them as bi-and trivariate B-splines of the input variables. The splines are then fitted directly on the measured B and λ. However, the difficulty is that the values of H x , σ xx and σ xy , on which B and λ are measured, are often chosen so that they lie on a regular (but not necessarily uniform) grid. This means that the invariants I 4 , I 5 and I 6 will be irregularly distributed, which makes it difficult to construct the splines. To avoid this problem, we define auxiliary variables which under the aforementioned H and σ become u = H x , v = σ xx and w = σ xy , which are regularly distributed. In addition, with these variables, we obtain We can thus identify the splines φ (2) (u, v) and φ (3) (u, v, w) directly from the measured flux density and magnetostriction. Note the factor 2 in the last term of (2.4), which results from the fact that the symmetry assumption λ yx = λ xy is done before differentiating the energy. It is emphasized, that after the splines are identified, the auxiliary variables (2.3) and the splines are well defined for arbitrary three-dimensional field and stress configurations.
To formulate the spline fitting problem, let us assume that the regular grid formed by H x , σ xx and σ xy consists of N u , N v and N w values of each variable, respectively. The energy densities are expressed as bi-and trivariate B-splines of u, v and w. The spline coefficients are denoted c lm and c lmn for the bi-and trivariate cases, respectively, where l = 1, . . . ,M u , m = 1, . . . ,M v and n = 1, . . . ,M w . To avoid overfitting, we typically want to have at most as many coefficients as measurement points: The bi-and trivariant energy densities in the grid points thus become where i = 1, . . . ,N u , j = 1, . . . ,N v and k = 1, . . . ,N w denote the indices of the measurement points, and the A-terms denote elements of collocation matrices A u , A v and A w which map the Bspline coefficients into the energy densities in the measurement points, and which depend only on the chosen grid and the order of the B-spline. In the notation, the repeated subindices are summed over. If the energy densities in the measurement grid were measured directly, the spline coefficients could be directly solved from (2.5). However, since only the partial derivatives (2.4) are measured, we end up with where the D-terms are elements of collocation matrices D u , D v and D w , which map the spline coefficients into the measured partial derivatives. The problem is overdetermined. In the bivariate case, we have 2N u N v equations for M u M v variables, and in the trivariate case 3N u N v N w equations for M u M v M w variables. Details of the solution will be discussed in §2d. The order of the B-splines can be chosen freely, but it should be greater or equal to 3 to make the second derivatives continuous, which is beneficial when solving nonlinear field problems using Newton-Raphson iteration. Outside of the range of measured data, the splines are extrapolated as quadratic functions, meaning linear extrapolation for the partial derivatives. This is a reasonable assumption if identification data is available close up to saturation, since both the B(H) and ε(σ ) relationships tend to become linear above saturation.
(c) Fitting of the thermodynamic model in terms of B and ε As mentioned in the previous section, measurements are typically available under known stress. Formulation of the model using the strain as the state variable thus needs some more considerations. Let us assume now that instead of H and σ , the state variables are B and ε, where ε = C -1 σ + λ is the total strain (C being the mechanical material stiffness matrix related to Hooke's Law and defined by Young's modulus E and Poisson's ratio ν). This choice of state variables is convenient, if the constitutive law is to be implemented in an FE formulation using the magnetic vector potential and mechanical displacement. The relevant invariants now become I 1 = tr ε, I 2 = tr ε 2 , I 4 = B·B, I 5 = B · eB and I 6 = B · e 2 B, where e = ε − 1 3 (tr ε)I. The total energy density function becomes ψ(I 1 , I 2 , I 4 , I 5 ,  H(B, σ )  and λ(B, σ ) where e xx is the xx-component of the deviatoric strain. Magnetostriction is assumed to be isochoric so that The diagonal strain components are 11) and the xx-component of the deviatoric strain is thus In addition, the shear component of the strain can be written as 13) and the invariants I 5 and I 6 thus become Using the same definitions as in (2.3) for the auxiliary variables, they become The notation v and w is used to denote that these variables are not yet regularly distributed due to the small magnetostriction terms λ xx and λ xy , which depend on B x , making v and w dependent on u. However, for each given value of u, we can easily interpolate new values H x , σ xx , σ xy , λ xx and λ xy for the field, stress and magnetostriction components using simply onedimensional interpolation so that v and w become regularly distributed: With these variables, we finally obtain (2.18)

(d) Implementation considerations
The fitting of the splines was implemented in Matlab 1 (v.R2017b). For the given measurement grid u, v and w and the desired spline order, the B-spline knots were calculated with the Matlab function optknt. Given the grid, the order and the knots, the collocation matrices A u , A v and A w in (2.5) were obtained with the function spcol. In order to construct the collocation matrices D u , D v and D w in (2.6), spcol was first used to obtain collocation matrices for mapping the partial derivative B-spline coefficients into the measured partial derivative values. These were multiplied with matrices obtained by function DerivBKnotDeriv (found in [41]), which map the B-spline coefficients into partial derivative B-spline coefficients. Different scaling factors should be chosen for different equations in the systems (2.6) and (2.18) in order to equalize their weighting and significance during the least-squares solution. For example, if in the bivariate case in (2.6), B x,ij are expressed in teslas and λ xx , ij in m m −1 , the magnetostriction has negligible effect on the least-squares solution, since its numerical values are in the range of 10 −6 . In our implementation, instead of solving (11), we divided the equations by the maximum corresponding measurement values: and thus equalized the weighting for the flux density and magnetostriction. The full equation system can be written as Ac = b, where the vector c includes the spline coefficients, A is the system matrix and vector b includes the scaled measurements. The sizes of A, in the trivariate case. However, since only the partial derivatives of the energy are measured, this system does not have a unique solution. A constant could be added to the energy without changing the partial derivatives. Thus, one of the spline coefficients has to be fixed to a constant to obtain a unique solution. In our implementation, this is done by writing the least-squares systems only for c lm with l + m > 2 and c lmn with l + m + n > 3 in the bi-and trivariate cases, respectively, fixing manually c 11  The ranks and condition numbers of the system matrices before (r(A) and κ(A)) and after (r(A ) and κ(A )) removing the columns corresponding to c 11 and c 111 will be compared in the results section to give insight on the uniqueness of the solution.

(e) Simplified multiscale model
It is important to verify if the thermodynamic constitutive law can be identified from measurements under H, σ xx and (possibly) σ xy in such a way that it can reasonably reproduce the magneto-elastic behaviour under arbitrary multiaxial loadings. Such a verification would require an extensive amount of measurements in complex loading conditions in three dimensions, which are not currently possible. Thus, the simplified three-dimensional multiscale material model developed in [17] is used as a reference. The SMS model is applied to produce simulation data for both identification and validation of the thermodynamic constitutive law. The SMS model is based on expressing a local potential energy for a domain α oriented along a given direction u α as a function of the field strength H and stress σ as are the magnetization and magnetostriction of the given domain, determined by the saturation values M s and λ s , respectively. The volume fraction of domains is expressed as where χ 0 is the initial magnetic susceptibility of the material. The magnetization vector and magnetostriction tensor are then obtained by weighting the domain magnetizations and magnetostrictions by the volume fraction and integrating over all possible orientations as In our implementation, the integrations over two direction angles defining a unit spherical shell are carried out numerically using adaptive two-dimensional quadratures. The magnetic flux density is obtained as B = μ 0 (H + M).

Application and results
where vector x contains the multiscale model results in the N u × N v = 440 or N u × N v × N w = 4840 measurement points, and vector x fit contains the corresponding results obtained by differentiating the fitted spline. Figure 4 shows the fitting results for the bivariate spline φ (2) (H, σ ). In this case, the fitting errors r fit for B x and λ xx were 0.0037% and 0.0059%, respectively, which means that a single energy function can represent both the flux density and magnetostriction consistently to the SMS model under uniaxial loading. The ranks and condition numbers of the system matrix are r(A) = 339 and κ(A) = 8.3 × 10 17 before and r(A ) = 339 and κ(A ) = 3.1×10 4 after removing the column corresponding to coefficient c 11 which is fixed manually to zero. The unchanged rank and the large drop in the condition number show that the system becomes well-conditioned after removing one column, making the solution unique. Figure 5 shows the fitting results in the trivariate case φ (3) (H, σ ) for three different values of the field strength H x . In this case, the errors for B x , λ xx and λ xy were 0.0041%, 0.0059% and 0.0017%, respectively. The ranks and condition numbers are r(A) = 4839 and κ(A) = 6.1×10 17 and r(A ) = 4839 and κ(A ) = 2.7 × 10 5 , again showing the uniqueness of the solution after fixing c 111 . Noteworthy is that the error does not increase from the bivariate case, which means that also the behaviour under shear stress is similarly predicted by the thermodynamic and SMS models. The agreement between the two models is surprisingly good.
From figure 4b, it can be seen that the SMS model produces magnetostriction under stress even when the magnetic field is zero. This behaviour is explained as the so-called E effect [42], which can be observed in ferromagnetic materials, and is often described as a change in the Young modulus E under stress. The spline-based thermodynamic model is able to reproduce this behaviour contrary to the analytical energy density functions (1.5) we have used earlier.
The same SMS results were next interpolated into the form H x (B x , σ xx , σ xy ), λ xx (B x , σ xx , σ xy ) and λ xy (B x , σ xx , σ xy ), and the bi-and trivariate splines were fitted in terms of B and ε, as described in §2c, again with an order of 3. The values used for the Young modulus and the Poisson ratio were E = 183 GPa and ν = 0.34. Figure 6 shows the fitting results for the bivariate spline φ (2) (B, ε). In this case, the fitting errors for H x and τ xx were 0.019% and 0.016%, respectively. Figure 7 shows the fitting results in the trivariate case φ (3) (B, ε) for three different values of the flux density B x . In this case, the errors for H x , τ xx and τ xy were 0.020%, 0.017% and 0.009%, respectively. The ranks and condition numbers behave similarly to the previous case, ensuring the uniqueness of the solutions. The errors grow slightly from the H-and σ -based model, mainly due to the required  numerical interpolations. Nevertheless, the spline-based approach is clearly able to cope with different combinations of input variables.

(b) Sensitivity analysis
It is important to study the robustness of the fitting procedure against measurement errors and noise, which are inevitable in real experiments. When measurements of magnetization and magnetostriction curves are obtained, the data are usually smoothed and denoised in a preprocessing stage before fitting constitutive models. Such preprocessing is likely to alter the relationship between the magnetization and magnetostriction curves in such a way that they cannot be exactly described by a single thermodynamic potential. To study the sensitivity of the   fitting procedure against such errors, artificial error is introduced to the SMS model results by changing the output flux density and magnetostriction as where the error r varies between 0.001% and 100%. λ xy is kept unchanged in the trivariate case. The fitting of the bivariate and trivariate splines φ (2) and φ (3) in terms of H and σ is repeated similarly to the previous section, and the fitting errors r fit from (3.1) are studied as a function of r.
The results are shown in figure 8. It is seen that the fitting errors r fit increase slowly with respect to r and remain below 10% for each case when r < 100%. The fitting procedure is thus robust and able to find a reasonable thermodynamic potential even for the case where the input data do not exactly fulfil the thermodynamic principles. These stresses have principal axes parallel to the x-, y-and z-axes, and thus the magnetic flux densities simulated with the two models become parallel to the field strength, which allows comparing the models in terms of the relative permeability. Figure 9a,b compares the relative permeability and magnetostriction λ xx (component parallel to the field) simulated with the SMS model and the bivariate spline. The results with the trivariate spline are almost equal to the bivariate ones and thus not separately shown. Both models behave very similarly, and the extrapolation of the spline appears to work sufficiently. However, when the magnetostrictions λ yy and λ zz (components perpendicular to the field) are compared in figure 9c,d in the cases of equibiaxial and pure shear loadings, rather significant differences are observed. The spline-based thermodynamic model produces exactly the same magnetostrictions in the yy-and zz-directions, despite the fact that the stress affects only in the xy-plane. This is a property of the invariant model and directly results from the fact that under a uniaxial field H = [H x 0 0] T and an arbitrary stress tensor σ , the invariants I 5 and I 6 become leading to ∂I 5 ∂σ yy = ∂I 5 ∂σ zz and ∂I 6 ∂σ yy = ∂I 6 ∂σ zz , (3.4) and thus λ yy = λ zz . The effect of the stress is thus transversely isotropic in the plane perpendicular to the field. This is not the case with the SMS model, in which the effect of the stress on the magnetostriction is anisotropic even in the plane perpendicular to the field. Threedimensional measurements under multiaxial loadings would be needed in order to validate the model predictions. Both models produce volume-preserving magnetostriction, such that λ xx + λ yy + λ zz = 0.  which can be rotated in the xy-plane by varying angle θ. The magnitude of the stress is σ xy = 50 MPa, and θ is varied from 0 to 180°. Figure 10a compares the magnetic flux density components B x and B y produced by the SMS and thermodynamic models at different angles θ. As expected, the bivariate spline model fails to accurately present the behaviour under shear stress, since the shear is not accounted for during the identification process. On the contrary, the flux density produced by the trivariate model agrees well with the SMS model. The magnetostrictions λ xx and λ yy are compared in figure 10b. Similarly to the results of figure 9c, λ yy differs clearly between the two models, and the bi-and trivariate models agree well. Finally, the magnetostrictions λ xy and λ zz are compared in figure 10c. The bivariate model produces zero λ xy , while the trivariate model agrees well with the SMS model. λ zz differs again similarly to figure 9d.

(e) Fitting to measurement results
The proposed thermodynamic approach was also tested against the uniaxial measurements discussed in §2a and figure 2. Following the idea described in §2c, the model was fitted in terms of B and ε as a bivariate spline φ (3) (B, ε). N u = M u = 100 and N v = M v = 17 values were used for B and variable v, and the spline coefficients, respectively. The results are shown in figure 11. Both the magnetization curves and the magnetostriction are satisfactorily fitted, the fitting errors for B x and τ xx being 2.5% and 1.9%, respectively. However, when comparing to the analytical   figure 3, the fitting is significantly better. Since the spline-based model does not make any assumptions on the shape of the energy density, the only factor affecting the quality of the fit is the correctness of the assumption, that a single thermodynamic potential is able to describe the magnetization and magnetostriction simultaneously.

Discussion and conclusion
An isotropic spline-based thermodynamic approach was presented for modelling coupled magneto-mechanical behaviour in ferromagnetic materials. The model is based on defining an energy-density functional, which depends on the invariant basis defined by the input vector and tensor, but assumptions on the exact form of the energy density are not made. The model is thus flexible and allows choosing freely the state variables, as a function of which the energy density is expressed. Flux density and strain are convenient choices, when the model is applied in finite-element analysis using formulations with magnetic vector potential and mechanical displacement. Magnetic field can be used in magnetic scalar potential formulations or if the model is to be coupled to a Jiles-Atherton hysteresis model, as in [36]. If the stress is known during the measurements, using it as a state variable simplifies the model derivations and the identification. The proposed identification procedure is not limited to magneto-elastic behaviour but can most likely be used also for other coupled problems, in which the constitutive laws are based on defining a free energy density and suitable measurements are available.
Comparison to the SMS model in the bi-and trivariate cases showed that both the multiscale model and the thermodynamic approach produce very similar magnetization curves. This has not been observed earlier when using analytical energy-density expressions, since these have prevented accurately fitting the thermodynamic model. The excellent agreement in the magnetization curves and magnetostrictions parallel to the magnetic field is rather surprising, considering that the two models are based on different modelling considerations. The observed differences in the magnetostrictions perpendicular to the field should be further studied and validated by measurements.
The main advantage of the multiscale models is that they rely on few intrinsic parameters, which have clear physical meaning and which are rather easy to obtain from measurements. However, deriving the models with input variables other than the magnetic field and stress is difficult, and it is thus complicated to implement the models in numerical tools. Since the proposed thermodynamic approach can be rather easily fitted against the multiscale model results, it could be used as a tool for simplifying the numerical implementation while accurately preserving the constitutive behaviour obtained from the multiscale model.
When fitting the spline-based model against measurement results, measurement error and noise obviously affect the quality of the fit. Experiments have shown that downsampling the measurement results and increasing the order of the B-spline make the fitted splines less prone to oscillations. In addition, the fitting can be improved by increasing the weight of, for example, the magnetization curve with respect to the magnetostriction curve before the least-squares solution.
If an isotropic material can be assumed, the three invariants I 4 , I 5 and I 6 are enough to describe the behaviour under arbitrary three-dimensional loadings based on the measurements in the regular grid described by the auxiliary variables u, v, w. Orthotropic anisotropy may be accounted for by considering two additional direction vectors describing two perpendicular symmetry planes and extending the scalar integrity basis of the free energy density to account for these vectors. This will result in 19 scalar invariants. In theory, a similar spline-based identification approach would be possible also in the anisotropic case, but a huge amount of measurement data would be required to identify the splines. The model could be further extended to plastic deformation by considering the equivalent plastic strain as an additional invariant, but identification of such a model would also require new measurements. To simplify the identification process, these measurements could also be replaced by multiscale models, which can account for anisotropy [43] and plastic deformation [44].
In this work, the Maxwell stress was neglected from both the thermodynamic approach and the SMS approach. The Maxwell stress could be considered in the momentum balance equations when solving coupled magneto-mechanical field problems. The Maxwell stress is obviously included in the measurements against which the thermodynamic model was fitted. However With a similar reasoning but taking into account that each shear component contributes twice to the energy due to symmetry of the stress and strain tensors, we get