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 $\lambda $–$H$ and $B$–$H$ 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 (Tb$0.3$Dy$0.7$Fe$1.9$) is the most commonly used magnetostrictive material that exhibits giant magnetostriction ($\u2248$1600 ppm).^{1} Galfenol (Fe$81.6$Ga$18.4$), on the other hand, exhibits moderate magnetostriction ($\u2248$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} Preisach^{11} 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} Armstrong^{9} 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 Sun^{15} 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 Dapino^{16} 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 $V$ enclosed by the boundary $\u2202V$ and surrounded by free-space $V\u2217$. The magnetomechanical system is governed by

- Navier’s equation, i.e., linear momentum conservation,where $T$ and $fB$ denote the stress tensor and the external body force, respectively; $u$ is the displacement field. The associated boundary condition is(1)$divT+fB=\rho \u22022u\u2202t2inV,$where the applied traction $ta$ is specified on $\u2202Vt$, a subset of the domain boundary $\u2202V$; $n$ denotes the outer normal to the domain boundary.$ta=Tnon\u2202Vt,$The infinitesimal strain $S$ is related to the displacement $u$ as(2)$S=12(gradu+graduT).$
- Magnetostatic equations, i.e., Gauss’s law for magnetism and Ampere’s law in the absence of electric fieldwhere $H$ is the magnetic field and $B$ is the magnetic flux density. The fields are governed by the following jump conditions at the material boundary:(3)$curlH=0,divB=0inV\u222aV\u2217,$where $[.]$ represents jump in the quantity across the boundary $\u2202V$.(4)$n\xd7[H]=0,n\u22c5[B]=0on\u2202V,$
- Coupled constitutive equations for stress and magnetic flux density in the general formwhere $C$ is the fourth-order stiffness tensor; $\lambda (T,H)$ and $M(T,H)$ represent the magnetostriction tensor and magnetization vector functions to be supplied from the material model.(5)$T=C(S\u2212\lambda (T,H)),B=\mu o(H+M(T,H))inV,$

### A. Energies of a magnetostrictive material

Energy per unit volume of a magnetizable and deformable medium is given by^{7}

where $Uanisotropy$ is the magnetic anisotropy energy, defined for a cubic crystal as

where $K1$ and $K2$ are the cubic anisotropy constants and $m1,m2,$ and $m3$ denote the direction cosines of magnetic moment orientation $m$ 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 $H1$, $H2$, and $H3$ are the scalar components of the applied magnetic field $Ha$ along the [100], [010], and [001] crystallographic directions, respectively. $Ms$ is the saturation magnetization (a material constant), and $\mu o$ is the permeability of free-space (universal constant).

Finally, the energy due to magnetoelastic interactions within a magnetic material medium is defined as

where $\lambda 100$ and $\lambda 111$ are the magnetostriction constants in the $\u27e8100\u27e9$ and $\u27e8111\u27e9$ directions, respectively; $\gamma \sigma $ is a scaling parameter (see, Atulasimha *et al.*^{7}); $T$ is the applied stress tensor whose components are defined in the Voigt notation as $T1=T11$, $T2=T22$, $T3=T33$, $T4=T12$, $T5=T23$, and $T6=T13$.

## 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 $m1,m2,$ and $m3$ and the Lagrange multiplier $\mu $. 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 $m1,m2,$ and $m3$ in terms of polar angles $\theta $ and $\varphi $ as

such that constraint (8) is satisfied for any combination of $\varphi $ and $\theta $ values. Equilibrium orientation of magnetic moments can be obtained by minimizing the free-energy (6) with respect to $\varphi $ and $\theta $, i.e., we solve for $\theta o$ and $\varphi o$ satisfying

such that

To avoid redundancy, only solutions in the range $0\u2264\varphi o<2\pi $ and $0\u2264\theta o<\pi $ 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 Armstrong^{13} 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 Armstrong^{13} 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 $n$ mutually independent easy axes directions, let $c\alpha $ denote the unit vector along the $\alpha $th easy axis ($\alpha =1,2,\u2026,n$). Owing to assumption (ii) above, a localized form of the free-energy function (6) in the $\alpha $th domain (corresponding to the $\alpha $th easy axis) is defined as

where in the standard index notation for summation over the variables $i$ and $j$ is used. $U|c\alpha $ denotes the free-energy function (6) evaluated at $m=c\alpha $. This localized expression is valid in a neighborhood of $c\alpha $, wherein the equilibrium configuration of the $\alpha $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 $\mu $. Thus, we continue to utilize our assumption (ii) to linearize the constraint $f(m1,m2,m3)=m12+m22+m32\u22121$ in the vicinity of $c\alpha $ 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 $c1\alpha 2+c2\alpha 2+c3\alpha 2=1$ is utilized. Substituting the localized free-energy expression (15) into (18), we obtain

where

and

The localized magnetic moment orientations $m\alpha $ can be solved using the system of Eqs. (19) and (20) as

where

Solution (21) is valid for invertible $\kappa \alpha $ such that $\kappa \alpha \u22121c\alpha \u22c5c\alpha \u22600$. This requirement may be violated for certain choices of $c\alpha $ and stress values. In the case of Galfenol, for instance, $\kappa \alpha $ becomes noninvertible when $T=0$ since easy axes $c\alpha $ are along $\u27e8100\u27e9$. This is dealt with by replacing $\kappa \alpha $ in Eq. (18) with

which is invertible in the case of Galfenol. It is noted that

due to constraint (19).

As the domain orientations $m\alpha $ move farther from easy axes $c\alpha $, 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 $m\alpha $ 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 $M$ can be calculated by summing magnetic moment orientations over all domains

where $Ms$ represents saturation magnetization, $n$ is the number of easy axes, and $m~\alpha $ is the magnetic moment orientation of the $\alpha $th domain obtained using (24). Also, $\xi an\alpha $ is the anhysteretic volume fraction of the $\alpha $th domain defined as the ratio of domain volume to sample volume and can be described using Boltzmann distribution^{23} as

The magnetic flux density $B$ 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 $\lambda m$ can be calculated using

where $\lambda m\alpha $ is the magnetostriction (strain due to magnetization) of the $\alpha $th domain and is dependent on $m~\alpha $. 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 $\lambda 1,\lambda 2$, and $\lambda 3$.^{27}

In the presence of external stresses (within the elastic regime), the total strain can be calculated as

where $s$ 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 Armstrong^{13} and Jiles *et al.*^{10} is utilized in this paper. The evolution equation for a particular domain can be expressed as (overset $\alpha $ ignored for notational convenience)

where $kp$ is a material constant quantifying the number of pinning sites per unit volume, $H$ is the magnitude of the applied magnetic field, $\xi an$ is the anhysteretic volume fraction calculated using (26), and $\xi hyst$ 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 $H$, we apply incremental field $\Delta H$ and calculate $\xi hyst$, in a stepwise manner as shown below,

where the subscript $k$ denotes the $k$th iteration. The bulk magnetoelastic properties are now calculated using $\xi hyst$ instead of $\xi an$ 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 literature^{28} and are listed in Table I. Easy axes of Galfenol are along the crystallographic directions $\u27e8100\u27e9$, i.e., the number of easy axes is given by $n=6$. The constitutive model response is categorized into (i) actuation characteristics, i.e., $\lambda $–$H$ and $B$–$H$ curves for different prestress values and (ii) sensing characteristics, i.e., $S$–$T$ and $B$–$T$ 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 . |
---|---|---|

$(3/2)\lambda 100$(–) | Magnetostrictive constant | $255\xd710\u22126$ |

$(3/2)\lambda 111$(–) | Magnetostrictive constant | $\u22127\xd710\u22126$ |

$\mu 0$ (H m$\u22121$) | Vacuum permeability | $4\pi \xd710\u22127$ |

$Ms$ (A m$\u22121$) | Saturation magnetization | $1.83/\mu 0$ |

$K1$,$K2$ (J m$\u22123$) | Anisotropy coefficient | $3.6\xd7104,0$ |

$E$ (GPa) | Young’s modulus | 59 |

$\rho $ (kg m$\u22123$) | Density | 7870 |

$\Omega $ (J m$\u22123$) | Smoothing constant | 625 |

$\gamma \sigma $(–) | Fit constant | 0.8 |

Parameters (unit) . | Name . | Value . |
---|---|---|

$(3/2)\lambda 100$(–) | Magnetostrictive constant | $255\xd710\u22126$ |

$(3/2)\lambda 111$(–) | Magnetostrictive constant | $\u22127\xd710\u22126$ |

$\mu 0$ (H m$\u22121$) | Vacuum permeability | $4\pi \xd710\u22127$ |

$Ms$ (A m$\u22121$) | Saturation magnetization | $1.83/\mu 0$ |

$K1$,$K2$ (J m$\u22123$) | Anisotropy coefficient | $3.6\xd7104,0$ |

$E$ (GPa) | Young’s modulus | 59 |

$\rho $ (kg m$\u22123$) | Density | 7870 |

$\Omega $ (J m$\u22123$) | Smoothing constant | 625 |

$\gamma \sigma $(–) | Fit constant | 0.8 |

### A. Galfenol actuation and sensing characteristics

#### 1. Actuation characteristics

The $\lambda $–$H$ and $B$–$H$ 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 ($H>Hsat$). 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 $S$–$T$ and $B$–$T$ 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 $H>Hsat$, the $S$–$T$ 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 $B$–$T$ 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 $B$–$H$ and $\lambda $–$H$ curves are calculated using

where $\lambda nonlineari$ and $\lambda locali$ represent the magnetostriction values calculated using the nonlinear model and the locally linearized model, respectively, and $N$ denotes the total number of data points.

The error percentage for the $B$–$H$ curves is observed to be within 2%, and the error in $\lambda $–$H$ curves is within 3.5% for all magnetic field orientations as shown in Fig. 5. The error in $\lambda $–$H$ curve is seen to be higher for directions farther from easy axis. The error is predominately higher for moderate/large magnetic field values ($H>Hsat$). At $H=$ 400 Oe, maximum percentage relative error calculated is along [211] for both $B$–$H$ and $\lambda $–$H$ 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 $B$–$H$ and $\lambda $–$H$ 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, $kp$, which is proportional to the density of pinning sites. This value has been reported for Galfenol as $kp=150$.^{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 $B$–$H$ and $\lambda $–$H$ 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 $Ha$ (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 ($AA\u2032$ 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 $M(T,H)$ and magnetostriction $\lambda (T,H)$, 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 ($\lambda $–$J0$) and magnetic flux density vs coil current density ($B$–$J0$) 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., $J0\u2265$ 1500 kA/m$2$. 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 $Fz$. No variation in magnetic induction is observed for zero applied magnetic field. Thus, an external field is specified in the form $Ha=Hak$ on the outer boundary as shown in Fig. 7, where $k$ 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 $B$–$Fz$ characteristics at various external field values as illustrated in Fig. 9(b). Also, the variation of magnetic flux density along the axis of symmetry $AA\u2032$ 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., $Ha$ $\u2265$ 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 $\lambda $–$H$ and $B$–$H$ curves accurately (average error $<$3.5%) for magnetic field orientations close to the easy axes $\u27e8100\u27e9$, the error margin increases as we move farther away from the easy axis, particularly for $H>Hsat$. 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.

## REFERENCES

*Smart Sensor Phenomena, Technology, Networks, and Systems 2008*(International Society for Optics and Photonics, 2008), Vol. 6933, p. 693314.

*Sensor Systems and Networks: Phenomena, Technology, and Applications for NDE and Health Monitoring 2007*(International Society for Optics and Photonics, 2007), Vol. 6530, p. 65300C.