Construction and Numerical Realization of a Magnetization Model for a Magnetostrictive Actuator Based on a Free Energy Hysteresis Model

Featured Application

This study constructs and analyzes the constitutive model of the giant magnetostrictive material actuator, which provides a theoretical basis for the development of a new intelligent structure of giant magnetostrictive material.

Abstract

Giant magnetostrictive actuators (GMA) driven by giant magnetostrictive material (GMM) has some advantages such as a large strain, high precision, large driving force, fast response, high reliability, and so on, and it has become the research hotspot in the field of microdrives. Research shows there is a nonlinear, intrinsic relationship between the output signal and the input signal of giant magnetostrictive actuators because of the strong coupling characteristics between the machine, electromagnetic field, and heat. It is very complicated to construct its nonlinear eigenmodel, and it is the basis of the practical process of giant magnetostrictive material to construct its nonlinear eigenmodel. Aiming at the design of giant magnetostrictive actuators, the magnetization model based on a free-energy hysteresis model has been deeply researched, constructed, and put forward by Smith, which combines Helmholtz–Gibbs free energy and statistical distribution theory, to simulate the hysteresis model at medium or high driving strengths. Its main input and output parameters include magnetic field strength, magnetization, and mechanical strain. Then, numerical realization and verification of the magnetization model are done by the Gauss–Legendre integral discretization method. The results show that the magnetization model and its numerical method are correct, and the research results provide a theoretical basis for the engineering application of giant magnetostrictive material and optimized structure of giant magnetostrictive material actuators, which have an important practical application value.

1. Introduction

In the field of intelligent manufacturing, micro-displacement and microactuation technology with an actuation stroke of less than $1\mathrm{m}\mathrm{m}$ and a resolution of less than $1\mathsf{\mu }\mathrm{m}$ has a wide range of applications [1,2]. At present, there are some functional materials used in the field of microactuation technology such as piezoelectric materials, shape memory alloys, and giant magnetostrictive material [3,4,5,6,7]. As a new type of functional material, giant magnetostrictive material have some advantages of a strong magnetostrictive effect, a high electromechanical coupling coefficient, high response speed, noncontact drive, and so on, and they are widely used in aerospace, precision control, and intelligent manufacturing fields.

Giant magnetostrictive material (GMM) is ferromagnetic material that were developed by Clark and Belson et al. in 1971[8]. Currently, the type of GMM used most commonly is ${\mathrm{TB}}_{\mathrm{X}}{\mathrm{Dy}}_{1-\mathrm{x}}{\mathrm{Fe}}_{\mathrm{y}}$, when $x=0.23$∼$0.5$ and $y=1.9$∼$2.0$. The trade name of GMM is Terfenol-D [8,9].

The giant magnetostrictive actuator (GMA) developed by the giant magnetostrictive material (GMM) has the characteristics of a large magnetostriction coefficient, large output power, fast response speed, and high magnetic(electric)–mechanical conversion efficiency, and it is widely used in the fields of ultraprecision machining, micromotors, vibration control, fluid machinery, sonar, and so on. It is one of the hotspots in the field of microactuation [10].

Kobayashi, who worked for Toshiba, had developed a new type of tool microfeeding device with a positioning accuracy on the nanometer scale. Tsinghua university [11] designed a linear peristaltic mechanism used the giant magnetostrictive material with a displacement accuracy of $\pm 15 \mathrm{n}\mathrm{m}$. Successful development of this actuator has promoted the development of research fields such as Micro Electro Mechanical System (MEMS).

Under the action of the external driving magnetic field, the giant magnetostrictive material rod will output displacement and force, and the input electrical signal and output displacement signal of the actuator have corresponding eigen relations. Due to the nonlinear strong coupling properties between the machine electromagnetic heat of the material, it is very difficult to construct a magnetic–mechanical coupling model. Clark [12] proposed a piezomagnetic constitutive equation of a giant magnetostrictive actuator based on a large volume of research, which became the basic equation for the study of magnetic–mechanical coupling characteristics of the giant magnetostrictive actuator. Benbouzid[13] used finite element simulations to calculate nonlinear dynamic characteristics of the giant magnetostrictive rod, which can guide the design and optimization of the magnetic circuit of the giant magnetostrictive actuator. However, the model is a two-dimensional planar model, and its dynamic characteristics are obtained from static characteristics of the giant magnetostrictive material, in addition, there are too many parameters, due to static modeling, which lead to complicated calculations and susceptibility to other field-related parameters, therefore, application of the model has certain limitations.

In the research of multiphysics dynamic magnetic–mechanical coupling models, Azoum [14] established a three-dimensional generalized finite element model based on the magnetic machine coupled constitutive equation, and simulations of the solenoid coil driving the GMM rod were made. At the same time, Benatar [15] used the multiphysics coupling software FEMLAB (Fimite Element Modeling Laboratory) (version: 3.1.0.157, COMSOL Inc., Stockholm, Sweden,1986) to establish a three-dimensional electromagnetic–machine coupling model of GMM transducers, which is a miscellaneous shareware software developed by COMSOL Inc.. Zhao of Zhejiang University [16] also used the weak solution of COMSOL software (version 5.3, COMSOL Inc., Stockholm, Sweden, 1986) to simulate the electromagnetic–mechanical coupling model of intelligent mast components designed with GMM. However, the coefficient matrix substituted by it is a constant coefficient matrix, and the simulation results have a certain error compared to the actual experimental results. None of the methods construct the model from the characteristics of the giant magnetostrictive material itself. A reasonable research method is to consider the intrinsic nonlinearity and hysteresis characteristics of giant magnetostrictive material when the output and input eigen relations of actuators is studied.

2. Magnetostrictive Mechanism of the Giant Magnetostrictive Material and Its Modeling Method

2.1. Ferromagnetic Properties of Material and Their Magnetostrictive Mechanisms

(1) Ferromagnetic principle of material

In an atom of matter, electrons orbit around a nucleus to produce orbital magnetic moments, at the same time, electrons also produce a spin magnetic moment by rotating around their own fixed axis [17], which produces ferromagnetic properties of the material.

Generally, the electrons of the ferromagnetic material atoms are not arranged in an ideal order. There are unfilled electron layers inside the atoms, the spin magnetic moment of the electrons is not neutralized, the atoms have a “permanent magnetic moment”, and the spin magnetic moment in the ferromagnetic material can be fixed in a small local area, this process is called spontaneous magnetization. The spontaneous magnetization zone is called a magnetic domain, and the boundary layer of adjacent magnetic domains is called a magnetic domain wall. If an external magnetic field is applied to the ferromagnetic material, the arrangement of some of the magnetic domains changes along the direction of the external magnetic field, which shows magnetic properties. In addition, the degree of magnetization of ferromagnetic materials is also related to temperature: when the temperature exceeds a certain critical value (Curie point), the ferromagnetism of the ferromagnetic material disappears, and the ferromagnetic material is changed to paramagnetic.

(2) Magnetostrictive mechanism

Specifically, there are three reasons to change the magnetization state of the ferromagnetic material by changing the shape and volume of ferromagnetic materials [18,19]. The main effect of magnetostrictive strain on giant magnetostrictive material is field-induced deformation [20].

① Spontaneous deformation caused by spontaneous magnetization

Spontaneous magnetization is caused by an exchange force independent of the external magnetic field. It is imaginary that there is a single domain crystal in the giant magnetostrictive material, which is spherical above the Curie temperature ($380 °\mathrm{C}$). When it is cooled from the Curie temperature, the exchange force causes the crystal to spontaneously magnetize, and the shape of the crystal changes (shown in Figure 1). However, since the directions of spontaneous deformation of the magnetic domains inside the giant magnetostrictive material are random, there is no obvious elongation or shortening in the macroscopic direction.

② Field-induced deformation caused by technical magnetization

When the strength of the external magnetic field gradually increases from zero, the magnetization of the giant magnetostrictive material gradually rises to a saturated state, then a magnetostriction phenomenon occurs, this process is called field-induced deformation. The field-induced deformation process can be divided into four stages: the reversible magnetostriction stage, the irreversible magnetostriction stage, the stage where the rotation of the magnetic moment of the magnetic domain plays a major role, and the magnetization asymptotic saturation phase.

After saturation of magnetostriction is reached, and the magnetic field strength $H$ of the external magnetic field is lowered, even if the applied magnetic field strength $H$ is reduced to zero, the magnetic domain cannot be completely restored to the initial state, and the magnetostriction rate $\lambda$ cannot be completely zeroed, this results in a certain residual magnetization inside the giant magnetostrictive material, which is called hysteresis.

2.2. Factors Affecting Magnetic Coupling Characteristics of Giant Magnetostrictive Material

The giant magnetostrictive material is magnetic anisotropy, and its saturation magnetostriction coefficient is changed along with the changing of the direction of the magnetic field. The magnetic domain distribution in the process of technical magnetization has dual effects of wall shift and domain rotation [21,22], which makes it extremely difficult to study the intrinsic properties of giant magnetostrictive material. According to references [23,24], a relationship between the magnetic field and magnetostriction coefficient of different crystal orientations for ${\mathrm{Tb}}_{0.3}{\mathrm{Dy}}_{0.7}{\mathrm{Fe}}_{1.95}$ has been obtained at 20 $°\mathrm{C}$ and 11 $\mathrm{M}\mathrm{p}$.

When magnetostriction is saturated, there is a certain hysteresis between the magnetostriction and the magnetic field strength, and its relationship is nonlinear. Because of the magnetostrictive effect and magnetostrictive inverse effect inside the giant magnetostrictive material, the two subsystems independent of each other between the magnetic system and the mechanical system are coupled [25,26]. At this time, if the giant magnetostrictive material is operated at a constant temperature, and only the strain in the longitudinal direction is considered, the magnetic field variables (magnetic field strength $H$, magnetic induction $B$) and mechanical field variables (stress $\sigma$, strain $\epsilon$) are correlated with each other. Assuming that the magnetic field strength $H$ and stress $\sigma$ are independent variables, the coupling relationship between the magnetic field and the mechanical field is expressed as a linear piezomagnetic equation [12].

$$\left\{\begin{array}{c}\epsilon =E^H\sigma +d_{33}H\\ B=d_{33}\sigma +\mu H\end{array}\right.$$

(1)

where $\epsilon$ is the total strain in the length direction of the GMM material. $E^H$ is Young’s modulus when the magnetic field strength $H$ is a constant. $d_{33}$ is pressure magnetic coefficient. $\sigma$ indicates the magnitude of stress. $\mu$ indicates the magnetic permeability when the stress is a constant.

According to Equation (1), the magnetostrictive effect is to make the magnetic system of the giant magnetostrictive material transfer energy to the elastic system. and the magnetostrictive inverse effect is to make the elastic system of the giant magnetostrictive material transfer energy to the magnetic system., simultaneously, energy exchange is achieved between the magnetic system and the elastic system of the material [27,28].

In the actual alternating magnetic field, the magnetic induction intensity $B$ of the giant magnetostrictive material is closer to one phase than the applied alternating magnetic field $H$ due to factors such as domain wall resonance, the hysteresis effect, and natural resonance. Thus, as the giant magnetostrictive material is continuously magnetized in the alternating magnetic field, the applied energy is continuously consumed, which causes hysteresis loss.

In addition, the giant magnetostrictive material also has ohmic loss, compressive stress characteristics [29], temperature characteristics [30], $\Delta E$ effects [31], electromechanical equivalent characteristics [32], and so on, these characteristics affect the hysteresis characteristics of giant magnetostrictive material and construction of the microactuator magnetic coupling model.

2.3. Comparative Study of Hysteresis Models of Giant Magnetostrictive Actuators

At present, there are two types of models mainly used to study the hysteresis phenomenon of giant magnetostrictive actuators [33]: hysteresis models based on mathematical methods (Preisach model) and hysteresis models based on physical methods (Jiles–Atherton model, free-energy hysteresis model).

(1) Preisach model

The Preisach hysteresis model [33] is a model for calculating magnetostrictive deformation by studying the hysteresis model of iron material proposed by Restorff and Clark. Subsequently, Tan proposed a Preisach hysteresis model describing frequency hysteresis behavior at $200 \mathrm{H}\mathrm{z}$ working range, which shows the dynamic characteristics of hysteresis behavior with increasing operating frequency [34]. At present, the research on the Preisach model has become comprehensive [35,36,37]. Because the model can only reflect the relationship between input and output, it cannot explain the magnetostrictive mechanism inside the giant magnetostrictive material, and it cannot reflect the relationship between magnetostriction and the external thermal and machinery fields, moreover, a large number of complex equations and nonphysical parameters need to be identified and processed in the model, which leads to low flexibility and a long processing cycle for the Preisach model in the actual application process.

(2) Jiles–Atherton model

The Jiles–Atherton hysteresis model [38] is founded based on the domain wall theory of ferromagnetic materials, this model is applied to build the dynamic hysteresis model of Terfenol-D by Dapino [39]. The Jiles–Atherton model is quite limited for the capture of hysteresis, and many physical causes of hysteresis have not been explored. It can only be applied to research a small number of ferromagnetic functional materials such as giant magnetostrictive material and shape memory alloys.

(3) Free energy hysteresis model

In 2003, Smith and M.J Dapino [40] proposed a free energy hysteresis model based on the Jiles–Atherton model, the model is combined Helmholtz–Gibbs free energy and statistical distribution theory to simulate the hysteresis model at medium or high driving strengths, and its main input and output parameters are included with magnetic field strength, magnetization, and mechanical strain. The Boltzmann statistical theory is adopted in the model, and the related energy relation function of the external thermal field and the mechanical field is substituted into the telescopic deformation of the material, furthermore, an intrinsic nonlinearity three-field model coupled by magnetic field, machanical field and thermal field is established, and it can be built a model accurately for the functional materials such as Terfenol-D and piezoelectric ceramics [41,42,43].

Since the free energy hysteresis model is a new hysteresis model proposed in recent years, it has the advantages of being a simple model, easily modifiable external influence factors (such as eddy current, temperature, etc.), and less parameters, and it is the main theoretical basis for studying the magnetic coupling model of giant magnetostrictive microactuators.

3. Construction of a Hysteresis Model of a Giant Magnetostrictive Microactuator Based on Free Energy

3.1. Research Process Based on the Free Energy Hysteresis Model

The conditions for the free energy hysteresis model are established when the temperature is constant and the influence of eddy current loss is ignored, it works in a low frequency range, the pre-pressure applied to the giant magnetostrictive microactuator is assumed to be large enough to ignore the magnetic anisotropy, and so on. to establish the free energy hysteresis model of the giant magnetostrictive actuator include three steps, the first, the Helmholtz–Gibbs free energy relationship for homogeneous materials with a constant internal field is established, then, the heterogeneity of the actual material and the non-constantness of the internal effective field are considered, and a random distribution function is substituted to derive the hysteresis relationship model between the applied magnetic field strength $H$ and the magnetization $M$. It is also important to study its magnetic machine coupling characteristics and construct the nonlinear relationship between the applied magnetic field strength $H$ and strain $\epsilon$. The research process is shown in Figure 2.

3.2. Theoretical Basis for the Establishment of the Free Energy Hysteresis Model

(1) Magnetization process of the giant magnetostrictive material

In order to establish the energy model of the material, it is necessary to analyze the lattice structure of the material by means of vector representation. In the paper, the symbol ‘()’ represents a plane, the symbol ‘[]’ represents a vector, and the symbol ‘〈 〉’ represents the general name of the series represented by the mark. Therefore, the various planes of the cube are represented as: $\left(100\right)$, $\left(010\right)$, $\left(001\right)$, $\left(\overline{1}00\right)$, $\left(0\overline{1}0\right)$, and $\left(00\overline{1}\right)$, and the direction vector is expressed as: [$100$], [$010$], [$001$], [$\overline{1}00$], [$0\overline{1}0$], and [$00\overline{1}$], among them,$\overline{1}$ indicates a negative direction.

In the current processing conditions, the lattice of Terfenol-D in the [$11\overline{2}$] direction grows in the form of a tree (as shown in Figure 3). In room temperature, the easy magnetization axis of Terfenol-D is generally in the $〈111〉$ direction series, and the maximum strain occurs when the magnetization $M$ is rotated from direction [111] to [$11\overline{1}$] or from direction [$\overline{111}$] to [$\overline{11}1$].

To take a sample that is not magnetized (as shown in Figure 4). Figure 4a shows that all regions in the specimen been not magnetized have nonzero spontaneous magnetization $M_0$. Figure 4b shows that the change in magnetization is mainly caused by the reversible domain wall motion when the input field is low. Figure 4c shows that two irreversible mechanisms are appearanced when the input field gradually increases, which a magnetic domain is increased caused by domain wall movement (consistent with the direction of the magnetic field), and the magnetic moment is rotated to the direction of easy magnetization [$11\overline{1}$]. In this phase, inputting a small change in the driving magnetic field produces a large change in magnetization and output mechanical strain. Figure 4d shows that the magnetization saturation state of the material is reached. At this time, the sample can be approximated as a single magnetic domain, and the magnetic moment is completely rotated from the easy magnetization direction to be consistent with the magnetic field direction.

In this study, it is assumed that the pre-stress of the material perpendicular to the direction of the magnetic moment is large enough to dominate the anisotropy of the lattice, then, the magnetization direction of the magnetic domain is easily changed from the original direction $〈111〉$ to the [111] and [$\overline{1}\overline{1}\overline{1}$] directions perpendicular to the axis [$11\overline{2}$]. This can be considered that the output strain is mainly caused by the rotation of the magnetic moment:

$$\lambda =\frac{3}{2}{\lambda }_s{(\frac{M}{M_s})}^2,$$

(2)

where ${\lambda }_s$ is the saturation magnetostriction rate. $M_s$ is the saturation magnetization. $M$ is magnetization. $\lambda$ is the magnetostriction rate.

In order to simulate the hysteresis behavior of a giant magnetostrictive actuator, the nonlinear relationship between the input magnetic field $H$ and the magnetization $M$ needs to be described quantitatively. Then, Equation (2) is introduced to derive the relationship between the input magnetic field $H$ and the output strain $\epsilon$.

(2) Description function of Helmholtz and Gibbs energy

To describe the relationship between magnetostatic and magnetoelastic properties of a giant magnetostrictive actuator, it is assumed that free energy is a function of temperature $T$ and other ordered parameters, where $e$ represents magnetization $M$ in magnetostatics and strain $\epsilon$ in magnetoelasticity. The external field is represented by $\tilde{\varphi }$ and is conjugated with $e$ in thermodynamics. When $M$ is represented by the ordered parameter $e$, the external field is represented as a magnetic field $H$. When $\epsilon$ is represented by the ordered parameter $e$, the external field is represented as stress $\sigma$. $\psi (e,T)$ is used to express the general relationship of Helmholtz’s free energy. When no external field is added, the condition of thermodynamic equilibrium is that $e$ is at a minimum, which is:

$$\varphi (e,T)=\frac{\partial \psi (e,T)}{\partial e}=0.$$

(3)

The energy response of the system can be expressed as $\varphi$. When the system is acted on by the external field $\tilde{\varphi }$, total free energy can be expressed as Equation (4):

$${\psi }_{\tilde{\varphi }}(e,T)=\psi (e,T)-\tilde{\varphi }e.$$

(4)

The equilibrium conditions applied by an external field are expressed by Equation (5):

$$\varphi (e,T)=\tilde{\varphi }.$$

(5)

In the case where the pre-pressure is sufficiently large, two magnetization directions in [111] and [$\overline{1}\overline{1}\overline{1}$] in the giant magnetostrictive material will be easily produced. It is necessary to consider both the internal energy and the magnetic moment energy of the material in these two magnetization directions. It is assumed that the interaction between magnetic moments is adiabatic, and the temperature is stable below the Curie temperature, the Helmholtz potential energy $\psi$ can be expressed by Equation (6).

$$\psi (M,T)=\frac{H_h{M}_s}{2}[1-{(M/M_s)}^2]+\frac{H_{h}T}{2T_c}[M\ln (\frac{M+M_R}{M_s-M})+M_s\ln (1-{(M/M_s)}^2)],$$

(6)

where $H_h$ is the biasing magnetic field, $M_s$ is saturation magnetization, and $T_c$ is the Curie temperature.

According to statistical mechanics, Helmholtz free energy can be represented by piecewise quadratic relations under isothermal conditions.

$$\psi (M)=\left\{\begin{array}{l}\frac{1}{2}\eta {(M+M_R)}^2,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.33em}\hspace{0.33em}M\le -M_I\hfill \\ \frac{1}{2}\eta {(M-M_R)}^2,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.33em}\hspace{0.17em} M\ge M_I\hfill \\ \frac{1}{2}\eta {(M_I-M_R)}^2(\frac{M^2}{M_I}-M_R),\hspace{1em}\left|M\right|<M_I\hfill \end{array}\right.,$$

(7)

where $M_I$ is the magnetization generated at the inflection point, $M_R$ is residual magnetization, $\eta$ is the converted slope, and $\eta =\frac{dH}{dM}$.

By Equation (6), Gibbs energy can be written by Equation (8):

$$G(H,M,T)=\psi (M,T)-HM.$$

(8)

(3) Hysteresis model of homogeneous material

For homogeneous materials, the effective field $H_e$ is equal to the application field $H$, and the average magnetization can be written by Equation (9).

$$\overline{M}=x{}_{+}〈M_{+}〉+x{}_{-}〈M_{-}〉,$$

(9)

where $x{}_{+}$/ $x{}_{-}$ is the probability of positive or negative magnetic moment, and $M_{+}$/ $M_{-}$ is the expected value of magnetization. $〈M_{+}〉$ and $〈M_{-}〉$ can be obtained by Equation (10):

$$\begin{array}{l}〈M_{+}〉={\displaystyle {\int }_{M_I}^{\infty }M\mu (G)dM}\\ 〈M_{-}〉={\displaystyle {\int }_{-\infty }^{-M_I}M\mu (G)dM}\end{array}$$

(10)

where:

$$\mu (G)=C{e}^{-GV/kT}.$$

(11)

The probability of Gibbs energy obtained is expressed quantitatively in the above equation, where $k$ is the Boltzmann constant, $C$ is a constant, and this constant can be selected according to the magnetization integral value 1, $V$ is lattice volume. Considering the Boltzmann energy balance, the approximate relaxation process characteristics of the material can be obtained in the lattice volume $V$, and then the constant $C$ can be estimated to obtain the average magnetization value.

$$\begin{array}{l}〈M_{+}〉=\frac{{\displaystyle {\int }_{M_I}^{\infty }M{e}^{-G(H,M)V/kT}dM}}{{\displaystyle {\int }_{M_I}^{\infty }{e}^{-G(H,M)V/kT}dM}}\\ 〈M_{-}〉=\frac{{\displaystyle {\int }_{-\infty }^{-M_I}M{e}^{-G(H,M)V/kT}dM}}{{\displaystyle {\int }_{-\infty }^{-M_I}{e}^{-G(H,M)V/kT}dM}}\end{array}$$

(12)

According to the normalized equation $x_{+}+x_{-}=1$, the matrix component can be quantitatively described by the differential Equation (13).

$${\dot{x}}_{+}=-p_{+-}{x}_{+}+p_{-+}(1-x_{+}).$$

(13)

The likelihood that the magnetic moment changes from the positive direction to the negative direction is expressed by Equation (14).

$$\begin{array}{l}{p}_{+-}=\frac{1}{\Gamma (T)}\frac{{\displaystyle {\int }_{M_I-\in }^{M_I}{e}^{-G(H,M)V/kT}dM}}{{\displaystyle {\int }_{M_I-\in }^{\infty }{e}^{-G(H,M)V/kT}dM}}\\ p_{-+}=\frac{1}{\Gamma (T)}\frac{{\displaystyle {\int }_{-M_I}^{-M_I+\in }{e}^{-G(H,M)V/kT}dM}}{{\displaystyle {\int }_{-\infty }^{-M_I+\in }{e}^{-G(H,M)V/kT}dM}}\end{array}$$

(14)

where $\in$ is constant, $\Gamma$ is the relaxation time of the material, and $\omega =\frac{1}{\Gamma }$, which indicates the frequency of magnetic moment conversion.

In order to qualitatively analyze the model and quantitatively analyze the thermal motion region, it is necessary to simplify the complex behavior of the average magnetization $\overline{M}$ in the region where the thermal motion is negligible, according to the equilibrium conditions by Equation (15).

$$\frac{\partial G}{\partial M}=0.$$

(15)

Combined with (8), Equation (16) is obtained.

$$\frac{\partial H}{\partial M}=\frac{{\partial }^2\psi }{\partial M^2}.$$

(16)

It can be seen from the above equation that the slope of the hysteresis kernel in this linear domain is $1/\eta$, and the inflection point $M_I$ is the key point. At this time, the resilience is the largest, and the original estimated parameter value of the free magnetic hysteresis model of the giant magnetostrictive actuator can be obtained according to the above equilibrium conditions.

In the case where magnetization occurs for a short time and the thermal motion is negligible, the jump will be completed in an instant, so the magnetization can be simplified by an asymptotic relationship (as shown in Figure 5). $M_{\min }$ can be obtained by solving Equation (15), and $\overline{M}=M_{\min }$, then, $\overline{M}$ is calculated.

Combining with the Helmholtz free energy expression (7), the magnetization can be obtained.

$$[\overline{M}(H;H_c,\xi )](t)=\left\{\begin{array}{l}[\overline{M}(H;H_c,\xi )](0)\hspace{1em},\hspace{1em}\tau (t)=\varphi \hfill \\ \frac{H}{\eta }-M_R\hspace{1em},\hspace{1em}\tau (t)\ne \varphi \hspace{0.17em}and\hspace{0.17em}H[\max \tau (t)]=-H_c\hfill \\ \frac{H}{\eta }+M_R\hspace{1em},\hspace{1em}\tau (t)\ne \varphi \hspace{0.17em}and\hspace{0.17em}H[\max \tau (t)]=H_c\hfill \end{array}\right.,$$

(17)

When the transition occurs, the set of times $\tau (t)$ can be expressed by Equation (18).

$$\tau (t)=\left\{t\in (0,T_f]|H(t)=-H_c\hspace{0.17em}or\hspace{0.17em}H(t)=H_c\right\}.$$

(18)

The initial value of the magnetic moment $[\overline{M}(H;H_c,\xi )](0)$ can be expressed by Equation (19).

$$[\overline{M}(H;H_c,\xi )](0)=\left\{\begin{array}{l}\frac{H}{\eta }-M_R\hspace{1em},\hspace{1em}H(0)\le -H_c\hfill \\ \xi \hspace{1em},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-H_c<H(0)<H_c\hfill \\ \frac{H}{\eta }+M_R,\hspace{1em}H(0)\ge H_c\hfill \end{array}\right..$$

(19)

Coercivity $H_c$ can be expressed by Equation (20).

$$H_c=\eta (M_R-M_I).$$

(20)

(4) Hysteresis model of heterogeneous materials

The magnetization model (Equation (17)) is established on the premise that the internal lattice and magnetic domain structure of the material are completely uniform. In fact, the non-uniformity of the lattice structure of the material is led by the defects of the material itself and the non-uniform free energy profile of its different regions, according to the hypothesis that the effective field $H_e$ is equal to the application field $H$ in the homogeneous material, the magnetic interaction and Weiss mean field effect are ignored. To assume that the parameter $H_c=\eta (M_R-M_I)$ is a normal distribution, with respect to the average coercivity ${\overline{H}}_c$, the total magnetization can be written by Equation (21).

$$M(H)={\displaystyle {\int }_0^{\infty }\overline{M}(H;H_c,\xi )}f(H_c)d{H}_c,$$

(21)

where $H_c$ is the density distribution function and can be expressed by Equation (22).

$$f(H_c)=C_1{e}^{-{(H_c-{\overline{H}}_c)}^2/b},$$

(22)

where $C_1$ and $b$ are coefficients, and $\overline{M}$ can be obtained from Equation (10) or (17). In Equation (21), its lower limit of integration is 0, which indicates that the hysteresis core width must be non-negative. In addition, it is defined as a logarithmic function to be also used to reflect the non-negative nature of the hysteresis kernel width.

In the Jiles–Atherton model, the Weiss field effect caused by the coupling between domains can be represented by Equation (23).

$$H_e=H+\alpha M,$$

(23)

where $\alpha$ is a constant mean field parameter.

Weiss fields quantitatively describe the interaction between atoms in giant magnetostrictive material. The effective field $H_e$ exhibits changes caused by uneven magnetic moment distribution, and it is assumed that the effective field $H_e$ is a normal distribution of the applied field. The magnetization can be expressed by Equation (24).

$$M(H)={\displaystyle {\int }_{-\infty }^{\infty }{C}_2\overline{M}({\rm H};H_c,\xi )}{e}^{-{(H-{\rm H})}^2/\overline{b}}d{\rm H}.$$

(24)

Thus, under the effective field of the low-frequency range without considering the eddy current loss, the complete magnetization model of the heterogeneous polycrystalline material can be expressed as:

$$[M(H)](t)=C{\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }[\overline{M}({\rm H}+H;H_c,\xi )](t)e^{-{{\rm H}}^2/\overline{b}}}{e}^{-{(H_c-{\overline{H}}_c)}^2/b}d{\rm H}}d{H}_c.$$

(25)

(5) Construction of a magnetic–mechanical coupling eigen model of giant magnetostrictive material

The Gibbs energy expression of Equation (8) combines the internal energy and magnetic moment energy of the isotropic material on the magnitude of the magnetic domain, but the magnetic coupling effect is ignored, and the magnetostrictive ability of the material cannot be exhibited. When anisotropy is exhibited for the lattice dominated by the actuator’s pre-stress, the Helmholtz free energy relation embodying magnetoelasticity can be expressed by Equation (26).

$${\psi }_e(M,\epsilon )=\psi (M)+\frac{1}{2}{Y}^M{\epsilon }^2-Y^M\gamma \epsilon M^2.$$

(26)

The corresponding Gibbs energy is expressed by Equation (27).

$$G(H,M,\epsilon )=\psi (M)+\frac{1}{2}{Y}^M{\epsilon }^2-Y^M\gamma \epsilon M^2-HM-\sigma \epsilon ,$$

(27)

where $\psi$ can be obtained from Equation (7). $Y^M$ is the Young’s modulus at a certain magnetization. $\gamma$ is the magnetic coupling coefficient.

In the case of intense thermal motion, the magnetization $\overline{M}$ can be obtained by Equation (9), and the Gibbs energy can be obtained by Equation (27). Magnetization is expressed by Equation (28) when the thermal motion is ignored.

$$[\overline{M}(H,\epsilon ;H_c,\xi )](t)=\left\{\begin{array}{l}[\overline{M}(H,\epsilon ;H_c,\xi )](0)\hspace{1em},\hspace{0.17em}\hspace{0.33em}\hspace{0.33em}\tau (t)=\varphi \hfill \\ \frac{H}{\eta -2Y^M\gamma \epsilon }-\frac{M_R\eta }{\eta -2Y^M\gamma \epsilon }\hspace{1em},\hspace{0.33em}\hspace{0.33em}\tau (t)\ne \varphi \hspace{0.17em}and\hspace{0.17em}H[\max \tau (t)]=-H_c\hfill \\ \frac{H}{\eta -2Y^M\gamma \epsilon }+\frac{M_R\eta }{\eta -2Y^M\gamma \epsilon }\hspace{1em},\hspace{0.33em}\hspace{0.33em}\tau (t)\ne \varphi \hspace{0.17em}and\hspace{0.17em}H[\max \tau (t)]=H_c\hfill \end{array}\right.,$$

(28)

where $H_c=\eta (M_R-M_I)$, $\tau (t)$ can be obtained by Equation (18), and Equation (29) is obtained.

$$[\overline{M}(H,\epsilon ;H_c,\xi )](0)=\left\{\begin{array}{l}\frac{H}{\eta -2Y^M\gamma \epsilon }-\frac{M_R\eta }{\eta -2Y^M\gamma \epsilon }\hspace{1em},\hspace{0.17em}\hspace{0.17em}\hspace{0.33em}H(0)\le -H_c\hfill \\ \xi \hspace{1em},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}-H_c<\hspace{0.33em}H(0)<H_c\hfill \\ \frac{H}{\eta -2Y^M\gamma \epsilon }+\frac{M_R\eta }{\eta -2Y^M\gamma \epsilon }\hspace{1em},\hspace{0.33em}\hspace{0.33em}H(0)\ge H_c\hfill \end{array}\right..$$

(29)

The equilibrium condition of the elastic constitutive relation is expressed by Equation (30).

$$\frac{\partial G}{\partial \epsilon }=0.$$

(30)

Then:

$$\sigma ={\left.\frac{\partial {\psi }_e}{\partial \epsilon }\right|}_M.$$

(31)

For undamped magnetostrictive material, the magnetic–mechanical coupling constitutive relation can be expressed by Equations (32) and (33).

$$\sigma =Y^M\epsilon -Y^M\gamma M^2,$$

(32)

$$M(H,\epsilon )=C{\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }\overline{M}({\rm H}+H,\epsilon ;H_c,\xi )e^{-{{\rm H}}^2/\overline{b}}{e}^{-{(H_c-{\overline{H}}_c)}^2/b}d{\rm H}d{H}_c}}.$$

(33)

When the thermal motion is neglected, $\overline{M}$ can be obtained from Equation (28). When the strong thermal motion or relaxation mechanism is considered, $\overline{M}$ can be derived from Equation (9) and can be obtained by Equation (27).

4. Numerical Implementation of a Magnetization Model Based on the Free Energy Hysteresis Model

According to reference [43], the model is constructed with the piecewise linear kernel function. In order to achieve the highest calculation accuracy and the fastest efficiency, during the numerical solution of the free energy hysteresis model, the Gauss–Legendre algorithm is used to discretize the integral, and the kernel function is realized by using the matrix representation.

4.1. Discretization of Integrals

Divide the integral limit of Equation (25) into a limited number of cells ${\Omega }_2$, and the Gauss–Legendre method is applied to discrete each integration interval. The general form of the Gauss–Legendre algorithm is:

$${\displaystyle {\int }_{-1}^{1}f(x)dx={\displaystyle \sum A_{k}f(x_k)}}.$$

(34)

Then, the numerical equation of the free energy magnetization model can be written by Equation (35).

$$[M(H)](t)={\displaystyle \sum _{i=1}^{N_i}{\displaystyle \sum _{j=1}^{N_j}{\upsilon }_1(H_{ci})}}{\upsilon }_2({{\rm H}}_j)[\overline{M}({{\rm H}}_j+H;H_{ci},{\xi }_i)]v_i{w}_j,$$

(35)

where ${\upsilon }_1(H_{ci})=c_1{e}^{-{(H_{ci}-\overline{H_c})}^2/b}$, ${\upsilon }_2({{\rm H}}_j)=c_2{e}^{-{{\rm H}}_j{}^2/\overline{b}}$, $H_{ci}$ and ${{\rm H}}_j$ are Gaussian integration points, and $v_i$ and $w_j$ are weight functions.

The original discretized equation can be transformed into Equation (36).

$$[M(H)](t)=C{\displaystyle \sum _{i=1}^{N_i}{\displaystyle \sum _{j=1}^{N_j}{e}^{-{(H_{ci}-\overline{H_c})}^2/b}}}{e}^{-{{\rm H}}_j{}^2/\overline{b}}[\overline{M}({{\rm H}}_j+H;H_{ci},{\xi }_i)]v_i{w}_j.$$

(36)

4.2. Kernel Function Implementation

According to Equation (17), the piecewise linear equation of the local magnetization can be described by Equation (37) in the case of where thermal motion is ignored.

$$\overline{M}=\frac{H}{\eta }+M_R\Delta ,$$

(37)

where, when $\Delta =1$, magnetization is located on the upper branch of the kernel function here, and when $\Delta =-1$, magnetization is located on the lower branch of the kernel function here (as shown in Figure 6).

In order to accurately represent the judgment conditions of the Equation (37), and the calculation speed is improved effectively, a matrix algorithm is used, as defined by Equation (38).

$$\begin{array}{cc}\hfill & {\Delta }_{int}={\left[\begin{array}{cc}\begin{array}{c}\begin{array}{ccc}-1& \cdots & -1\end{array}\\ \begin{array}{ccc}⋮& & ⋮\end{array}\\ \begin{array}{ccc}-1& \cdots & -1\end{array}\end{array}& \begin{array}{c}\begin{array}{ccc}1& \cdots & 1\end{array}\\ \begin{array}{ccc}⋮& & ⋮\end{array}\\ \begin{array}{ccc}1& \cdots & 1\end{array}\end{array}\end{array}\right]}_{N_i\times N_j}\hfill \\ \hfill & H_c={\left[\begin{array}{c}\begin{array}{ccc}{H}_{c1}& \cdots & H_{c1}\end{array}\\ \begin{array}{ccc}⋮& & ⋮\end{array}\\ \begin{array}{ccc}{H}_{c{N}_i}& \cdots & H_{c{N}_i}\end{array}\end{array}\right]}_{N_i\times N_j}\hfill \\ \hfill & h_k={\left[\begin{array}{c}\begin{array}{ccc}{H}_k+{}_1& \cdots & H_k+{}_{N_j}\end{array}\\ \begin{array}{ccc}⋮& & ⋮\end{array}\\ \begin{array}{ccc}{H}_k+{}_1& \cdots & H_k+{}_{N_j}\end{array}\end{array}\right]}_{N_i\times N_j}\hfill \end{array}$$

(38)

The weight function vector is defined by Equation (39).

$$\begin{array}{l}{V}^T={\left[\begin{array}{ccc}{v}_1{\upsilon }_1(H_{c1})& \cdots & v_{N_i}{\upsilon }_1(H_{c{N}_i})\end{array}\right]}_{1\times N_i}\\ W^T={\left[\begin{array}{ccc}{w}_1{\upsilon }_2({{\rm H}}_1)& \cdots & w_{N_j}{\upsilon }_1({{\rm H}}_{N_j})\end{array}\right]}_{1\times N_j}\end{array}$$

(39)

According to the magnetic field strength value $H_k$, magnetization $M_k=M(H_k)$, and it can be calculated by the following algorithm:

$$\left\{\begin{array}{l}dH=H_k-H_{k-1}\\ h_k=h_k+dH\\ \Delta =sign(h_k+H_{c\cdot }\ast \Delta )\\ \overline{M}=\frac{h_k}{\eta }+M_R\Delta \\ M_k=V^T\overline{M}W\end{array}\right.,$$

(40)

where $H_k$ is the value of the magnetic field strength corresponding to the $k$ element after discretization, $h_k$ is the effective field, the $ij$ element in the $\Delta$ constructed matrix indicates whether the $j$ effective field just crosses the $i$ critical magnetic field strength, and $H_c$ is the critical magnetic field strength.

4.3. Verification Based on the Free Energy Hysteresis Model

To verify the correctness of the numerical implementation method based on the free energy hysteresis model, the parameters of [41] are used to analyze: $M_R = 3.7 \times {10}^{4}A/m$, $\eta =14$, $H_c=300A/m$, $b = 1 \times {10}^8{A}^2/m^2$, $\overline{b} =\hspace{0.33em}\hspace{0.33em}8 \times {10}^8{A}^2/m^2$, $C = 2.52 \times {10}^{-8}$, and $\gamma = 4.5 \times {10}^{-15}{m}^2/A^2$ (the specific parameters of different giant magnetostrictive actuators can be estimated by the least-squares method). Assuming the pre-pressure is constant, the magnetic anisotropy considered in Equation (36) can be ignored. Then, $\epsilon =0$ in Equation (28), and the measured strain can be solved.

In the period of the input field, it is proposed to take 200 breakpoints for calculation, the relationship between the output magnetization $M$, the output displacement $\sigma$, and the input magnetic field strength $H$ is shown in Figure 6 and Figure 7, respectively. The calculation results of the magnetic machine coupling model in the original reference [41] are shown in Figure 8 and Figure 9. When the two figures are compared, it can be seen that the calculation results are consistent with the simulation results in reference [41].

5. Conclusions

The dynamic magneto-mechanical coupling model of a giant magnetostrictive actuator based on its intrinsic nonlinearity and hysteresis effect is an important theoretical basis for designing GMAs, predicting the displacement output of GMAs, and to provide a reliable basis for GMA system design. The hysteresis characteristics and nonlinearity of giant magnetostrictive actuators have been researched intensively based on the theoretical study of the free energy hysteresis model in this paper. The magnetization model of homogeneous material is derived according to the expression of Helmholtz and Gibbs energy, and two extensions are carried out to obtain the magnetization model of the heterogeneous material. Finally, the magnetic–mechanical coupling constitutive relation of the magnetic field strength $H$ and strain $\epsilon$ is obtained, which is driven from magneto-mechanical coupling effects. Then, the Gauss–Legendre method is used to discretize the integral, the kernel function is realized by using the matrix representation, and the numerical realization of the free energy magnetization model is derived. Reliable experimental data were selected as parameters to verify the free energy magnetization model, and the simulation results are consistent with the existing experimental results. The results show that they are correct for the magnetization model based on the hysteresis characteristics of free energy and its numerical implementation, and the research results provide a theoretical basis for the engineering application of giant magnetostrictive material and the optimization of micro-actuator structures of giant magnetostrictive material.