Micromagnetic modeling of magnetostrictive materials
under intrinsic stress
Y.C. Shu
*, M.P. Lin, K.C. Wu
Institute of Applied Mechanics, National Taiwan University, Taipei 106, Taiwan, ROC Received 16 September 2002; received in revised form 25 March 2003
Abstract
We have developed a framework based on micromagnetics to explore the effect of stress on the magnetostrictive behavior in ferromagnetics. Our approach is different from the conventional one which simply replaces the total strain by magnetostrain. Question arises for such an approach because of the loss of strain compatibility. Here, we have included the kinematic constraints in our micromagnetic model and developed a modified boundary integral formalism to calculate the intrinsic stress induced by incompatible magnetostrain. We have shown that for small magnetostriction of the order of 105, the results predicted by the present approach are slightly different from those predicted by the
conventional method. But we have found that for large magnetostriction around 103 order of magnitude, the
con-ventional approach is insufficient to predict magnetic domain patterns and hysteresis precisely, and the effective magnetic field induced by intrinsic stress cannot be neglected.
2003 Elsevier Ltd. All rights reserved.
Keywords: Incompatible magnetostrain; Micromagnetics; Magnetostriction
1. Introduction
A magnetostrictive material is one which is able to change its physical dimensions in response to a change of state of magnetization. In other words, it exhibits a change in length accompanied by an inverse change in girth when it is subjected to a magnetic field. Conversely, the state of magneti-zation changes if an external force is applied causing strain. The coupling between magnetic and mechanical energies gives rise to the capability
in transduction that allows a magnetostrictive material to be used in actuator and sensor appli-cations. All ferromagnetic materials exhibit vary-ing degrees of magnetostrictive strain, however only some of them exhibit sufficient strain for practical use. The highest room temperature magnetostrictive strain of a pure element is that of Co which is extremely small and is around 60 microstrain at saturation. Fe and Ni are common ferromagnetic materials and also have magneto-striction with magnitudes around this range. On the other hand, by alloying elements ‘‘giant’’ magnetostriction under relatively small magnetic fields can be achieved (Clark, 1980). Among these, Terfenol-D (TbxDy1 xFe2, x 0:3) is at present
the most widely used magnetostrictive material. It
*
Corresponding author. Tel.: 2-3366-5627; fax: +886-2-2363-9290.
E-mail address:[email protected](Y.C. Shu). 0167-6636/$ - see front matter 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechmat.2003.04.004
is capable of strains as high as 1600 microstrain and, since the 1980s, has enjoyed the greatest commercial success for application in a great many fields (Jiles, 1994).
Magnetostrictive materials at no applied fields often display extremely complex domains on which magnetization and strain are approximately con-stant. These domains, occurring at a microscopic scale, evolve under the action of applied fields and loads which in turn give rise to the overall response of the materials. Therefore, studying the mecha-nism that governs domain evolution is a key step toward understanding the macroscopic behavior of magnetostrictive materials. In the past several decades domain analysis and evolution under external fields have been studied extensively for rigid ferromagnetic materials within the framework of micromagnetic theory developed by Brown (1963) and his followers (see Aharoni, 1996 for an extensive review and James and Kinderlehrer, 1990 for a modern treatment). In spite of many useful results predicted by this theory, very few analytic solutions of BrownÕs micromagnetic equations have been found in literature as these equations are inherently nonlinear and nonlocal in nature. Numerical micromagnetics, on the other hand, provides an alternative to apply micromagnetic methods to large-scale domain structures. Nowa-days computational micromagnetics has led to a deeper understanding of hysteresis behavior by visualization of magnetization reversal (Fidler and Schrefl, 2000).
Complicated magnetomechanical coupling and time consuming computation, however, have pre-vented many useful methods proposed originally for rigid ferromagnetics from being applied di-rectly to deformable ferromagnetics. Much has been done for theoretical development of magne-tostrictive materials. It includes early pioneering works such as Brown (1966), Tiersten (1965), Pao and Yeh (1973), Hutter and van de Ven (1978) and Jiles and Atherton (1984) (see also Eringen and Maugin, 1990; Maugin, 1988 and references therein). James and Kinderlehrer (1993) and DeSimone and James (2002) have proposed an energetic approach and used very sophisticated modern mathematical techniques to analyze do-mains in giant magnetostrictive materials.
How-ever, little or none has provided a convenient way for direct simulation of domain patterns and evolution under external fields and loads for a wide class of magnetostrictive materials. Fabian and Heider (1996) have used the continuum theory of defects to include magnetostriction in micro-magnetic calculation with some success in study-ing titanomagnetites. Recently Voltairas et al. (2000a,b) have studied nonuniform magnetization reversal in stressed ferromagnetic films. Their simulation of magnetization reversal, however, is confined to be one-dimensional with severe restrictions on elastic states.
To overcome difficulties in the calculation of magnetostrictive energy, the conventional ap-proach to include magnetostriction is to replace the total compatible strain by magnetostrictive strain in micromagnetic models (Cullity, 1972; Zhu, 1993; Callegaro and Puppin, 1997). As a consequence, the energy of applied loads becomes another source of magnetic anisotropy energy while the magnetostrictive self-energy vanishes. The intrinsic stress resulting from incompatible magnetostrictive strain as magnetization rotates is not taken into account while the external stress causes the induction of a magnetostrictive uniaxial anisotropy, usually called stress-induced anisotropy (Izawa, 1984; Jeong and Bogy, 1995). The material is then treated as a rigid ferromagnetics with varying magnetic easy axes depending on external stress. Recently Zhu et al. (2001) have used this approach to study the influence of external stress on coercivity in magnetic thin films and obtained qualitatively general agreement with experiment (Callegaro and Puppin, 1996). Although the re-sults obtained by this approximation have been found to be useful in various applications, ques-tion arises because of loss of strain compatibility and negligence of intrinsic stress. This paper is primarily concerned with this issue and studies the effect of intrinsic stress on the magnetostrictive behavior in various ferromagnetic materials. In addition, a new method for the direct simulation of magnetic domains and their evolution under external fields and loads is proposed here.
Our framework is based on micromagnetics accounting for magnetostrictive effect. Micro-magnetics is a nonlocal, nonconvex variational
problem. We write in Section 2 the total energy of the crystal as a sum of the exchange energy be-tween the spins, the total anisotropy energies including magnetostrictive energy, the energy of the applied fields and loads, and the energy of the induced magnetic field. The state of the ferro-magnetic crystal is obtained as minimizers of this total energy subject to two constrained equations: MaxwellÕs and elastic equilibrium equations. The anisotropy energy density encodes phenomeno-logically the information that the crystal prefers certain magnetized and distorted states, and can be nonconvex with multiple wells due to the presence of multiple such states. It can be shown that energy minimization based on this theory automatically leads to domain patterns in equilibrium based on very few parameters (James and Kinderlehrer, 1993; James and Wuttig, 1998). Instead of pro-viding further mathematical analysis on magnetic domains, we study the solutions of micromagnetic equations deduced from the present theory. We use the Landau–Lifshitz–Gilbert (LLG) equation to obtain both equilibrium and transient magne-tization configurations. The LLG equation is de-rived based on the microscopic physics and is able to describe the relaxation process of the magneti-zation. We emphasize this point by showing in Section 2.4 that the total free energy is decreasing as magnetization evolves under the LLG equation. The key issue in performing simulation of do-main evolution is to solve the elastic equilibrium equation for each magnetization configuration. This method, which we call the constrained ap-proach here, is in contrast to the relaxed apap-proach that neglects the effect of intrinsic stress. In Section 3, we calculate the magnetostrictive energy which requires the solution of a boundary value problem in elasticity theory. Since such a calculation has to be performed at every step of the iterative proce-dure, a time consuming computation seems to be inevitable. Fortunately, the improved availability of large-scale computer power has enabled us to include magnetostrictive effect in the LLG solver and provides us an opportunity for deeper understanding of detailed and subtle physical behavior of magnetization reversal. As the stress-induced magnetic field (see (2.14)) involves vari-able strain rather than displacement, we develop a
strain-based boundary integral formulation in Section 3.3. This formalism, originally proposed by Wu et al. (1992), has been shown to be more accurate in the calculation of stress as well as have faster rate of convergence than conventional dis-placement-based numerical schemes. We extend their formulation to account for the effect of magnetostrictive stress and implement the associ-ated numerical algorithm in Section 4.
The simulation results for various magneto-strictive films are presented in Section 5. We study magnetization in equilibrium and reversal using constrained and relaxed methods and make com-parison between these two approaches. We use Ni and Terfenol-D as the model materials to repre-sent for ferromagnetics with small and large magnetostrictive strains, respectively. Our simula-tion demonstrates that the relaxed approach turns out to be a reasonable approximation in predicting equilibrium domains and their evolution in Ni films while it becomes a crude estimation for Terfenol-D films. It shows that the effective mag-netic field induced by intrinsic stress cannot be ignored for materials with large magnetostrictive strain. We conclude in Section 6 with a discussion.
2. Framework and formulation 2.1. Kinematics and magnetostriction
Magnetostriction is the spontaneous deforma-tion of a magnetic crystal caused by a change of its state of magnetization. Certain materials such as Ni2MnGa, one of ferromagnetic shape-memory
materials, have extremely large magnetostriction of the order of 101(Ullakko et al., 1996; Tickle and
James, 1999; Wu et al., 2002). They are not taken into account here because of difficulty in handling physical quantities defined in deformed configura-tion. We refer to James and Kinderlehrer (1993) and James (2002) for detailed discussion. Instead, materials under consideration here such as Fe, Ni or Terfenol-D have magnetostriction with orders ranging from 103 to 105. Such a small change in
shape suggests us to choose the geometrically linear framework of micromagnetics, which Brown has called the conventional theory of magnetostriction
(Brown, 1966). However, we are aware that, as pointed out by Brown (1966), this conventional linear strain formulation misses certain terms that should be present according to a direct lineariza-tion starting from the geometrically nonlinear theory. As a result, the magnetic body force and body couple are neglected, and the anti-symmetric part of stress vanishes. With these reservations, we consider a ferromagnetic single crystal shown in Fig. 1. The crystal occupies the region X R3in its reference configuration. The displacement of the crystal is described by the function u : X! R3
and the linear strain is denoted by e : X! M33
s where
Mmms is the set of all symmetric second-order ten-sors defined in Rm. The displacement–strain rela-tion is given by
e½u ¼1 2 ru
n
þ ðruÞTo: ð2:1Þ
The conventional theory of ferromagnetics as-sumes that below the Curie temperature the mag-netization M varies with point x in X while maintaining a fixed magnitude:
jMðxÞj ¼ Ms 8x 2 X; ð2:2Þ
where Ms is called the saturation magnetization
which is a material constant depending on tem-perature only. Let
mðxÞ ¼ 1 Ms
MðxÞ : X ! S2; ð2:3Þ
where S2 denotes the three-dimensional sphere of
unit radius. The free or spontaneous magnetostrain tensor corresponding to the magnetization m is denoted as e0ðmÞ : X ! M33
s . In the case of cubic
materials the most general quadratic form of e0ðmÞ
is e0ðmÞ ¼3 2 k1 0 0 m ( m 1 3I þ ðk1 1 1 k1 0 0Þ X i6¼j mcimcjðec i e c jÞ ) ; ð2:4Þ where I2 M33
s is the identity tensor,fe c 1; e c 2; e c 3g is
an orthonormal set of crystal basis, mc
i are
com-ponents of m in the crystal basis, k1 0 0and k1 1 1are
independent, dimensionless material parameters indicating the strength of the magnetoelastic interaction. Above a b is the tensor product of two vectors a and b with components ða bÞij¼
aibj in some orthonormal basis. Note that
mag-netostrain e0ðmÞ is an even function in m; i.e.,
e0ðmÞ ¼ e0ðmÞ.
2.2. Magnetoelastic energy
Let oX denote the boundary of the domain X and vXðxÞ the characteristic function of X such that vX¼ 1 in X and vX¼ 0 outside X. Assume
oX¼ oX1[ oX2[ oX3 where oXi\ oXj¼ ; a.e.
for i6¼ j. Let u0:oX
1[ oX3! R3 and t0:oX2[
oX3! R3 be the applied displacement and
trac-tion on the boundary oX. The mechanical boun-dary conditions are for i; j¼ 1; 2; 3,
u¼ u0; x2 oX 1; rn¼ t0; x2 oX 2; u li¼ u0i; ðrnÞ lj¼ tj0; i6¼ j; x2 oX3; ð2:5Þ
where r : X! M33s is the stress tensor and
fl1; l2; l3g is an orthonormal set of vectors such
that l1and l2are tangent to the boundaryoX3and
l3¼ n is the unit outward normal to the boundary
oX3. Note that the boundary condition (2.5)3is the
mixed type and u0 t0¼ 0 there. Let H0:
R3! R3
be the applied magnetic field that would be present were the ferromagnetic crystal absent.
We postulate that below the Curie point we can obtain the magnetization as the one that minimizes m(x) H0 t0 x Ω
Fig. 1. A ferromagnetic crystal subjected to external magnetic field H0and mechanical loading t0.
EðmÞ ¼ Z X rm Arm þ WaðmÞ þ1 2 e e0ðmÞ C e e0ðmÞ dx Z X l0H 0 M dx Z oX2[oX3 t0 u dS þl0 2 Z R3 jr/j2dx ð2:6Þ subject to two constraints
r Hdþ MvX
¼ 0 Hd¼ r/ in R3; r r ¼ 0 r¼ C e e½ 0ðmÞ in X;
ð2:7Þ where / is the magnetic potential resulting from r Hd¼ 0, Hd :
R3! R3
the demagnetization or stray vector field, and the right-hand side of (2.7)2
is the linear constitutive assumption with C as the elastic modulus. We assume that C is the positive-definite symmetric fourth-order tensor such that for i; j; k; l¼ 1; 2; 3,
Cijkl¼ Cijlk¼ Cklij; ð2:8Þ
where Cijklare components of C in an orthonormal
basis. Above in (2.6) A2 M33
s is a
positive-defi-nite tensor, l0¼ 4p 107 N/A2 the coefficient of
permeability of free space, and Wa: S2! R the
anisotropy energy density whose expression in terms of lowest-order terms for cubic and uniaxial crystals is WaðmÞ ¼ Kc 0þ K c 1ðm c2 1m c2 2 þ m c2 2m c2 3 þ m c2 3m c2 1Þ cubic crystals; Ku 0þ K1u 1 ðm eÞ 2 n o uniaxial crystals; 8 > > > < > > > : ð2:9Þ where Kc 0, K u
0 are irrelevant constants such that
WaðmÞ P 0 for all m 2 S2, K1c, K u
1 are anisotropy
constants depending on materials and e is certain unit vectors.
Let us clarify the notation. Given M : X! MsS2, we solve (2.7) to obtain / and e
and plug them back to (2.6) to obtain EðmÞ. MaxwellÕs equation (2.7)1 is of course solved over
all space R3 with M¼ 0 outside X; we emphasize
this by writing MvXwhere vX is the characteristic
function of X explained above. Thus, by (2.7)1, we
mean:
r ðr/ þ MÞ ¼ 0 in X; r2/¼ 0 in R3n X;
s r/ þ Mt n ¼ 0 onoX; /! 0 as jxj ! 1;
where s t denotes the jumpacross the interface. Next, elastic equilibrium equation together with constitutive assumption in (2.7)2 can be solved
were the virtual distribution of body forces gen-erated by f¼ r Ce0ðmÞ applied to the body as
well as the constitutive relation changed as r¼ Ce. In other words, we mean
r rþ f ¼ 0; r¼ Ce; f¼ r Ce0ðmÞ:
ð2:10Þ Note that (2.10) is similar to the formulation of thermoelastostatics.
Each of the terms in the functional (2.6) has a physical interpretation. The first term, called the exchange energy, penalizes changes in the magne-tization, and thus is interpreted as the energy of forming a magnetic domain wall. The second and third terms, the (total) anisotropy energy, is the energetic cost that the crystal must pay if the magnetization and strain deviate from the pre-ferred states at that temperature; thus this builds in the information that the crystal prefers certain spontaneous magnetization and strain at a given temperature. Note that the zeros of Wa define the
easy axes, i.e., the directions along which the crystal is magnetized most easily. For example, for cubic crystals, the case Kc
1>0 in (2.9)1 indicates
the easy axes are along h1 0 0ic while K c 1<0 in
(2.9)1 indicates the easy axes are along the body
diagonals,h1 1 1ic. Similarly a positive Ku
1 in (2.9)2
describes an easy axis along the e direction while a negative Ku
1 shows an easy plane perpendicular to
the hard axis e for uniaxial crystals. The fourth, called the Zeeman energy, is the potential energy of the applied magnetic field, and this enforces the desire of the magnetization to align with the ap-plied magnetic field. The fifth is of course the po-tential energy of the mechanical loading device.
The final term, called the demagnetization energy or stray field energy, is the magnetostatic self-en-ergy associated with the magnetic field generated by the ferromagnetic crystal itself.
Finally, we consider a common case where the exchange tensor is isotropic; i.e., A¼ AI where A is called the exchange constant. In addition, we as-sume the constrained displacement boundary condition if it exists; i.e., u0¼ 0 8x 2 oX
1[ oX3. In
this case, let r be the auxiliary stress state
satis-fying r r ¼ 0; x2 X; r n¼ t0; x2 oX 2; ðr nÞ l j¼ t0j; x2 oX3; ð2:11Þ
where ljis the same as that defined in (2.5). Thus,
the total free energy of the ferromagnetic crystal (2.6) can be replaced by EðmÞ ¼ Z X Ajrmj2 þ WaðmÞ þ1 2 e e0ðmÞ C e e0ðmÞ dx Z X ðl0H 0 M þ r eÞ dx þl0 2 Z R3jr/j 2 dx ð2:12Þ
subject to the constraints (2.7). Suppose r ¼ r
0n n where r0>0 and n is the loading
direction. It follows that magnetization is pre-ferred to be aligned to maximize n en which is the projection of strain on the loading direction n n. Eq. (2.12) with the helpof the auxiliary stress r
is the micromagnetic formulation accounting for the effect of magnetoelastic interaction. It has been used by James and Wuttig (1998) and DeSimone and James (2002) to study various types of mag-netic domains and will be adopted in the rest of this paper. Note that we have used the total strain e in the fifth term of (2.12) instead of choosing magnetostrain e0ðmÞ proposed by Hubert and
Sch€afer (1998) (see Eq. (3.61) in p. 148 there). The idea of using the total strain comes from the ori-ginal formulation (2.6) which results from the principle of minimum potential energy in funda-mental mechanics (Dym and Shames, 1973). The
distinction between these two formulations is not significant for magnetic materials with extremely small magnetostriction (see Section 3.2 for expla-nation). However, it may lead to inconsistent results in certain materials with large magneto-striction (see Section 5.2).
2.3. Effective magnetic fields
The task in micromagnetics is to determine the magnetization mðxÞ which minimizes the free en-ergy (2.12) subject to the constraints (2.7). A number of techniques to solve this minimization problem for rigid ferromagnetics have been pro-posed in the past decades which all are based on the variational method proposed by Brown (1962). Due to the constraint equation given by (2.2), Brown has considered a small variation of the direction of the magnetization vector instead of a small variation of the magnetization by an arbi-trary function. The variational principle leads to modified BrownÕs equations accounting for the magnetoelastic effect mðxÞ HeffðxÞ ¼ 0 8x 2 X; om onðxÞ ¼ 0 8x 2 oX; Heff¼ 1 l0Ms dE dm ¼ Heþ Haþ H0þ Hsþ Hd; ð2:13Þ
where Heffis the effective magnetic field defined by the variational derivative of the free energy 1
l0Ms dE dm, H
e
the exchange magnetic field, Ha the anisotropic magnetic field, H0 the external mag-netic field, Hsthe stress-induced magnetic field and Hdthe demagnetization or stray field. The field H0 has been discussed in Section 2.2, Hd is obtained by solving MaxwellÕs equation (2.7)1, and the rest
of magnetic fields are defined by He¼ 2A l0Ms r2m; Ha¼ 1 l0Ms oWaðmÞ om ; Hs¼ 1 l0Ms C e e0ðmÞ oe 0ðmÞ om ; ð2:14Þ
where e is implicitly dependent on magnetization and is obtained by solving (2.7)2. Note that (2.13)
can be derived using the variational method simi-lar to that for rigid ferromagnetics (Aharoni, 1996) while the stress-induced magnetic field Hs is ob-tained by considering Z X d 1 2C½e½u e0ðmÞ ½e½u e0ðmÞr e½u dx ¼ Z X C½e½u e0ðmÞ ½de½u de0ðmÞ r de½udx ¼ Z X r ð r Þ de½u r de0ðmÞ dx ¼ Z X roe 0ðmÞ om dm dxþ Z oX ðr r Þn du ds ¼ Z X roe 0ðmÞ om dm dx
because ofr r ¼ 0, r r ¼ 0, boundary
condi-tions (2.5) and (2.11) and the fact that C is sym-metric.
The effective magnetic field Heff provides an
interpretation of the micromagnetic equations. From (2.13)1 either the magnetization vector m is
aligned along the effective magnetic field Heff at each point in the body X or the effective magnetic field Heff vanishes by itself since m6¼ 0 from (2.2).1 In both cases the torque exerted on any magnetization vector has to be zero in static equilibrium. Finally, we are aware that (2.13) is the necessary condition for the minimization of total free energy. Magnetic domains satisfying the mic-romagnetic equations given by (2.13) are the local equilibrium magnetization.
2.4. Magnetization dynamics
If BrownÕs equations (2.13) do not hold at some points, the nonzero torque m Heff
provides a precession of the magnetization around the effec-tive field. Indeed, the dynamic description of micromagnetic processes in a ferromagnetic material can be described by the famous LLG evolution equation (Landau and Lifshitz, 1935; Gilbert, 1955)
dm
dt ¼ cm H
eff acm ðm HeffÞ; ð2:15Þ
where t is the time. Above the first term on the right-hand side is the gyromagnetic term, with c 2:21 105m/A/s being the gyromagnetic ratio.
The second on the right-hand side is the damping term, with a being the dimensionless damping coefficient. The damping term allows the magneti-zation to turn towards the effective field until both vectors are parallel in the static solution. It is obvious that BrownÕs micromagnetic equations can be viewed as a particular case of (2.15), giving the static equilibrium when there is no change in time. The boundary conditions here are the same as that in the static equilibrium given by (2.13).
Clearly from (2.15) the torque m Heff
pro-vides the driving force for the evolution of mag-netization, and the small quantity of jdm
dtj implies
the (local) equilibrium magnetization configura-tion. In fact, the evolution of magnetization de-scribed by LLG equation (2.15) also implies the decrease of the total free energy (2.12), and this can be understood by considering
1 l0Ms dE dt ¼ Z X ðHeff _mÞ dx ¼ Z X Heff cm Heff acm ðm HeffÞ dx ¼ ac Z X
Heff ðm H effÞm HeffÞdx
¼ ac Z
X
f m H eff2 H eff2gdx 6 0
ð2:16Þ due to the famous Cauchy–Schwartz inequality andjmj ¼ 1. The more general form of dissipation
1For example, if the crystal is large enough to reach the large-body limit (DeSimone, 1993), the exchange energy can be neglected and in this case the closure four domains results in Heff¼ 0 in rigid magnetic cubic crystals (Hubert and Sch€afer,
mechanism other than LLG equation in rigid micromagnetics can be found in the recent work of Podio-Guidugli (2001).
As the result of (2.16), LLG equation cannot only provide the information of time evolution of magnetization, but also offers a convenient tool to study the static equilibrium magnetization. As a matter of fact, the integration algorithm of LLG and LLG-like equations of motions have been used widely in literature to study domain walls and equilibrium magnetization distributions by assuming dm
dt
is less than the prescribed small value (see Nakatani et al. (1989) and Berkov et al. (1993) for a detailed description). Thus from now on we will use (2.15) as our fundamental tool to investigate magnetic domain patterns and their evolution for various magnetostrictive crystals. Other methods to solve (2.13) directly can be found in Aharoni and Jakubovics (1986), LaBonte (1969) and Miltat et al. (1989).
3. Solutions of constrained equations
In magnetostatics the solution of the magnetic potential / in (2.7)1can be derived using the Green
function method introduced by the potential the-ory (Kellogg, 1969). There are no problems with boundary conditions or material properties as the Green function is independent of both of them. Unfortunately the analogous procedure meets certain difficulties for magnetoelastic interactions. Essentially we try to solve
r ðCe½uÞ ¼ r ðCe0ðmÞÞ 8x 2 X ð3:1Þ
which seems to be similar to r ðr/Þ ¼ r ðMvXÞ 8x 2 R
3
ð3:2Þ once the magnetization MðxÞ is given at each point in X. However, (3.1) is in general much more dif-ficult to be solved than (3.2) because of the com-plications associated with the elastic tensor C and boundary conditions. It is in principle not possible to extract the fourth-order tensor C out of the differential equation in (3.1). In addition, this tensor is material dependent, so is the solution of (3.1). Finally, the Green function for (3.1) is dependent on the domain X itself and the
associ-ated boundary conditions given by (2.5) while the Green function for the magnetic equation (3.2) is independent of domains and materials. As a result, there are no generally applicable procedures for the solution of the elastic case (3.1), at least, ex-plicit expressions similar to the magnetic solution. We are then forced to resolve such a difficulty ei-ther by a crude approximation (see Section 3.2) or by certain sophisticated numerical schemes (see Section 3.3).
3.1. Maxwell equation
The solution of the magnetic potential / in (2.7)1 can be easily obtained by Green function
method. We refer to Jackson (1962) and Kellogg (1969) for the detailed derivation and only list the formulation. Suppose the ferromagnetic crystal is a finite body occupying a domain X in the three-dimensional space, the magnetic potential is given by /ðxÞ ¼ 1 4p Z X rx0 Mðx0Þ jx x0j dx 0 þ Z oX Mðx0Þ n jx x0j dSx0 ; ð3:3Þ
while for an infinitely long ferromagnetic cylinder along the e3 direction, the magnetic potential is
given by /ðxpÞ ¼ 1 2p Z X rx0 p Mðx 0 pÞ lnjxp x0 pj dx 0 p Z oX Mðx0pÞ n lnjxp x0pj dSx0p ; ð3:4Þ where xp¼ ðx1; x2Þ and X here is the finite in-plane
cross section of the long cylinder. Note that the Green functions for (3.3) and (3.4) are the New-tonial potential 1
4p 1
jxj and logarithmic potential 1
2plnjxpj. Both functions are independent of solid
boundaries and material properties.
3.2. Equilibrium equation by approximation There are in general no explicit expressions similar to (3.3) or (3.4) for the solution of elastic
equilibrium equation (2.7)2. To overcome such a
difficulty and complexity, a very common ap-proach in literature is to replace the total strain e by magnetostrain tensor e0ðmÞ (Cullity, 1972; Cullen
et al., 1997; Callegaro and Puppin, 1997); i.e.,
e e0 ð3:5Þ in (2.12). As a result, (2.12) becomes EapproxðmÞ ¼ Z X Ajrmj2 n þ eWaðmÞ l0H 0 Mo dx þl0 2 Z R3jr/j 2 dx; ð3:6Þ where eWaðmÞ is the modified anisotropic energy
density defined by e
WaðmÞ ¼ WaðmÞ r e0ðmÞ: ð3:7Þ
Eq. (3.6) is the same as that for rigid ferromag-netics except that the anisotropy energy density WaðmÞ is replaced by eWaðmÞ. Note that (3.7) shows
that the energy of applied loads becomes another source of magnetic anisotropy energy while the magnetostrictive self-energy vanishes; i.e., the external stress causes the induction of a magne-tostrictive uniaxial anisotropy, usually called stress-induced anisotropy (Ross et al., 1996; Deng et al., 1997). The material is then treated as a rigid ferromagnetics with varying magnetic easy axes depending on external stress. This approximation, which we call relaxed approach, neglects the intrinsic stress due to the assumption of (3.5). Question arises for such an approximation because of loss of strain compatibility. Kinematic com-patibility requires (Gurtin, 1972)
r r e ¼ 0 ð3:8Þ while (3.8) may not be satisfied if eðxÞ is simply replaced by e0ðmðxÞÞ for arbitrary magnetization
at each point x.
The stress-induced magnetic field defined by (2.14) becomes Hrs¼ 1 l0Ms r oe 0ðmÞ om ; ð3:9Þ
where (3.9) is derived by the variational derivative of the relaxed energy given by (3.6) and the superscript r in Hrs emphasizes that this field is obtained by the relaxed approach. So clearly from
(3.9) Hrs is zero if the external stress r is absent;
and therefore, the ferromagnetic crystal is treated as a rigid material regardless of the possible exis-tence of intrinsic stress. Hence, this relaxed ap-proach is expected to be a crude approximation for magnetic materials with large magnetostriction at zero applied loads.
Finally, we remark on the current development of active materials such as shape-memory alloys, ferroelectrics and ferromagnetics. In the modern theory of these materials, the study of micro-structure requires the total strain eðxÞ to be re-placed by the stress-free strain e0ðxÞ at each point
if e0belongs to a class of strains S in which strains
are symmetry related and compatible with one another. The magnitudes of these strains do not need to be small, and in some materials such as shape-memory materials, they are very large. The intrinsic stress vanishes since strains in S are compatible. As a result, the effect of external stress on microstructure redistribution can indeed be approximated by considering the following po-tential energy:
r he0i;
where hi denotes the average over the body. The favorable microstructure is expected to the one which minimizes the above potential energy. While the results predicted by this theory have been found to be very useful in many applications, one drawback of this approach is that the energy barrier caused by intrinsic stress as microstructure evolves may prevent favorable states from being achieved, and the estimate of this barrier is in general not revealed by this approach. The sys-tematic description of this theory for various ac-tive materials can be found in Bhattacharya (1993), Shu and Bhattacharya (2001), Li and Bhattacharya (in preparation) and DeSimone and James (2002).
3.3. Equilibrium equation by modified boundary element method
As described in the introduction of Section 3, it is not an easy task to solve the elastic equilibrium equation given by (2.7)2 for arbitrary
two common cases: one is for a thin film and the other for a long cylinder with arbitrary finite cross section. For both cases, we may assume
m¼ mðx1; x2Þ 2 S2; u¼ uðx1; x2Þ 2 R3: ð3:10Þ
The problem is now simplified to be two-dimen-sional in space, but the magnetization m is still free to rotate in three dimensions. The simplification from (3.10) does not lead to the direct solution of (2.7)2 in closed forms. Instead, it provides an
opportunity to solve the elastic equilibrium equa-tion by certain sophisticated numerical schemes. There are two major numerical methods: finite element method (FEM) and boundary element method (BEM). Both methods have their own merits and are applied to a variety of scientific problems according to their specific needs. Tradi-tionally both methods have better numerical pre-cision in displacement u than strain e (Oden and Reddy, 1976). However, in our present case, the main variable we need is strain e rather than dis-placement u (see (2.14)3). We thus need to look for
a new strain-based formalism. Indeed, Wu et al. (1992) and Wu (2000) have proposed a new boundary integral equation formalism for two-dimensional anisotropic elasticity without consid-ering the effect of body force (i.e., f ¼ 0 in (2.10)). This formulation is presented in complex-variable using Stroh formalism for anisotropic elasticity (Stroh, 1958). The advantage of it is that in numerical implementation the boundary integrals can be integrated analytically along each element. In addition, it can be shown that this formalism has better numerical precision in variable strain as well as faster rate of convergence than conven-tional displacement-based schemes. We thus pro-pose to solve (2.7)2by extending their approach to
the case accounting for the effect of virtual body force resulting from the divergence of magneto-stress (i.e., f 6¼ 0 in (2.10)).
Let an elastic body occupy a long cylindrical domain with finite in-plane cross section denoted by X as shown in Fig. 2. Note that we have slightly abused the notation X which was used to denote the domain occupied by the whole ferromagnetic body. Suppose there are no body force densities and stress is related to strain by r¼ Ce where C is the elastic modulus satisfying (2.8). Let ak and pk
be the eigenvectors and eigenvalues of the follow-ing equation: Q þ ðR þ RTÞp kþ Tp2k ak ¼ 0
ðno sum over kÞ; ð3:11Þ where
Qik¼ Ci1k1; Rik¼ Ci1k2; Tik¼ Ci2k2;
i; k ¼ 1; . . . ; 3;
and Cijkl are the components of the elasticity
ten-sor C in some reference basis. As pk can be shown
to be nonreal (Stroh, 1958; Ting, 1996), we thus take pk, k¼ 1; 2; 3, with positive imaginary part for
definiteness. Let A¼ ða1ja2ja3Þ be the 3 · 3 matrix
whose ith column is replaced by the column vector ai defined by (3.11). Further, we define the matrix
B related with A by
B¼ RTAþ TAP; ð3:12Þ
where P¼ diag½p1; p2; p3; i.e., a matrix with
diag-onal entries replace by p1; p2; p3and zero in the rest
of entries. Set
z¼ x1þ ix2; f¼ n1þ in2; ð3:13Þ
where z and f are field and source points in the complex plane and i¼pffiffiffiffiffiffiffi1 is the imaginary number. Now letðx1; x2Þ 2 X, ðn1;n2Þ 2 oX and l; n
are unit tangent and outward normal vectors to the boundary oX as shown in Fig. 2(a). Letou
osðzÞ
and tðzÞ be the displacement gradient and traction at point z along an arbitrary contour s. Both can be related by the tangential displacement gradient
Fig. 2. (a) An elastic body occupies a long cylindrical domain with in-plane cross section denoted by X. (b) The strain e0 is assumed to be piecewise uniform in subdomains Xi X, i¼ 1; . . . ; B. Only three subdomains are shown here.
dl¼ouol and traction tnon the boundary oX by the
following dual set of boundary integral equations proposed by Wu et al. (1992): bou osðzÞ ¼ Z oX oU os ðz; fÞtnðfÞ oW os ðz; fÞdlðfÞ dl; ð3:14Þ btðzÞ ¼ Z oX oW T os ðz; fÞtnðfÞ þoV osðz; fÞdlðfÞ dl; ð3:15Þ where dl¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dn21þ dn22 q , b¼ 1 in X and b ¼1 2 at
the smooth boundaryoX, and Uðz; fÞ ¼ R½AGðz; fÞAT; Wðz; fÞ ¼ R½AGðz; fÞBT;
Vðz; fÞ ¼ R½BGðz; fÞBT:
Above R stands for the real part, A and B are defined by (3.11) and (3.12), and the matrix func-tion Gðz; fÞ is defined by Gðz; fÞ ¼ 1 ip lnðz1 f1Þ 0 0 0 lnðz2 f2Þ 0 0 0 lnðz3 f3Þ 0 @ 1 A; ð3:16Þ where zk¼ x1þ pkx2; fk ¼ n1þ pkn2: ð3:17Þ
Eqs. (3.14) and (3.15) with b¼1
2provide a pair
of boundary integral equations for the tangential displacement gradient and traction if the contour s is chosen to coincide with the boundary. Either (3.14) or (3.15) can be solved to obtain the un-known tangential displacement gradient or trac-tion. Once the unknown boundary data is determined, the displacement gradient and traction in any directions inside the domain X can therefore be computed from (3.14) and (3.15) with b¼ 1.
Now suppose that the constitutive relation be-comes rðxÞ ¼ C½eðxÞ e0ðxÞ where e0 is the
ei-genstrain. Assume there is a set of disjoint domains Xi X, i ¼ 1; . . . ; B, such that e0ðxÞ ¼ 0
if x2 X n ðX1[ [ XBÞ and is uniform in Xi as
shown in Fig. 2(b). The eigenstrain e0 may take
different values for different subdomains Xi. As the
traction is balanced (continuous) across the boundaryoXi, we have lim x2Xi x!oXi rni lim x2XnXi x!oXi rni¼ r0ni; ð3:18Þ
where ni is the unit outward normal to the
boundaryoXi, and
r¼ Ce; r0¼ Ce0: ð3:19Þ
Physically, it means there is a prescribed body force density r0ðxÞn
iðxÞvoXi supported on the
boundary of each Xi if the constitutive law were
replaced by r¼ Ce throughout the whole domain X. To account for such an effect, note that Wu et al. (1992) have pointed out that the displacement gradient (3.14) and traction (3.15) along an arbi-trary contour s in a finite body X can be viewed as those due to the dislocation dipoleu dl and body force density tn supported on the boundaryoX in
an infinite space. As a consequence, the displace-ment gradient and traction along any arbitrary contour s due to piecewise uniform eigenstrain e0ðxÞ can be obtained from (3.14), (3.15), (3.18),
and (3.19) and are given by bou osðzÞ ¼ Z oX oU osðz; fÞtnðfÞ oW os ðz; fÞdlðfÞ dl þX i Z oXi oU os ðz; fÞt 0 niðfÞ dl; ð3:20Þ b tðzÞþX i t0ðzÞvXi ! ¼ Z oX oW T os ðz; fÞtnðfÞ þ oV osðz; fÞdlðfÞ dl X i Z oXi oWT os ðz; fÞt 0 niðfÞ dl; ð3:21Þ where t0ðzÞ ¼ ðCe0ðzÞÞn0; t0 niðfÞ ¼ ðCe 0ðfÞÞn i; where n0and n
iare the unit normals to the contour
s and the boundary oXi, respectively. The
un-known tangential displacement gradient or trac-tion can be obtained by solving (3.20) and (3.21) if the contour s is chosen to coincide with the boundary. Once the unknown boundary data is
determined, the displacement gradient and trac-tion in any directrac-tion inside X can therefore be computed from (3.20) and (3.21) with b¼ 1.
In magnetostrictive material, the magnetostrain e0ðmÞðxÞ is in general not piecewise constant inside
X. However, we may assume that it is uniform in each small cell in which the magnetization m is approximated to be constant. In this situation, we can use (3.20) and (3.21) to solve (2.7)2numerically
and the details are provided in the next section. Finally, the previous analysis is applied to a long cylinder with finite cross section. In the case of thin films, (3.20) and (3.21) are still valid with slight modification. Under the plane stress assumption, the elastic modulus C is replaced by the generalized plane stress modulus Cp, and only the in-plane components of magnetostrain e0 are
taken into account in the formulation of both (3.20) and (3.21). To make it clear, let the film occupy the domain X¼ S ð0; hÞ where h is the film thickness and S the in-plane cross section of the film. For simplicity, assume that S is a square domain; i.e., S¼ ð0; LÞ2. Further, decompose u¼ upþ ~bx 3 where up¼ ðu p 1ðx1; x2Þ; u p 2ðx1; x2Þ; 0Þ and ~b¼ ð2~b1ðx1; x2Þ; 2~b2ðx1; x2Þ; ~b3ðx1; x2ÞÞ. Set b 1¼ ~ b1 e031, b 2¼ ~b2 e032 and b 3¼ ~b3 e033 and ðb
1; b 2; b 3Þ is the minimizer of infb2R3 12E CE
where b¼ ðb1; b2; b3Þ and E¼ ep11 e0 11 e p 12 e012 b1 ep21 e0 21 e p 22 e022 b2 b1 b2 b3 0 @ 1 A: ð3:22Þ Let u¼u L and x¼ x
L. It follows from Shu (2002)
that the total magnetostrictive energy per unit volume can be shown to be
1 hL2 Z X 1 2ðe e 0Þ Cðe e0Þ dx ¼ U0þ h L U1þ o h L ; ð3:23Þ whereoðaÞa ! 0 as a ! 0 and
U0¼ Z 1 0 Z 1 0 1 2 e p e0p Cp ep e0pdx 1dx2; U1¼ Z 1 0 Z 1 0 2 eb Cp ep e0p dx1dx2: ð3:24Þ
Above, Cpis the generalized plane stress modulus, e0p2 M22
s the in-plane magnetostrain tensor with
components e0pab ¼ e0
ab, ep2 M 22
s the in-plane
strain tensor with components defined by
epab¼1 2 oup a oxb þou p b oxa ; and eb2 M22
s with components defined by
eb ab¼ 1 2 ob a oxb þob b oxa
for a; b¼ 1; 2. As a consequence from (3.23), we only need to consider the energy U0if the ratioLhis
small, and minimizing it yields the reduced elastic equilibrium equation rp rp¼ 0 where rp¼
Cp½ep e0p and r
pis the in-plane gradient. Notice
that U1 can be completely determined once U0 is
known.
4. Numerical algorithm
4.1. Integration of LLG equation
The integration of LLG equation (2.15) is based on the algorithm provided by the National Insti-tute of Standards and Technology (NIST) of the US (Donahue, 2000) which has been used exten-sively in a variety of applications and fundamental research (Donahue, 1998; Crew and Lewis, 2001; Thiaville et al., 2002). It is based on the finite difference method by replacing the continuous solution of magnetic domains by a discrete set of lattice. At each lattice point the differential oper-ators are replaced by finite difference operoper-ators using a second-order predictor–corrector tech-nique of the Adams type, and the conditions on the boundary of magnetic domains are replaced with their discrete counterparts.
The exchange energy is evaluated using the eight-neighbor bilinear interpolation (Donahue and McMichael, 1997), the anisotropy and Zee-man energy terms are computed assuming con-stant magnetization in each mesh. The stray field is calculated as the convolution of the magnetization against a kernel describing the mesh to mesh magnetostatic interaction. The convolution is computed using Fast Fourier Transform. The
stress-induced magnetic field is based on the modified BEM described in Section 3.3 and the algorithm is provided in Section 4.2.
The modeling system is a two-dimensional thin film with film thickness much smaller than its lat-eral dimensions. Thus the plane stress condition is used throughout the calculation. We also assume that the film thickness is on or smaller than the order of magnetic exchange length defined by p ffiffiffiffiA
K1
q
where A; K1 are exchange and anisotropy
parameters (see (2.12) and (2.9). So we may pos-tulate that the magnetization is the function of in-plane variables x1, x2only. The film is divided into
a large number of small square cells. Inside each cell the magnetic moment m is free to rotate in three dimensions but its magnitude is kept to be unity. The Neumann type of magnetic boundary condition given by (2.13) is assumed throughout the calculation.
4.2. Stress-induced magnetic field
We recall that the stress-induced magnetic field defined by (2.14) is Hcs¼ 1 l0Ms roe 0ðmÞ om ; ð4:1Þ
where the superscript ÔcÕ denotes that this field is obtained by the constrained approach in which the constrained (equilibrium) equation (2.7)2has to be
solved for each magnetization configuration. Note that Hcs is different from H
r
sdefined by the relaxed
approach (3.9). We have implemented two algo-rithms to compute Hc
sand H r
sin order to compare
the differences between these two methods. As Hr s
is quite easy to be implemented, we now focus on the numerical implementation of Hcs next.
We first concentrate on the term oeom0ðmÞ. Recall thatfec
1; e c 2; e
c
3g is an orthonormal crystal basis and
fe1; e2; e3g the orthonormal basis in the reference
configuration. They are related by eci ¼ Rei; Rij¼ ei Rej¼ ei ecj;
i; j¼ 1; 2; 3; ð4:2Þ where R is the proper orthogonal tensor and Rij
are the components of the tensor R in the basis ei.
Let mc
i and mi be the components of the
magneti-zation vector in the bases ec
i and ei, respectively.
They are related by mc
j¼ Rijmi. Further, let e0 c ij and
e0
ij be the components of strain tensor e
0 in the
bases ec
i and ei, respectively. The transformation
rule for tensors in different bases yields e0 ij¼ Ripe0 c pqRjq. It follows that oe0 ij omk ¼ RipRjq oe0c pq omk ¼ RipRjq oe0c pq omc s omc s omk ¼ RipRjqRks oe0c pq omc s ; ð4:3Þ
where the general expression of e0 in terms of
components of m in the crystal basis is given by (2.4) for cubic materials. Notice that summation over repeated indices is implied in (4.3) unless noted otherwise.
We next focus on calculating the stress r¼ C½e e0ðmÞ obtained by solving the equilibrium
equation (2.7)2. In Section 3.3 we propose a dual
set of equations (3.20) and (3.21) to calculate the displacement gradient and traction along an arbi-trary contour s inside a body X in which the magnetostrain is assumed to be piecewise uniform. The first stepis to calculate the boundary dis-placement gradient and traction by letting the contour s coincide with the boundaryoX. Let the boundaryoX and subboundaries oXi, i¼ 1; . . . ; B,
be approximated by N1and N2ðiÞline segments over
which the normal tractions and tangential dis-placement gradients are assumed to be constant. The resulting discretized equations of (3.20) and (3.21) with the midpoints of the line segments chosen as the collocation points are
XN1 k¼1 ðW ÞjkðdlÞk ¼X N1 k¼1 ðU ÞjkðtnÞkþ XB p¼1 XN2ðpÞ q¼1 ðU Þjqðt 0 nÞ p q; j¼ 1; 2; . . . ; N1; ð4:4Þ XN1 k¼1 ðW TÞjkðtnÞk ¼X N1 k¼1 ðV Þ jkðdlÞk XB p¼1 XN2ðpÞ q¼1 ðW T Þjqðt0 nÞ p q; j¼ 1; 2; . . . ; N1; ð4:5Þ
where ðdlÞk and ðtnÞk are the tangential
displace-ment gradients and normal tractions at the kth element and the elements are numbered consecu-tively along the tangential direction l of the boundaryoX, andðt0
nÞ p q¼ ½ðCe
0Þnp
qare the normal
tractions at the qth element of the subboundary oXp. If j6¼ k, the matrices U , W and V for
anisotropic materials are given by the following analytic expressions: ðU Þjk ¼ R½AðG ÞjkAT; ðW Þ jk¼ R½AðG Þ jkB T; ðV Þjk¼ R½BðG Þ jkB T ; ð4:6Þ where ðG Þ jk¼ Diag ½d 1 jk;d 2 jk;d 3 jk are diagonal
matrices with d1jk, d2jk and d3jk given by
d1jk¼1 ip ^fðjÞ 1 ^ fðkÞ1 ln fðjþ12Þ 1 f ðkþ1Þ 1 fðjþ12Þ 1 f ðkÞ 1 ! ; d2jk¼1 ip ^fðjÞ 2 ^ fðkÞ2 ln fðjþ12Þ 2 f ðkþ1Þ 2 fðjþ12Þ 2 f ðkÞ 2 ! ; d3jk¼1 ip ^fðjÞ 3 ^ fðkÞ3 ln f ðjþ1 2Þ 3 f ðkþ1Þ 3 fðjþ12Þ 3 f ðkÞ 3 ! : ð4:7Þ
Above in (4.7) ^fðjÞa ¼ cos hjþ pasin hjwith hjbeing
the angle between the jth element and the x1-axis,
fðkÞa ¼ nðkÞ1 þ panðkÞ2 withðn ðkÞ 1 ;n ðkÞ 2 Þ and ðn ðkþ1Þ 1 ;n ðkþ1Þ 2 Þ
as the end points of the kth element, fðjþ12Þ a ¼ n ðjþ1 2Þ 1 þ pan ðjþ1 2Þ 2 and ðn ðjþ1 2Þ 1 ;n ðjþ1 2Þ 2 Þ is the
midpoint of the jth element (see Fig. 3 for the detailed notation). If the closure of subdomain Xp
is completely contained in X, the matrices ðU Þjq
andðW Þjq for anisotropic materials are the same
as those in (4.6) except ^fðqÞa and fðqÞa denote points in
the qth element of subboundaryoXp. On the other
hand, for j¼ k or if some parts of the subboun-daries oXp coincides with the boundary oX (and
j¼ q in this case), ðU Þjj¼ ðV Þ jj¼ 0; ðW Þ jj¼ 1
2I ðno sum over jÞ: ð4:8Þ Finally, if the elastic material is isotropic, the ex-plicit expression of these three matrices U , W and V can be found in Wu et al. (1992).
Either (4.4) and (4.5), which are referred to as type 1 and type 2 equations, can be used to solve unknown boundary variables numerically. Once all the boundary data including dl and tn are
ob-tained, we can use (3.20) and (3.21) together with (4.6) and (4.7) to calculate the stress at any points inside X. Examples to show the applicability of the formalism given by (4.4) and (4.5) are given in Appendix A.
5. Results
In this section we explore the effect of intrinsic stress on the equilibrium magnetic domains and hysteresis calculated using the constrained ap-proach. We also compare our results with those obtained based on the relaxed approach (i.e., the negligence of the influence of internal stress). We use Ni and Terfenol-D as the representative materials.
5.1. Extremely small magnetostriction
Ni is one of the common ferromagnetic mate-rials with 105order of magnitude of spontaneous
magnetostriction. The associated magnetic mate-rial parameters are Ms¼ 4:8 105 A/m, A¼
9 1012 J/m, Kc
1¼ 5:7 103 J/m3, and k1 0 0¼
4:6 105, k
1 1 1¼ 2:4 105. The cubic
aniso-tropic elastic constants in the crystal basis are C11¼ 2:5 1011 N/m2, C12¼ 1:6 1011 N/m2,
C66¼ 1:18 1011N/m2in terms of Voigt notation.
We consider a single crystal Ni film with 100 nm· 100 nm · 10 nm dimensions. Assume the crystal basis coincides with the reference basis. There are no external stress and applied field, and
(
ξ
1(k+1/2),
ξ
2(k+1/2))
(
ξ
1(k+1),
ξ
2(k+1))
x
1θ
ktherefore, r ¼ 0 and H0¼ 0 in this case. The
initial magnetization is mðxÞ ¼ e2 and the final
results of magnetization evolution described by LLG equation are shown in Fig. 4(a) and (b) un-der the conditionjm heffj < where ¼ 104and
heff ¼ 1 MsH
eff. We have also used the smaller
parameter ¼ 105as the criterion and found that
the results are almost the same as those obtained using the original criterion. Note that Fig. 4(a) and (b) are magnetization patterns obtained using the relaxed and constrained approach. Obviously the results show that there are no significant differ-ences between these two methods; in other words, the effect of intrinsic stress on magnetic domain patterns is not significant for magnetostrictive materials with extremely small magnetostriction.
We next investigate the hysteresis property such as coercivity of Ni thin films under intrinsic stress. The external field H0 is applied along the e2
direction with cycling magnitudes. Fig. 5(a) shows the simulated hysteresis loops calculated based on the relaxed (dashed line) and constrained (contin-uous line) approaches. The overall loops appear to be similar while there is a difference in coercive fields. Let Hr
c and H c
c be the coercive magnetic
fields obtained by relaxed and constrained meth-ods. According to Fig. 5(a) Hr
c ¼ 7160 A/m and
Hc
c ¼ 8750 A/m. The error estimate is
jHr c Hccj Hc c 100 18%: It is expected that Hc c is larger than H r c since more
energy input is needed to overcome the increase of elastic energy due to incompatible magnetostrain during magnetization rotation against reversal magnetic field.
(a) (b)
(c) (d)
Fig. 4. Magnetic domain patterns for Ni films simulated based on the relaxed (a) and constrained (b) approaches. Magnetic domain patterns for Terfenol-D films simulated based on the relaxed (c) and constrained (d) approaches.
-80000 -40000 0 40000 80000 -1 -0.5 0 0.5 1 R C -80000 -40000 0 40000 80000 -1 -0.5 0 0.5 1 R C H0 (A/m) H0 (A/m) m m (a) (b)
Fig. 5. Hysteresis loops calculated based on the constrained (C) and relaxed (R) approaches: (a) Ni thin films; (b) Terfenol-D thin films.
5.2. Large magnetostriction
Terfenol-D (TbxDy1xFe2, x 0:3) is one of
giant magnetostrictive materials with 103order of
magnitude of magnetostriction. The relevant magnetic properties of Terfenol-D are Ms¼ 8105
A/m, A¼ 9 1012 J/m, Kc
1 ¼ 6 10
4 J/m3, and
k1 0 0 0, k1 1 1¼ 1:64 103 (Clark, 1986; Lord
and Harvey, 1994). The bulk cubic elastic con-stants in the crystal basis are C11¼ 1:41 1011N/
m2, C
12¼ 6:48 1010 N/m2 and C66¼ 4:87 1010
N/m2 in terms of Voigt notation (Dewar, 1997).
5.2.1. Uniform magnetization
Consider a single crystal Terfenol-D thin film with 100 nm· 100 nm · 10 nm dimensions. The crystal orientation is assumed to be½1 1 1cjje1 and
½1 1 2cjje2 as shown in Fig. 6(a). There are no
ap-plied traction and field in this case; therefore, we take r ¼ 0 and H0¼ 0 in (2.13). The initial
magnetization pattern is m¼ e1 and the final
re-sults of magnetization evolution described by LLG equation are shown in Fig. 4(c) and (d) under the same criterion described in the case of Ni. Note that Fig. 4(c) and (d) are magnetization patterns simulated using the relaxed and constrained ap-proaches.
In contrast to the case of Ni films, results sim-ulated based on these two methods are completely different as can be seen in Fig. 4(c) and (d). In the former case, the elastic energy is not taken into
account due to the negligence of intrinsic stress. So magnetization vectors are free to rotate and prefer to be oriented in one of the easy axes. As Kc
1<0,
h1 1 1ic is the easy direction for Terfenol-D.
How-ever, only in-plane easy axes including½1 1 1c and
½1 1 1c are taken into account since Gioia and
James (1997) have shown that the in-plane mag-netization is the preferred state for very thin fer-romagnetic films to reduce demagnetization energy. Combining all these facts leads to domain patterns shown in Fig. 4(c): magnetization vectors are aligned to ½1 1 1c direction on both right and left sides of the film because of reduction of the stay field energy, and to ½1 1 1c in the middle to reduce the anisotropy energy.
On the other hand, the elastic energy plays an important role in equilibrium domain patterns for magnetic materials with large magnetostriction. The consideration of intrinsic stress due to mag-netization rotation and boundary conditions pre-vents magnetization vectors from deviating too much from their initial easy direction ½1 1 1c
dur-ing evolution. As a result, most of magnetization vectors in Fig. 4(d) are oriented to½1 1 1cwhile the
rest of them on both right and left sides of the film deviate a little from the easy direction to reduce the stay field energy.
Finally we study hysteresis of Terfenol-D thin films. The external field H0is applied along the e1
direction with cycling magnitudes. Fig. 5(b) shows the simulated hysteresis loops calculated based on
35.3
˚
35.3˚
19.5˚
c ce
2e
3 [111] [112] c ce
1e
2e
3 [001] [110] [111] c [111] c _e
1 __ [111]c [110]c _ [110]c _(a)
(b)
the relaxed (dashed line) and constrained (contin-uous line) approaches. In contrast to the results shown in Fig. 5(a) and (b) shows a wide difference in hysteresis loops simulated using these two dis-tinct methods. Let Hr
c and Hcc be the coercive
magnetic fields obtained by relaxed and con-strained approaches. According to Fig. 5(b) Hr
c¼ 23,900 A/m and Hcc¼ 49,300 A/m, and the
error estimate is jHr c H c cj Hc c 100 52%:
In this case we have found that the coercive magnetic field Hc
c is twice larger than H r
c obtained
without consideration of intrinsic stress.
5.2.2. Twin
The magnetic domain patterns simulated in Section 5.2.1 are not commonly observed in Ter-fenol-D crystals at zero applied field. Instead, the observed microstructure is very complicated in the MFM image (Lord et al., 1997; Schmidt et al., 1998). Most of domain patterns appear in the way which we call ‘‘twins’’ (James and Kinderlehrer, 1993). We thus investigate the effect of intrinsic stress on this special microstructure.
Consider a single crystal Terfenol-D thin film with 200 nm· 200 nm · 10 nm dimensions. The crystal orientation is assumed to be½1 1 0cjje1 and
½0 0 1cjje2 as shown in Fig. 6(b). Assume r ¼ 0
and H0¼ 0 in (2.13). There are two easy axes in the plane of the film: one is along ½1 1 1c and the other is along½1 1 1c. The angle between the easy axis and e1 direction is about 35.3 as also shown
in Fig. 6(b). The twinning plane is ð0 0 1Þc plane associated with the½1 1 1cand½1 1 1 magnetization pair. Let m1¼ 2 ffiffiffi 6 p e1þ 1 ffiffiffi 3 p e2; m2¼ 2 ffiffiffi 6 p e1þ 1 ffiffiffi 3 p e2 ð5:1Þ
be aligned to the easy directions½1 1 1cand½1 1 1c,
respectively. So the anisotropy energy is zero. From (2.4) and the crystal orientation shown in Fig. 6(b), the associated components of magneto-strain in the reference basis are
e0ðm 1Þ ¼ k1 1 1 2 1 pffiffiffi2 0 ffiffiffi 2 p 0 0 0 0 1 0 B @ 1 C A; e0ðm 2Þ ¼ k1 1 1 2 1 pffiffiffi2 0 pffiffiffi2 0 0 0 0 1 0 B @ 1 C A: ð5:2Þ
Set ^n¼ e2. The conditionðm1 m2Þ ^n¼ 0 implies
that there are no internal ‘‘magnetic charges’’ on the twin plane and the demagnetization energy is re-duced to zero if the boundary effect is neglected. In addition, choosing a¼pffiffiffi2k1 1 1e1, we find e0ðm1Þ
e0ðm
2Þ ¼ a n þ n a satisfying the compatibility
condition (Silhav y, 1997). So e0ðm
1Þ and e0ðm2Þ are
compatible strain pair and, therefore, the intrinsic stress vanishes, so is elastic energy.
The initial magnetization pattern is m¼ m1
for 0 < x2<100 nm and m¼ m2 for 100 nm <
x2<200 nm. The inclusion of exchange energy
and boundary effect gives the final magnetic do-main patterns shown in Fig. 7(a) and (b) based on the relaxed and constrained methods. Note that we have found the identical simulation results for different initial configurations as long as the devi-ations of initial magnetization vectors from the easy directions are kept within 15. In order to see the differences between these two methods, let hr and hc be the angles of the magnetization vectors from the x1-axis in Fig. 7(a) and (b). Fig. 8(a) is the
contour plot showing the error estimate based on the formula jhrhhccj 100. The maximum error
estimate is upto 30% concentrated in the middle lower and upper parts of the film. However, this error estimate would be meaningless if hc is close to 90 since we can take 90 hc as another ref-erence angle. To see this, Fig. 7(a) and (b) shows the significant differences in angles in the domain wall region on the left-hand side of the film. But such a change in magnetization angles cannot be reflected by the original error estimate shown in Fig. 8(a). Therefore, we define another error esti-mate based on the formulamaxjhjhrhrhcjcj 100, and the
result is shown in Fig. 8(b). The new comparison diagram not only shows the trend of the error estimate but also reflects the change inside the wall region simulated by different approaches.
The average angle from the x1-axis away from
the domain wall region is around 30.9 in Fig. 7(a) while it is around 34.8 in Fig. 7(b). Clearly the magnetization vectors in Fig. 7(b) stay at much closer to the easy axes than those in Fig. 7(a) and it is also true inside the wall region. The negligence of elastic energy in Fig. 7(a) results in significant deviation from the easy directions for magnetiza-tion vectors near the boundary to reduce
demag-netization energy. On the other hand, elastic energy becomes dominant in the case of large magnetostriction. The inclusion of elastic energy in Fig. 7(b) keeps magnetization vectors be aligned to the easy directions as close as possible at the cost of increase of stay field energy. This agrees with the ‘‘constrained theory’’ proposed by DeSimone and James (2002). In that theory the total energy is minimized over the special class of functions in
(a) (b)
Fig. 7. Equilibrium twins for Terfenol-D films: (a) simulation based on the relaxed approach; (b) simulation based on the constrained approach.
(a) (b)
Fig. 8. (a) Error estimate for the twin structure based on the formula defined byjhrhhccj 100. (b) Error estimate for the twin structure based on the formula defined bymaxjhjhrhrhcjcj 100.
which magnetization vectors are assumed to be oriented only along the easy directions. Our sim-ulation confirms that this assumption turns out to be a reasonable approximation for materials with high anisotropy and large magnetostriction.
Finally, we investigate the hysteresis property. Let an external magnetic field be applied cyclically along the diagonal (45) direction of this square film. The average magnetization projected on the diagonal direction versus applied field is shown in Fig. 9(a) calculated based on the relaxed (dashed line) and constrained (continuous line) approaches. Again we see the wide difference in the overall loops calculated based on these two methods. Fig. 9(b) shows the hysteresis of the average shear strain versus the applied magnetic field. The simulation based on the relaxed approach shows an interme-diate stepduring the exchange of stability while this intermediate stepdisappears in the constrained approach. When the applied field decreases to zero from positive value, most of magnetization vectors are aligned to ½1 1 1c easy direction, resulting in
positive shear strain as shown in Fig. 10(a). From
(5.2) the angle due to shearing is about c 2¼
e2 e0ðm1Þe1¼p1ffiffi2k1 1 1¼ 0:11% for this Terfenol-D
film. During the exchange of stability, the relaxed approach which neglects the intrinsic stress shows a stable intermediate step: most of magnetization vectors stay close to½1 1 1ceasy direction, causing
the large negative shear strain as shown in Fig. 10(b). However, there is no such an intermediate step in the constrained approach. Instead, the magnetization is immediately rotated to ½1 1 1c direction when the reversal applied field exceeds to coercive magnetic field as shown in Fig. 10(c). This gives rise to positive magnitude of shear again as magnetostrain e0ðmÞ is an even function of
mag-netization (see (2.2)).
6. Conclusion
We have explored the effect of intrinsic stress on the magnetostrictive behavior of ferromagnetic thin films. We start with a theoretical framework based on micromagnetics accounting in detail for the effect of stress. The micromagnetic formulation is a nonlocal, nonconvex variational problem subject to two constrained equations: MaxwellÕs and elastic equilibrium equations. The conven-tional approach which relaxes the elastic equilib-rium equation takes into account only the external stress. The intrinsic stress caused by incompati-ble magnetostrain as magnetization rotates, how-ever, is neglected by this method. In contrast, we
-100000 -50000 0 50000 100000 -1 -0.5 0 0.5 1 R C -100000 -50000 0 50000 100000 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 R C m Shear Strain (%) H0 (A/m) H0 (A/m) (a) (b)
Fig. 9. (a) Hysteresis loops of average magnetization versus applied magnetic field calculated using the constrained (C) and relaxed (R) approaches. (b) Hysteresis loops of the average shear strain versus applied field simulated using the constrained (C) and relaxed (R) approaches.
−γ
γ γ
(a) (b) (c)
Fig. 10. The average shear deformation caused by magnetiza-tion rotamagnetiza-tion from (a) to (b) to (c). Note that magnetostrain e0ðmÞ is an even function of magnetization m.
propose a modified boundary integral formalism to solve the elastic constrained equation. This formulation has the advantage of higher accuracy in the calculation of stress as well as faster rate of convergence than conventional displacement-based numerical schemes. We further implement it into the LLG solver to study both equilibrium and transient magnetization configurations.
We use Ni and Terfenol-D as the model mate-rials. The former has small magnetostriction of the order of 105 while the latter has large
mag-netostriction of the order of 103. Figs. 4 and 5
highlight the striking contrast in equilibrium mag-netization and hysteresis calculated based on the constrained and relaxed methods. Our simulation shows that the relaxed approach turns out to be a reasonable approximation in predicting domain patterns and their evolution in the case of Ni films while it becomes a crude estimation for Terfenol-D films. In addition, we consider an example of twin which is the most observable microstructure in magnetostrictive materials. Fig. 8 shows the sig-nificant error estimate for the simulation of twins using the relaxed method. The result shows that away from the domain wall region the inclusion of magnetostrictive energy keeps magnetization vec-tors to be aligned to the easy axes as close as possible. This agrees with the ‘‘constrained the-ory’’ of DeSimone and James (2002) applied for materials with high anisotropy and large magne-tostriction.
Acknowledgements
The authors are glad to acknowledge the finan-cial support of the National Science Council of Taiwan under the grant NSC 90-2622-E-002-006.
Appendix A
We consider two examples to demonstrate the accuracy of the formalism given by (3.20) and (3.21). First we recall that the constitutive relation is r¼ Cðe e0Þ where e0 is the piecewise uniform
eigenstain. A new kind of patch test of compati-bility with geometry and boundary conditions is introduced in Fig. 11(a). The plane stress elastic constants used here are C11¼ 1:4767 1011N/m2,
C12¼ 0:5766 1011 N/m2, C66¼ 1:1848 1011 N/
m2 in terms of Voigt notation. The eigenstrain is
given by e0ðxÞ ¼ e0 1¼ 1 0 0 0 ; x2 X1; e0ðxÞ ¼ e0 2¼ 0 2 2 0 ; x2 X2: ðA:1Þ
Note that one can check that e0 1 and e
0
2 are
com-patible each other. So there is no internal stress across the boundary shared by oX1 and oX2. In
fact, the exact solution is
uðxÞ ¼ ðx1;0Þ; x2 X1; ð0:05; 4x1 0:2Þ; x2 X2; ðA:2Þ 0.1 m 0.05 m 0.05 m Ω1 ε0 1 Ω2 ε0 2 cut 0.1 m 0.1m 0.02 m 0.02 m Ω1 Ω (a) (b)
Fig. 11. (a) A patch test for strain compatibility. Note that the constrained boundary condition u2¼ 0 is enforced on the left-hand side of the bottom and is represented by circles. (b) A model problem to test the accuracy of the modified BEM formalism.
which results in r 0. The computational result also shows that the computed stress r is extremely close to zero which agrees with the exact solution. Finally, note that the total boundary elements used in the simulation are N1¼ 6, N2ð1Þ¼ 4 and
N2ð2Þ ¼ 4 where N1 is the number of boundary
ele-ments on the whole boundary, N2ð1Þand N ð2Þ 2 are the
number of boundary elements used to describe the boundariesoX1 andoX2(see (3.20) and (3.21)).
The next example with geometry and boundary conditions is shown in Fig. 11(b). The plane stress elastic constants used here are the same as those given in the previous example. The eigenstrain is given by e0ðxÞ ¼ 3:9 104 0 0 3:9 104 ; x2 X1; 0; x2 X n X1: 8 < : ðA:3Þ The total line elements used in the simulation are N1¼ 80 and N2ð1Þ¼ 16 on the boundary oX and
oX1, respectively. The stresses computed along the
line cut at 1/4 of the sideline of the body are shown in Fig. 12. We also use the commercial NA-STRAN FEM code to check the accuracy of the BEM formalism given by (3.20) and (3.21). Note that we have used the eight-node square mesh in the FEM simulation and the total number of meshes is 400. From Fig. 12 we find that the numerical results obtained by BEM simula-tion with few line elements on the boundaries agree well with those obtained by FEM simula-tion.
References
Aharoni, A., 1996. Introduction to the Theory of Ferromag-netism. Clarendon Press, Oxford.
Aharoni, A., Jakubovics, J.P., 1986. Cylindrical magnetic domains in small ferromagnetic spheres with uniaxial anisotropy. Philos. Mag. B 53, 133–145.
Berkov, D.V., Ramst€ock, K., Hubert, A., 1993. Solving micromagnetic problems. Phys. Stat. Sol. (a) 137, 207–225. Bhattacharya, K., 1993. Comparison of the geometrically nonlinear and linear theories of Martensitic transformation. Continuum Mech. Therm. 5, 205–242.
Brown, W.F., 1962. Magnetostatic Principles in Ferromagne-tism. In: Wohlfarth, E.P. (Ed.), Selected Topics in Solid State Physics, vol. 1. North-Holland, Amsterdam. Brown, W.F., 1963. Micromagnetics. Wiley, New York. Brown, W.F., 1966. Magnetoelastic Interactions. Springer,
Berlin.
Callegaro, L., Puppin, E., 1996. Stress dependence of coercivity in Ni films: thin films to bulk transition. Appl. Phys. Lett. 68, 1279–1281.
Callegaro, L., Puppin, E., 1997. Rotational hysteresis model for stressed ferromagnetic films. IEEE Trans. Magn. 33, 1007– 1011.
Clark, A.E., 1980. In: Wohlfarth, E.P. (Ed.), Ferromagnetic Materials, vol. 1. North-Holland, Amsterdam (Chapter 7). Clark, A.E., 1986. Magnetostriction in twinned [1 1 2] crystals of Tb0:27Dy0:73Fe2. IEEE Trans. Magn. Mag-22, 973–975. Crew, D.C., Lewis, L.H., 2001. Effect of grain alignment on
magnetic structure in nanoscale material. IEEE Trans. Magn. 37, 2512–2514.
Cullen, J.R., Hathaway, K.B., Clark, A.E., 1997. Critical behavior of cubic magnetostrictive materials under stress. J. Appl. Phys. 81, 5417–5419.
Cullity, B.D., 1972. Introduction to Magnetic Materials. Addison-Wesley, Reading, MA.
Deng, H., Jarratt, J.D., Minor, M.K., Barnard, J.A., 1997. Artificially controlled stress anisotropy and magnetic properties of FeTaN thin films. J. Appl. Phys. 81, 4510– 4512. 0.02 0.04 0.06 0.08 0.1 0 2 4 6 8 FEM BEM FEM BEM FEM BEM 0.02 0.04 0.06 0.08 0.1 -2 -1 0 1 2 0.02 0.04 0.06 0.08 0.1 -4 -2 0 2 4 x1 x1 x1
σ
11 (MP a)σ
12 (MP a)σ
22 (MP a) (b) (a) (c)Fig. 12. Comparison of the modified BEM formulation with commercial NASTRAN FEM code for stresses computed along the line cut at 1/4 of the sideline of the body shown in Fig. 11(b).
