Essential Measurements for Finite Element Simulations of Magnetostrictive Materials

We discuss which magnetoelastic material properties are essential to measure in order to model magnetostrictive materials in finite element simulations. We show knowing the magnetic constitutive relation is sufficient, if the elastic behavior without magnetic field is known a priori. We neglect hysteresis, and our starting point is to express the effect of mechanical deformation on the magnetic constitutive relation with a small strain tensor and magnetic flux density. It follows that the (energetic) state of a magnetostrictive material is independent of its history. Then, a certain choice of history allows us to keep magnetism and elasticity distinct. We demonstrate with open source software Elmer, how one can set up such magnetoelastic simulations. These simulations rely on data obtained from magnetostrictive measurements. Finally, it is discussed how a measurement setup and the finite element model should be combined in order to verify the approach with experiments.


I. INTRODUCTION
M AGNETOSTRICTION was discovered by J. Joule and its applications [1] include rotating motors [2], [3], magnetic actuators [4], [5], and energy-harvesting [6]. The magnetostrictive materials offer many potential advantages. For instance, in the case of electric linear or rotating motors, magnetostrictive materials allow precise control and comparatively high torque at low speeds [7].
The purpose of this paper is to present a model for finite element simulations of magnetostrictive systems. The main focus is on the properties of magnetostrictive materials. In particular, we consider which material properties are essential to measure and how they can be exploited in the model.
There are two main approaches to finite element modeling of magnetostrictive materials. These are the force-based approach and the strain-based approach [8]. This paper relies on the force-based approach, where the forces are obtained by the virtual work principle.

II. MAGNETOELASTIC ENERGY
Let us consider here the geometrical magnetoelastic configuration shown in Fig. 1. The primary field quantities for describing magnetic and elastic quantities are the magnetic flux density b and the elastic displacement field d.Ifthe system is nonhysteretic, then the particular state of the system is independent of the history. Correspondingly, the energy required to change the system from one particular state to another can be represented as a unique real number.
Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. as arguments. For the physics of magnetostriction and its relation to energy, see for instance [9] and [10].
The functional is always a measurable quantity. It is the work done by external mechanical forces or electrical power sources into the system. In this paper, we restrict the focus on static magnetoelastic systems, for which such an energy functional exists.
We define the displacement field d on the magnetoelastic material M only. We assume that there is enough air around M so that the deformation do not change the topology of our configuration. We say that (d) is the overall deformed configuration, whereas the deformation only occurs on M. With this notation, = (0) is the reference configuration, i.e., the configuration when d = 0.

III. VARIATIONS OF MAGNETIC FLUX DENSITY
Let us consider the small changes of the energy functional . We first focus on the variations δb of the argument b while d is held fixed. Notice, the fixed configuration implies that there is no elastic deformation work during the variation of b. The change δ b is defined as for small variations δb of the magnetic flux density.
According to the electromagnetic theory [11], the electric work δ W done by the current source W in changing the magnetic field from b to b + δb is Here, h is the magnetic field as it is defined classically in a rigid configuration, which now happens to be (d).
We require that the external electrical work corresponds to the change in the internal energy. That is, we have Using the magnetic relations where a is the magnetic vector potential and j s is the known source current, we are able to write the change δ b as We defined the h-field on the deformed configuration, the same way as in [12]. Therein, however, magnetostrictive effects were not considered. In general, the magnetic properties of M may change due to the deformation. This is to say, the h(b)-relation may depend on the mechanical deformation.
In [13] and [14], the dependence of the h(b)-relation on deformation is introduced via stress tensor. On the other hand, in [15] and [16], the dependence is due to the small strain tensor , that is In this paper, we assume that the magnetic constitutive relation can be expressed in the form of (7). We emphasize that this is a hypothesis to be verified experimentally, as will be discussed later on.
According to the well-established infinitesimal strain theory [17], the small strain tensor is properly defined in the reference (0) configuration, as However, (2) and (6), where the constitutive relation is needed, are written in the deformed configuration. We get around this difficulty by assuming so small displacements that the measurement of magnetic field can be assumed to take place in the reference configuration (0). Accordingly, we replace (d) by (0) in (2) and (6).

IV. VARIATIONS OF DISPLACEMENT
Let us then consider the change δ d , which we define as When the variation δd of the displacement occurs, external mechanical forces can do mechanical work to the system [18]. There are two kinds of external mechanical forces, the external surface forces f s and external body forces f b . The virtual work δ M done by the external mechanical forces can be written as We require that the change in the internal energy corresponds to the virtual work, that is As we have ignored hysteresis, we are free to choose any history for our configuration. From now on, we assume that the deformation of the system takes place first, and only then the magnetic field rises. The total energy is then where Here, the integrands ω e , ω m and ω m0 are defined as where ν 0 is the the reluctivity of vacuum. Object σ () is the stress tensor in the absence of magnetic field.
We call e and m by the name elastic and magnetic energy functionals, respectively. We emphasize that this is not the only possible choice. According to the discussion in Section III, we measure magnetic field in the reference configuration (0). We assume also that the plain elastic constitutive relation σ () of M is known. Now, we are able to write The change in the elastic energy δ d e can be written as Exploiting the Frechet derivative [19], the change δ d m can be written as when higher order terms are ignored. Operator ∂ d is the Frechet derivative with respect to the displacement field d. Combining (10), (11), (18), (19), and (20) gives which together with (8) [19]. If there are no external mechanical forces, then (21) reduces to

V. MAGNETIC FORCES
In order to solve (22), we need to get an expression for the magnetic force f m inside M. This can be given by Again, in principle, the magnetic forces should be evaluated in the deformed configuration (d). However, we assume small displacements and consequently the magnetic forces can be evaluated in the reference configuration, that is, Applying the virtual power principle to the magnetic energy is technically rather tricky. In particular, it requires taking into account the dependence of configuration on the displacement.
The virtual power principle for a general magnetostrictive material is applied in [20]. Accordingly, the magnetic force f m consists of three terms The first term is called the inhomogeneity term, the second one is the anisotropy term, and the third term is called the magnetostrictive term. According to [21], the magnetostrictive stress tensor σ M inside M is defined as and the divergence of the magnetostrictive stress tensor is the magnetostrictive force term.

VI. HELMHOLTZ FREE ENERGY DENSITY
We defined the stress tensor for the elastic system in the standard way. According to [21], it is possible to define an augmented Maxwell stress tensor σ T , whose divergence yields all the magnetic force components given in (24). It follows that the equilibrium equation of linear elasticity can be written as: This makes it possible to interpret the object σ me = σ + σ T as the stress tensor of the coupled magnetoelastic system. In [22], the Helmholtz free energy density ψ is introduced to define the coupled constitutive relations of magnetoelastic systems. The idea is further extended in [23]. The constitutive relations are constructed such that they can be expressed as partial derivatives of the (local) energy density. Finite element implementations employing this approach can be found in [24]- [26].
For a comparison, in [25], the energy density is chosen such that the magnetic constitutive relation can be expressed by where the double bracket notations denotes tensors of rank 2 and I refers to the identity tensor. Parameters μ 1 and μ 2 describe the magnetostrictive effects to be obtained by measurements. Equations (27) and (28) can be recognized as a special case of (7). Once the h(, b)-relation is measured, it is straightforward to construct the free energy density object that corresponds to the magnetic constitutive relation (7) by integration.
It is our strategic choice not to define the elastic constitutive relation for a magnetoelastic system. The reason is that it is not obvious how the coupled constitutive relation models the effects of the inhomogeneity force component in (24) on the boundary ∂ M. This is particularly an issue, when one employs (21) to solve the elastic equations. Consequently, we define only the elastic constitutive relation σ () for elastic systems without magnetization. Evidently, it is not clear whether the elastic material parameters are measurable in the absence of magnetic field for all materials.

VII. MEASUREMENT SETUP
In order to run finite element simulations for a magnetostrictive material M, we need to choose a discretization and a proper amount of test functions δd i and δa i in (6) and (21). Since we assume that the elastic constitutive relation σ () is measured a priori, all what is needed is to measure relation h(, b).
In principle, the measurement of h(, b) is straightforward. Employing external mechanical stress and current sources, one is able to change the values of (x), b(x), and h(x) on point x ∈ M. One then measures the triplet ((x), b(x), h(x)). After this, one changes the external stress or currents and measures the triplet again. Repeating the process iteratively, one should obtain enough data to construct the constitutive relation h(, b).
Notice that it is not necessary to obtain homogeneous b or fields within the sample M. If all the quantities , b, and h can be measured even at one point x, then this should be enough. In practice, however, homogeneous b-field may be favorable.
Measurement setups for simultaneous measurements of magnetostrictive strain, magnetic field, and flux density are presented for instance in [27]- [30].

VIII. EXAMPLES
In this section, we demonstrate with examples the work flow from measuring the h(, b) relation to the finite element analysis.

A. Measurement Data
The magnetostrictive measurement data was provided by Aalto University (A. Belahcen). The data are measured using a measurement setup similar to the one presented in [27]. The measurement data contain the magnetostrictive strain in the direction of the magnetic field inside some electrical steel. In addition, the values for magnetic field and magnetic flux density are included in the measurements.
From now on, we employ two Cartesian coordinate systems, the local unprimed and global primed systems, which are explained in Fig. 2. During the measurements, external mechanical stress is varied, and for each value of external stress, a hysteresis loop is measured in the (b 3 , 33 , h 3 )-space,  as demonstrated in Fig. 3. The positive external stress is tensile, while negative is compressive.
During the measurement of each hysteresis loop, there is an offset strain due to the mechanical external stress, which is not measured. Only the magnetostrictive strain relative to the externally varied configuration is included in the measurements. However, the constitutive relation h(, b) requires the total strain , which is the sum of the mechanical and magnetostrictive strains.
We can recover the mechanical strain by employing the hypothesis presented in [30]; the magnetization reaches M s at high magnetic field, the domain configuration, and thus the values of magnetostriction are strictly identical whatever the stress level.
The magnetization M as a function of h is plotted in Fig. 4. We set M s = 1.35 × 10 6 A/m. As a next step, we plot (M, )-graphs in a way that they coincide at M = M s (see Fig. 5). From Fig. 5, we recover the elastic offsets for strain. Finally, we are able to plot all the measurements in the (b e , 33 , h 3 )-space (see Fig. 6).
We plot the reluctivity ν 33 = h 3 /b 3 as a function of 33 and b 3 (see Fig. 7). We observe that the reluctivity ν 33  does not depend on b 3 when |b 3 | < 0.3 T, but it depends on 33 . Moreover, when | 33 | < 0.5 μ m/m the dependence ν( 33 ) can be expressed as The values for the constants are ν = −10.2 × 10 6 m/A and ν M = 20.5 m/A. Since the measurements only include the 33 -component of the strain tensor, we have to assume that the other components vanish. Then (28) reduces to 33 (30)  and consequently, the measurements seem to define the sum of the parameters μ 1 and μ 2 . Equations (29) and (30) look very similar, but unfortunately (29) is defined for the reluctity tensor, whereas (30) is defined for the permeability tensor, which makes their direct comparison bit difficult.

B. Magnetostrictive Stress Tensor
The measurements do not give any information on the possible anisotropy of the sample. Hence, we are forced to use an isotropic model for the reluctivity tensor. In the unprimed coordinate system, the reluctivity tensor ν can be written as According to [20], the magnetostrictive stress tensor σ M defined in (25) can be written using the Einstein summation convention for magnetic linear materials by Accordingly, we have in the local unprimed coordinates, where b 2 is the square of the norm of the b-field.

C. Change of Coordinates
The transformation tensor S i j between the primed and unprimed coordinate system depends on b and it is defined as In particular, the transformation of the relucticity and magnetostrictive stress tensors from local to global coordinates are and generally, every component can have nonzero values. For the inverse transformations, the transpose of S i j is employed, for instance The divergence of magnetostrictive stress tensor in the global coordinates is

D. Numerical Solution
We use Elmer [31] for numerical solutions of the coupled magnetomechanical problem. Elmer is an open-source FEM software for multiphysical problems and it contains solvers for linear elasticity and 3d-quasistatic magnetic fields [32]. We use the MagnetoDynamics and StressSolve-solvers in order to solve magnetic and the displacement fields of linear elasticity. Elmer uses Gauss-Seidel iterations to solve coupled multi-physical problems. The number of steady state iterations required to solve a particular multi-physical problem depends on how tightly the physical models are coupled together [33].
The MagnetoDynamics solver uses edge elements w i to represent the vector potential a, such that a = a i w i , where the coefficients a i are constants. Consequently, the magnetic As fields between solvers are usually passed as nodal fields in Elmer, we wish to represent b with linear nodal basis functions λ i , such that This can be done in the weak sense by requiring which leads to a linear system, from which the coefficients b λ i can be solved. For elastic analysis, we assume that the sample is supported on the faces F 1 and F 2 , whose normals are out of and into the page. Correspondingly, the boundary condition d(x) = 0 for x ∈ F 1 ∪ F 2 is used for the elastic solution. Fig. 9. Magnetic flux density is plotted in the solid parts on the symmetry plane. The white lines represent the direction of b. We observe that the condition |b| < 0.3 T holds within the sample M. We also note that the b-field is not homogeneous inside M.
We now have the magnetic flux represented as nodal field within the magnetostrictive material M. It follows immediately that each component of the transformation tensor S i j has nodal representation. Furthermore, the magnetostrictive stress and reluctivity tensors have nodal representation in primed and unprimed coordinate systems for each component. The strain field is represented as nodal field by default in Elmer.
On one hand, the magnetic field calculation is coupled to the elastic solution by the strain-dependent reluctivity tensor. On the other hand, the elastic solution is affected by the magnetic forces. The strain-dependent reluctivity tensor is straightforward to implement into Elmer, but the magnetic forces require some care.
Elmer includes nodal magnetic forces [29] for isotropic nonlinear magnetic materials. These correspond to the inhomogeneity forces in (24). As we use the isotropic model, the anisotropy term of (24) is zero. However, we need to add  Deformation of the magnetostrictive sample is shown, when the contribution of the magnetostrictive stress tensor is turned OFF. The displacements are now scaled by 30 × 10 6 . The local strain 33 is shown as colormap on ∂ M with the units μ m/m. the contribution of the magnetostrictive stress tensor when the magnetic forces are passed to the elasticity solver. Alternatives to obtain the nodal forces can be found for instance in [25].
We need also to evaluate the divergence of σ M . According to (38), the divergence is well defined within the elements in the finite element mesh, as we employ the nodal representation for the magnetostrictive stress tensor components. The StressSolve-solver requires integrals of the body force, which are straightforward to evaluate.

E. Test Case
Our test material sample M is a cube. The sample is part of a magnetic circuit, as described in Fig. 8. The value of the total coil current is adjusted such that conditions |b| < 0.3 T and | 33 | < 0.5 μ m/ m hold. We set the the value 193.053 GPa for Young's modulus and the value 0.29 for the Poisson ratio for M. Fig. 9 shows the magnetic field of the magnetic circuit and Fig. 10 shows the deformation due to the magnetic forces. The solution of the problem took 16 steady state iterations. For a comparison, the deformation without the effects of magnetostrictive stress tensor is presented in Fig. 11. In this case, the solution took three iterations, because also the inhomogeneity forces depend on the magnetostriction.
We observe that magnetostrictive stress tensor dominates the deformation. The magnetostrictive force in the sample is tensile in the direction of the magnetic flux, as expected from measurements.

IX. CONCLUSION
We assumed the magnetic constitutive relation is of the  form h(, b). We also assumed that the purely elastic constitutive relation σ () of the material is known or can be specified from measurements. The analysis shows that the simultaneous measurements of h, b, and are sufficient for complete finite element analysis of the magnetostrictive material.
The approach presented relies on a certain decomposition of magnetoelastic energy. In other words, since hysteresis is ignored, we were able to choose any magnetoelastic history for the material. This implies the strain-dependent magnetic constitutive relation needs to be measured, whereas the elastic constitutive relation can be defined in a standard way without the magnetic field inside the magnetostrictive material.
We demonstrated corresponding finite element simulations with examples employing the open source finite element software Elmer. The solution process converged in all cases and the computed strain was in the expected direction and in align with the measurements.
The constitutive relation h(b, ), condensing all microscopic interactions onto mesoscopic level, has to be specified with experiments. Accordingly, to further test and verify the presented approach a measurement setup and its finite element model is needed. Then, for a given material the values of , b, and h should be specified iteratively between measurements and computations seeking for a converging sequence. If the material parameters obtained and the model of magnetostriction can be successfully exploited to model new corresponding cases, this then verifies both the measurements and the employed approach to magnetostriction.