This paper presents a computationally efficient constitutive model for magnetostrictive materials. High computational efficiency is achieved through the use of local linearization (about easy axes) and discrete energy-averaging techniques. The model is applied to iron-gallium alloys (Galfenol) and tested for different magnetic field orientations relative to the easy axes. It is observed that the model accurately predicts both sensing and actuation characteristics while reducing the computation time by a large factor (1000 times) when compared to the nonlinear energy minimization models. Furthermore, the average error observed in – and – curves is less than 3.5% with the error increasing at magnetic field orientations farther from easy axes, particularly at large magnetic field values. Finally, the model is integrated with a finite element framework to predict the response of a Galfenol rod transducer system, and parametric studies are performed for different current and prestress conditions to optimize the device performance.
I. INTRODUCTION
Magnetostrictive materials can deform in response to an externally applied magnetic field (Joule effect) or can be magnetized by applying stress (Villari effect). Terfenol-D (TbDyFe) is the most commonly used magnetostrictive material that exhibits giant magnetostriction (1600 ppm).1 Galfenol (FeGa), on the other hand, exhibits moderate magnetostriction (450 ppm) but is ductile compared to Terfenol-D and has higher load-carrying capacity.2 Some applications of these materials include energy harvesting,3 torque sensing,4 active vibration control,5–7 and structural health monitoring.8 In order to optimally utilize these materials toward device design, there is a need to develop computationally efficient mathematical models and simulation tools that are capable of predicting their coupled nonlinear response.
Some of the seminal works in ferromagnetic material modeling include the models developed by Armstrong,9 Jiles and Atherton,10 and Preisach.11 Preisach11 developed a hysteresis model for magnetic materials using the mathematical construct of a hysteron which was extended to magnetostrictive materials by Adly and Mayergoyz.12 Armstrong9 calculated bulk magnetization and magnetostriction using an energy-based probability density function and integrating overall possible magnetic moment orientations. To improve computational efficiency, the magnetic moment orientations were subsequently restricted to the easy axes, and the probability density function was used in discrete form.13 Jiles and Atherton developed a hysteresis model incorporating impedances in domain wall motion due to phenomena like domain wall pinning.10 Atulasimha et al.14 developed a constitutive model for Galfenol by summing the contributions over 98 crystallographic directions to predict the material response for magnetic field orientations farther from the easy axes. Zheng and Sun15 developed a phenomenological constitutive model to incorporate prestress and temperature effects in giant magnetostrictive materials. However, most of these models are computationally expensive due to the nonlinearity in the problem and summation of magnetomechanical response over a large number of orientations.
Evans and Dapino16 proposed a discrete energy-averaged (DEA) model for Galfenol by summing over the six easy axes orientations and using a localized form of anisotropy energy defined around each easy axis. This model significantly improved the computational efficiency and showed a close match with experiments for magnetic field orientations close to easy axes. Tari et al.17 extended the model to avoid singularities arising in matrix inversion by using eigenvalue decomposition. Computational efficiency of the DEA model was further improved by Tari et al.18 and Deng et al.19 While the existing DEA models are computationally efficient, they have not been validated for directions farther from easy axes.
Models were also utilized in design and analysis of magnetostrictive material based transducer devices. Datta et al.20 developed a model to predict the quasistatic response of a magnetomechanical cantilever beam for sensing applications. Wun-Fogle et al.21 developed a planar magnetomechanical rotation model to predict sensing performance of an amorphous magnetostrictive material adhered to an aluminum cantilever beam and subjected to various loading conditions. More recently, Shu et al.22 utilized the DEA model to simulate the 1D dynamic response of Galfenol-driven unimorph actuators. The DEA model was also utilized by Santapuri et al.23 to model the dynamic response of composite laminate plate actuators with embedded magnetostrictive materials.
In this work, a computationally efficient yet accurate constitutive model is developed to calculate the magnetostrictive material response for different external magnetic field orientations relative to the easy axes. The model is developed by locally linearizing the total free-energy about each easy axis and subsequently using discrete energy-averaging techniques to obtain the net magnetization and magnetostriction. The model is formulated to avoid any singularities arising at zero magnetic field and stress values without added mathematical complexity. Finally, the model is applied to a Galfenol rod transducer system, and parametric studies are performed to demonstrate its utility toward device design. The modeling framework presented in this work is general and can also be applied to other ferroic materials exhibiting domain microstructure.
This paper is structured as follows: The 3D magnetoelastic governing equations along with different energies in a magnetostrictive material are described in Sec. II. As a first step, a constraint-free anhysteretic model is developed in Sec. III, which is equivalent to the existing energy minimization constitutive models.7 Computational efficiency of the model is significantly improved in Sec. IV by making microstructurally motivated simplifications and locally linearizing the material model presented in Sec. III about each easy axis. The model is applied to single crystal Galfenol in Sec. V to obtain the sensing and actuation characteristics. The results are compared with the existing literature to assess the accuracy and computational efficiency of the model. To further demonstrate the utility of this model toward transducer design, it is integrated with a finite element framework to analyze an axisymmetric rod transducer. Finally, the results are presented with a discussion on the overall efficiency and accuracy of the model.
II. 3D MAGNETOMECHANICAL MODEL
Consider a magnetostrictive medium occupying a volume enclosed by the boundary and surrounded by free-space . The magnetomechanical system is governed by
- Navier’s equation, i.e., linear momentum conservation,where and denote the stress tensor and the external body force, respectively; is the displacement field. The associated boundary condition is(1)where the applied traction is specified on , a subset of the domain boundary ; denotes the outer normal to the domain boundary.The infinitesimal strain is related to the displacement as(2)
- Magnetostatic equations, i.e., Gauss’s law for magnetism and Ampere’s law in the absence of electric fieldwhere is the magnetic field and is the magnetic flux density. The fields are governed by the following jump conditions at the material boundary:(3)where represents jump in the quantity across the boundary .(4)
- Coupled constitutive equations for stress and magnetic flux density in the general formwhere is the fourth-order stiffness tensor; and represent the magnetostriction tensor and magnetization vector functions to be supplied from the material model.(5)
A. Energies of a magnetostrictive material
Energy per unit volume of a magnetizable and deformable medium is given by7
where is the magnetic anisotropy energy, defined for a cubic crystal as
where and are the cubic anisotropy constants and and denote the direction cosines of magnetic moment orientation along the [100], [010], and [001] crystal directions, respectively, such that
Energy due to the applied magnetic field, i.e., the Zeeman energy, is defined as
where , , and are the scalar components of the applied magnetic field along the [100], [010], and [001] crystallographic directions, respectively. is the saturation magnetization (a material constant), and is the permeability of free-space (universal constant).
Finally, the energy due to magnetoelastic interactions within a magnetic material medium is defined as
where and are the magnetostriction constants in the and directions, respectively; is a scaling parameter (see, Atulasimha et al.7); is the applied stress tensor whose components are defined in the Voigt notation as , , , , , and .
III. CONSTITUTIVE MODEL DEVELOPMENT VIA FREE-ENERGY MINIMIZATION
The anhysteretic response of magnetostrictive materials is obtained by minimizing the resulting free-energy (6) constrained by (8). This can be carried out by minimizing the Lagrange function
in terms of variables and and the Lagrange multiplier . However, this technique is computationally intensive as the minimization results in a nonlinear system of equations.
Alternatively, a constraint-free approach may be adopted by parameterizing and in terms of polar angles and as
such that constraint (8) is satisfied for any combination of and values. Equilibrium orientation of magnetic moments can be obtained by minimizing the free-energy (6) with respect to and , i.e., we solve for and satisfying
such that
To avoid redundancy, only solutions in the range and are considered. In the absence of any external magnetic field or prestress, these energy-minimizing orientations correspond to crystal easy axes.
This model is equivalent to the existing models by Armstrong13 and Atulasimha et al.,14 and the results obtained match closely with experimental data. Although better than the minimization of the Lagrange function (11), this methodology is also computationally expensive due to the nonlinearity of the system of Eq. (13). This model is thus used as a benchmark to compare the results obtained from the locally linearized model, to be presented in Sec. IV, to verify its accuracy and computational efficiency.
In what follows, microstructurally motivated simplifications will be made to this model to further improve computational efficiency while maintaining accuracy.
IV. COMPUTATIONALLY EFFICIENT LOCALLY LINEARIZED CONSTITUTIVE MODEL DEVELOPMENT
Noting that magnetic moment orientation remains uniform within each domain, we develop localized expressions for energy in each domain and subsequently average over the total volume to obtain net magnetization and magnetostriction. This approach, also referred to as “discrete energy-averaging,” has been used by Armstrong13 and Evans and Dapino.16 To obtain the energy-minimizing magnetic moment orientations, the following simplifying assumptions are made:
Contributions of magnetic exchange energy and magnetostatic energy are neglected. This assumption is valid for bulk material systems.
Magnetic moment orientations remain close to crystal easy axes. This assumption allows us to develop approximate expressions for energy in the neighborhood of easy axes.
Given a magnetic material with mutually independent easy axes directions, let denote the unit vector along the th easy axis (). Owing to assumption (ii) above, a localized form of the free-energy function (6) in the th domain (corresponding to the th easy axis) is defined as
where in the standard index notation for summation over the variables and is used. denotes the free-energy function (6) evaluated at . This localized expression is valid in a neighborhood of , wherein the equilibrium configuration of the th domain can be calculated by minimizing this localized expression subjected to constraint (8).
However, nonlinearity of constraint (8) results in a nonlinear 6th order algebraic equation in terms of the Lagrange multiplier . Thus, we continue to utilize our assumption (ii) to linearize the constraint in the vicinity of to obtain
The constraint is thus reduced to
upon linearization. The reduced Lagrangian function is thus defined as
in which the minimization is performed by solving the following system of equations:
along with the linearized constraint rewritten as
wherein the property is utilized. Substituting the localized free-energy expression (15) into (18), we obtain
where
and
where
Solution (21) is valid for invertible such that . This requirement may be violated for certain choices of and stress values. In the case of Galfenol, for instance, becomes noninvertible when since easy axes are along . This is dealt with by replacing in Eq. (18) with
which is invertible in the case of Galfenol. It is noted that
due to constraint (19).
As the domain orientations move farther from easy axes , the reduced model may not be able to describe the material behavior accurately due to the linearization of the constraint. In such cases, the magnetic moment is renormalized as
to ensure that the magnitude of magnetic orientation does not become greater than unity.24
A. Calculation of bulk magnetostriction, magnetization, and magnetic flux density
Bulk magnetization can be calculated by summing magnetic moment orientations over all domains
where represents saturation magnetization, is the number of easy axes, and is the magnetic moment orientation of the th domain obtained using (24). Also, is the anhysteretic volume fraction of the th domain defined as the ratio of domain volume to sample volume and can be described using Boltzmann distribution23 as
The magnetic flux density can now be calculated using
Also, energy per unit volume of the magnetostrictive material can be calculated by summing magnetic anisotropy, magnetomechanical energy, and Zeeman energy across all domains, i.e.,
Similarly, the bulk magnetostriction can be calculated using
where is the magnetostriction (strain due to magnetization) of the th domain and is dependent on . Particularly, for a cubic system,25,26
If the material is sufficiently prestressed such that all magnetic moments are perpendicular to the direction of magnetization at the beginning of the magnetization process, it is common practice to omit the negative one-third term appearing in the normal strain components , and .27
In the presence of external stresses (within the elastic regime), the total strain can be calculated as
where represents the elastic compliance matrix for the magnetostrictive material.
B. Incorporation of hysteresis
The constitutive model discussed thus far does not incorporate hysteretic effects arising due to material defects like pinning sites. To incorporate hysteresis, an incremental model adopted by Armstrong13 and Jiles et al.10 is utilized in this paper. The evolution equation for a particular domain can be expressed as (overset ignored for notational convenience)
where is a material constant quantifying the number of pinning sites per unit volume, is the magnitude of the applied magnetic field, is the anhysteretic volume fraction calculated using (26), and is the hysteretic volume fraction.
A stepwise finite-difference scheme is adopted to calculate the hysteretic volume fraction. Starting from a known initial volume fraction at a specified , we apply incremental field and calculate , in a stepwise manner as shown below,
where the subscript denotes the th iteration. The bulk magnetoelastic properties are now calculated using instead of in (25)–(29).
V. RESULTS AND DISCUSSION
The locally linearized constitutive model outlined in Sec. IV is applied to Galfenol. The results are compared against the nonlinear model discussed in Sec. III (which is equivalent to the model presented by Atulasimha et al.28). The material parameters used here are taken from the literature28 and are listed in Table I. Easy axes of Galfenol are along the crystallographic directions , i.e., the number of easy axes is given by . The constitutive model response is categorized into (i) actuation characteristics, i.e., – and – curves for different prestress values and (ii) sensing characteristics, i.e., – and – curves for different magnetic field values. The model is subsequently applied to a Galfenol rod actuator system, and the magnetomechanical response of the actuator is analyzed.
Parameters (unit) . | Name . | Value . |
---|---|---|
(–) | Magnetostrictive constant | |
(–) | Magnetostrictive constant | |
(H m) | Vacuum permeability | |
(A m) | Saturation magnetization | |
, (J m) | Anisotropy coefficient | |
(GPa) | Young’s modulus | 59 |
(kg m) | Density | 7870 |
(J m) | Smoothing constant | 625 |
(–) | Fit constant | 0.8 |
Parameters (unit) . | Name . | Value . |
---|---|---|
(–) | Magnetostrictive constant | |
(–) | Magnetostrictive constant | |
(H m) | Vacuum permeability | |
(A m) | Saturation magnetization | |
, (J m) | Anisotropy coefficient | |
(GPa) | Young’s modulus | 59 |
(kg m) | Density | 7870 |
(J m) | Smoothing constant | 625 |
(–) | Fit constant | 0.8 |
A. Galfenol actuation and sensing characteristics
1. Actuation characteristics
The – and – characteristics of single crystal Galfenol are presented in Figs. 1 and 2 for different applied magnetic field orientations relative to the easy axes. Results generated using our locally linearized constitutive model are compared to the existing nonlinear model,28 which is experimentally validated. The present model shows a close match with existing results while improving the computational efficiency. A slower saturation is observed as the angle between magnetic field direction and easy axes increases. Furthermore, it is observed that the accuracy of the present model reduces for magnetic field directions farther from the easy axis, especially for larger applied fields (). This is to be expected as our model assumes magnetic moments within each domain to remain close to the easy axes, and this assumption does not hold well as the angle between them increases.
The magnetization process in any magnetic material can be attributed to the following two phenomena: (i) domain growth (for small/medium external magnetic fields), wherein the domains aligned favorably to the field direction grow at the expense of domains that are aligned farther from the field; (ii) domain rotation (for large external magnetic fields), wherein the favorable domains that have survived start rotating toward the field direction and are fully aligned when the applied magnetic field is sufficiently large.10 However, in our model, magnetic moments are assumed to remain close to easy axes, i.e., the effect of domain rotation is not incorporated, resulting in reduced accuracy for larger field values [see Figs. 1(c) and 2(a), for instance].
2. Sensing characteristics
The sensing characteristics, i.e., the – and – characteristics of Galfenol single crystal with prestress along different directions are presented in Figs. 3 and 4. The nominal Young’s modulus value of 59 GPa along the [100] direction is used.7 Due to the magnetostriction of Galfenol, the slope of strain vs stress curve, which quantifies the material compliance, varies with both applied magnetic field and prestress values, as shown in Fig. 4.
It can be observed from Fig. 4 that the stress–strain behavior within elastic limit is linear at very low magnetic fields but becomes nonlinear for moderate to large values. However, for , the – curve becomes linear again. In other words, for large magnetic field inputs, even large compressive stresses do not produce a change in the orientation of the magnetic moments, i.e., the strain produced here is purely elastic.
Slope of the – curve is known as sensitivity. It can be observed from Fig. 3 that sensitivity is much larger along [100] compared to the other two cases.
3. Error analysis and computational efficiency
To assess the accuracy of our framework, the locally linearized constitutive model is compared with the constraint-free nonlinear model (presented in Sec. III). Normalized root mean square percentage error of both – and – curves are calculated using
where and represent the magnetostriction values calculated using the nonlinear model and the locally linearized model, respectively, and denotes the total number of data points.
The error percentage for the – curves is observed to be within 2%, and the error in – curves is within 3.5% for all magnetic field orientations as shown in Fig. 5. The error in – curve is seen to be higher for directions farther from easy axis. The error is predominately higher for moderate/large magnetic field values (). At 400 Oe, maximum percentage relative error calculated is along [211] for both – and – curves with a magnitude of 2.4% and 9.7%, respectively.
However, the model presents a major computational advantage compared to the fully nonlinear minimization problem. The computation time elapsed for plotting – and – curves for a specified prestress (along [100]) using the locally linearized model is 0.239 077 s, whereas the nonlinear constraint-free model takes around 37 min (on a Windows 7 desktop, Intel Xeon 3.4 GHz CPU, and 32 GB RAM).
Thus, our locally linearized model has a reduced computation time by more than 1000 times while maintaining accuracy for most magnetic field orientations. The present model is formulated to avoid any singularities arising at zero magnetic field and stress values without added mathematical complexity. Furthermore, the model utilizes the material parameters already reported in the literature and does not require any parameter optimization. Owing to its overall efficiency, it is very convenient to integrate this model with a finite element framework to perform transducer design and analysis.
4. Hysteresis effects
The hysteresis model needs one additional parameter, , which is proportional to the density of pinning sites. This value has been reported for Galfenol as .14 Hysteresis plots are generated along [100] by combining the locally linearized model (21) with incremental domain evolution equation (33) to incorporate losses due to pinning sites. The hysteretic – and – curves thus obtained are shown in Fig. 6 and compared with the results obtained from the nonlinear model.
In what follows, utility of the locally linearized model is demonstrated by applying it toward the analysis of a Galfenol rod transducer system.
B. Analysis of a Galfenol rod transducer
In this section, the locally linearized constitutive model is utilized to analyze the response of a Galfenol rod under different input magnetic field and prestress conditions. Response of the transducer is analyzed for both sensing and actuation modes.
The transducer setup described in Fig. 7 consists of a cylindrical Galfenol rod fixed at the bottom, a steel casing to provide a flux return path, and the surrounding air domain (modeled as free-space) sufficiently large to ensure zero magnetic potential at its outermost boundary.
The Galfenol rod needs to be magnetically excited for both actuation and sensing applications. This can be achieved using a current carrying coil axisymmetrically wound around the magnetostrictive rod. Alternatively, the external field (or background field) may be directly specified at the outer free-space boundary as shown in Fig. 7. The presence of a magnetic medium (such as Galfenol) results in perturbation of the applied field, which can be calculated by solving the magnetomechanical governing equations (1)–(5). Axisymmetry in the problem is utilized ( is the axis of symmetry) to further improve computational efficiency.
In actuator mode, the deflection of the magnetostrictive rod is calculated as the magnetic field induced in the rod is varied (by varying the coil current) for a fixed prestress value. In the sensor mode, magnetic induction in the rod is studied by varying the input axial force at a fixed external magnetic field.
1. Solution methodology
The magnetomechanical transducer system is analyzed using governing equations (1)–(5). Finite element analysis is used to solve the resulting system of equations along with the locally linearized constitutive model. The weak form of the governing equations are discretized in COMSOL Multiphysics (version 5.3a), and the bulk magnetization and magnetostriction , obtained from (25) and (29), are supplied to COMSOL as external MATLAB functions.
2. Simulation results
1. Actuator analysis
A Galfenol rod of length 50 mm and diameter 6 mm is actuated using the field produced by a current carrying coil. Parametric simulations are performed to obtain magnetostriction vs coil current density (–) and magnetic flux density vs coil current density (–) plots at various prestress values as illustrated in Fig. 8. The maximum saturation magnetostriction is observed to be 170.2 ppm for zero prestress and 255 ppm for all other prestress values. Saturation is observed at large current density values, i.e., 1500 kA/m. The 2D axisymmetric element mesh used here consists of 265 domain elements and 93 boundary elements.
2. Sensor analysis
In sensor mode, the magnetic flux density induced in the magnetostrictive rod is studied by varying applied axial load . No variation in magnetic induction is observed for zero applied magnetic field. Thus, an external field is specified in the form on the outer boundary as shown in Fig. 7, where is the unit vector along the rod’s axis. Parametric simulations are performed on a rod of length 100 mm and diameter 10 mm to obtain – characteristics at various external field values as illustrated in Fig. 9(b). Also, the variation of magnetic flux density along the axis of symmetry is studied in Fig. 9(a). It is observed that the magnetic flux density is uniform in the free-space but reaches a maximum at the center of the rod. Furthermore, it is observed that the magnetic flux density decreases with increasing compressive load. Saturation is observed for large external fields, i.e., 877.5 Oe, beyond which force sensing will not be feasible. The 2D axisymmetric element mesh used in these simulations consists of 108 domain elements and 40 boundary elements.
The computation time of the COMSOL model with integrated locally linearized model for one simulation is in the range of 1–5 min. The fully nonlinear model takes at least a few hours for each simulation and has greater convergence issues (all simulations were performed on a Windows 7 desktop with Intel Xeon 3.4 GHz CPU and 32 GB RAM).
VI. CONCLUSIONS
A computationally efficient locally linearized model has been developed to predict the coupled magnetomechanical behavior of magnetostrictive materials. This model accurately predicts the sensing and actuation characteristic while improving computational efficiency by more than a factor of 1000. While the model describes the – and – curves accurately (average error 3.5%) for magnetic field orientations close to the easy axes , the error margin increases as we move farther away from the easy axis, particularly for . Another advantage of this model is that no further parameter optimization is needed for this model as it utilizes the standard magnetostrictive material constants reported in the literature. Furthermore, the model is formulated to avoid singularities at zero field values without added mathematical complexity.
To demonstrate the utility of this model in transducer design, our locally linearized model was integrated with the finite element framework to analyze an axisymmetric rod system. As expected, the finite element implementation using the locally linearized model demonstrated a major computational advantage compared to the well established nonlinear models while maintaining accuracy. Due to its relatively simple execution and significant computational advantage, this framework can be effectively used as a design tool.
ACKNOWLEDGMENTS
S. Santapuri, S. K. Wahi, and M. Kumar gratefully acknowledge the financial support received from Science and Engineering Research Board, DST (File No. ECR/2017/002161) and IIT Delhi Seed Grant. M. J. Dapino would like to acknowledge financial support by the member organizations of the Smart Vehicle Concepts Center, a National Science Foundation Industry-University Cooperative Research Center (www.SmartVehicleCenter.org) established under NSF Grant Nos. IIP-1238286 and IIP-1738723.