1-s2.0-S0022509624002333-main

J. Mech. Phys.
Solids 191 (2024) 105767
Contents lists available at ScienceDirect
Journal of the Mechanics and Physics of Solids

journal homepage: www.elsevier.com/locate/jmps
A micromagnetic-mechanically coupled phase-field model for

fracture and fatigue of magnetostrictive alloys
Shen Sun a ,1 , Qihua Gong a ,1 , Yong Ni b ,∗, Min Yi a ,∗
a State Key Laboratory of Mechanics and Control for Aerospace Structures & Key Laboratory for Intelligent Nano Materials and Devices of
Ministry of Education & Institute for Frontier Science & College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics
(NUAA), Nanjing 210016, China
b
CAS Key Laboratory of Mechanical Behavior and Design of Materials, Department of Modern Mechanics, CAS Center for Excellence in
Complex System Mechanics, University of Science and Technology of China, Hefei 230026, China
ARTICLE INFO ABSTRACT
Keywords: Magnetostrictive alloys are usually brittle materials with micromagnetic structures. Their
Micromagnetic-mechanical coupling structural reliability and durability depend on the complex micromagnetic-mechanical coupling
Phase-field model at smaller length scales encompassing the evolution of micromagnetic structures. Herein
Fracture and fatigue
we propose a micromagnetic-mechanically coupled phase-field model for fracture and fa-
Crack propagation
tigue behavior of magnetostrictive alloys with evolution of the micromagnetic structure. The
Micromagnetic domain
Magnetostrictive alloys
thermodynamically-consistent model is derived from microforce theory, laws of thermodynam-
ics, and Coleman–Noll analysis. The evolution of crack phase-field and magnetization-vector
order parameters that are fully coupled is governed by history field dependent Allen–Cahn and
Landau–Lifshitz–Gilbert equations, respectively. The model is extended to fatigue by introducing
a degradation prefactor for the fracture energy as a function of positive elastic energy. One-
dimensional analyses are then presented to anatomize the crack driving forces in terms of fully
coupled micromagnetic-mechanical and pure mechanical driving force. We demonstrate the
model capabilities by finite-element numerical studies on the micromagnetic domain evolution
during the crack propagation and the influence of external magnetic field for type-I, type-
II, and three-point bending fracture, as well as for the fracture of a single-edge notched
specimen with an elliptical inclusion. The simulation result shows that depending on how
micromagnetic domains are switched under micromagnetic-mechanical coupling, the magnetic
field can enhance or decrease the critical load. In the presence of inclusion with larger fracture
toughness, a crack is found to nucleate in the tri-junction of multi-domain micromagnetic
structure owing to the high elastic strain around the tri-junction point. It is further found
that a suitable magnetic field promoting magnetization-vector rotation around the crack tip
could remarkably improve the fracturing load and fatigue life. The results demonstrate the
model promising for the study of micromagnetic-mechanically coupled fracture and fatigue in
magnetostrictive alloys.
1. Introduction
Magnetostrictive alloys have received considerable attention in the field of sensors, actuators, and energy harvesting de-
vices (Eerenstein et al., 2006; Lei et al., 2013; Beveridge et al., 2011; Aguilar-Arteaga et al., 2010; Jiles, 2003). Magnetostrictive
∗ Corresponding authors.
E-mail addresses: yni@ustc.edu.cn (Y. Ni), yimin@nuaa.edu.cn (M. Yi).
1 Authors contributed equally.
https://doi.org/10.1016/j.jmps.2024.105767
Received 2 April 2024; Received in revised form 10 June 2024; Accepted 3 July 2024
Available online 6 July 2024
0022-5096/© 2024 Elsevier Ltd. All rights are reserved, including those for text and data mining, AI training, and similar technologies.
S. Sun et al. Journal of the Mechanics and Physics of Solids 191 (2024) 105767
alloys are characterized by the spontaneous magnetization and magnetic domain structure at the micromagnetic scale. The
magnetization vectors in the domains will be parallel to the external magnetic field if a strong magnetic field is applied. With
a composition near Tb0.27 Dy0.73 Fe2 , Terfenol-D is one of the leading magnetostrictive alloys which can sensitively response to
the changes of external magnetic field. It is brittle, but displays low magnetocrystalline anisotropy and high magnetostriction at
ambient temperature (Clark, 1980). The macroscopic characteristics of magnetostrictive alloys are the results of the micromagnetic
structures and their interaction with the external field. To clarify the macroscopic performances of magnetostrictive alloys, it is
necessary to explore the evolution of micromagnetic structure under external fields. Landau and Lifshitz (1935) first reported
the antiparallel magnetic domain structure using the continuum physics description of matter. Then Brown (1963) proposed the
micromagnetism theory, which is based on the magnetization vector that conquers the limit of the magnetic domain theory. Gilbert
(2004) reformulated the Landau–Lifshitz equation and then applied the famous Landau–Lifshitz–Gilbert (LLG) equation to describe
the temporal evolution of magnetization vector.
Magneto-mechanical properties are sophisticated when magnetostrictive alloys are placed in the electromagnetic field. Since
1960s, macroscopic phenomenological constitutive models for magnetostrictive alloys have been initiated. Based on enormous
multiphysics experiments, phenomenological magneto-mechanical coupling models at the macroscale were derived, including the
linear theory of Pao and Yeh (1973) under the constant static magnetic field and nonlinear magneto-elastic coupling models (Zhou
and Zheng, 1997; Zheng and Liu, 2005; Zhou et al., 2008, 2009). On the other aspect, the magnetic domain switching, inhomogeneity
and hysteresis were theoretically analyzed (Maugin et al., 1992; Maugin, 1995). Nevertheless, these models are difficult to explain
the mutual coupling between the micromagnetic structure and the mechanical response.
Since magnetostrictive alloys are usually brittle, the failure criterion of this kind of materials still needs accurate interpretations.
Scholars usually combine theory and experiment to promote their research. Shindo (1977, 1978, 1980) not only solved the stress
fields around the crack tip considering different crack forms, but also focused on the changes of magnetic fields that are influenced
by the dynamic stress concentration and the crack (Shindo et al., 1997; Shindo and Tohyama, 1998; Shindo et al., 1999). Liang
et al. (2000) derived the complex function solutions for plane isotropic magnetoelasticity and plane cracks from the linear theory
of Pao and Yeh (1973), and then proposed the fracture criteria. So far, the magnetoelastic fracture problem has been solved under
various boundary conditions. The path independent integral and energy release rate of ferromagnetic crack problem based on the
extension of fracture mechanics method is accessible as well.
However, whether the external magnetic fields will cause the variations in the fracture properties remains to be discussed.
The properties may be connected to the geometry of the crack. On the one hand, soft ferromagnetic materials have a very small
magnetostriction. As for manganese–zinc ferrite ceramics, the magnetic field has no obvious effects on fracturing load (Fang et al.,
2004). On the other hand, Shindo et al. (2004, 2005, 2006, 2008) demonstrated that nickel–iron soft ferromagnetic alloys in
magnetic field fail at a lower load compared to the tests without magnetic field. Similarly, Colussi et al. (2016) combined experiments
and simulations to study the precracked Terfenol-D specimens and demonstrated that the presence of magnetic field decreases the
fracturing load.
The above fracture performance and its dependence on external magnetic field are intrinsically related to the micromagnetic-
mechanical coupling. Nevertheless, both the micromagnetic domain structure evolution and the crack topology are incisively
challenging topics for magnetostrictive alloys. Phase-field method is supposed to handle the two key problems appropriately. On the
one hand, phase-field models usually use the magnetization vector as the order parameters, in which micromagnetics phase-field
simulation quantitatively predicted the domain structure, magnetization switching and hysteresis processes (Fidler and Schrefl,
2000; Fischbacher et al., 2007). Zhu et al. (2001) introduced the magnetoelastic energy term into the LLG equation so that the
stress field can be included in the model. Landis (2008) modeled the evolution of magnetic domain and martensite structures
in ferromagnetic shape memory alloys with the application of microforce theory. Based on the time dependent Ginzburg–Landau
(TDGL) equation, Wang and Zhang (2013) proposed a real-space phase-field model suitable for the microstructure evolution of
ferromagnetic materials, in which the strain and the magnetization vector were taken as the independent variables. Yi and Xu
(2014) proposed an unconstrained phase-field model to simulate the magnetic domain evolution in ferromagnetic materials.
On the other hand, the phase-field method is emerging as an efficient route for tracking the crack topology. Francfort and
Marigo (1998) put forward the phase-field idea for quasi-static brittle fracture through the study of Griffith (1920) fracture model
and variational formula. Amor et al. (2009) then proposed the volumetric–deviatoric decomposition of the elastic strain tensor
to distinguish the tensile part and the compressive part for the elastic energy density, while Miehe et al. (2010a) elaborated the
spectral split for the brittle fracture. Steinke and Kaliske (2019) proposed an alternative way to decompose the stress with respect
to the crack orientation, in which the convergence is better but the determination of the crack orientation is intricate. Then the
ductile fracture phase-field model was proposed in which a plastic part of strain energy density was added to the stored energy
functional (Ambati et al., 2015; Borden et al., 2016; Alessi et al., 2018a). Due to the significance of fatigue problems in engineering
practice, researches on the prediction of fatigue properties have been conducted widely (Mei et al., 2020; Dang et al., 2022; Yi et al.,
2024). The phase-field fracture model is further extended to fatigue (Alessi et al., 2018b; Carrara et al., 2020; Seiler et al., 2020;
Schreiber et al., 2020; Lo et al., 2019; Tang et al., 2024). With the wide application of multi-functional materials, the multi-field
fracture problems emerge. Bourdin et al. (2014) and Miehe et al. (2015b,a) connected the fracture with the thermal effect. Sridhar
and Keip (2019) focused on the piezoelectric effect during the fracture progress. Xu et al. (2010) and Abdollahi and Arias (2011)
simulated the failure in a ferroelectric materials. Shi et al. (2017) analyzed the domain-switching in a central cracked ferromagnetic
plate. Phase-field models have been applied to the fracture of soft elastomers which contain hard magnetics (Moreno-Mateos et al.,
2023) or piezoelectric particles (Moreno-Mateos et al., 2024) as well. Confirmed by experiments, the models have been proved
reliable in simulating the enhancement of fracture property by magnetic or piezoelectric particles. Moreno-Mateos and Steinmann
2
Fig. 1. Sharp and diffusive crack topology. (a) Sharp crack surface 𝛤 embedded in the solid 𝛺. (b) Diffusive crack surface 𝛤𝑙 (𝑐) and the micromagnetic domain
structure.
(2024) also provided a possible idea for the combination of configurational force method and phase-field model. However, these
models use the macroscopic constitutive relations to include the multi-field coupling, and thus are powerless for ascertaining the
microstructure evolution (e.g., ferroelectric domain in ferroelectrics, micromagnetic domain in ferromagnetics, etc.) during the crack
propagation or fracture process.
In order to clarify the role of microstructure evolution in the fracture behavior of magnetostrictive alloys, herein we propose
a micromagnetic-mechanically coupled phase-field fracture and fatigue model that resolves the micromagnetic domain structure
and its temporal/spatial evolution, as well as the crack. Specifically, the magnetization vector and crack phase field are utilized as
order parameters to describe the micromagnetic domain structure and crack topology, respectively. In contrast to the macroscopic
constitutive relations for the multi-field coupling, constitutive relations derived from a micromagnetic-mechanically coupled free
energy functional based on micromagnetic theory is adopted. The model is thermodynamical consistently derived from microforce
theory, laws of thermodynamics, and Coleman–Noll analysis. The model capabilities are demonstrated by finite-element numerical
simulations on the micromagnetic domain evolution during the crack propagation and the influence of external magnetic field on
fracture behavior of magnetostrictive alloys.
This article is organized as follows. In Section 2, the fully micromagnetic-mechanically coupled phase-field fracture model for
ferromagnetic materials is presented, with a focus on the thermodynamical consistent model derivation, evolution equations for
micromagentic structure and crack phase field, and free energy density. Section 3 presents the extension of model in Section 2 to
fatigue. Section 4 briefly illustrates the finite-element implementation of the model. Section 5 discusses the crack driving force in a
one-dimensional analysis. Then several numerical examples on fracture and fatigue of magnetostrictive alloys are shown in Section 6
to demonstrate the model capability. Section 7 summarizes the work.
2. Micromagnetic-mechanically coupled phase-field model for fracture
For exploring the fracture behavior of magnetostrictive alloys, both the micromagnetic domain evolution and crack propagation
have to be simultaneously considered. Considering a region 𝛺 ⊂ 𝐑𝑛 , 𝑛 = 1, 2, 3 with the boundary 𝜕𝛺, there is a crack in the body,
as shown in Fig. 1. The structure is divided into different micromagnetic domains which are presented by the magnetization vector
𝑀𝑖 that could be the phase-field order parameter. The vector is represented by the product of the constant saturation magnetization
magnitude 𝑀s and the magnetization unit vector 𝑚𝑖 at a constant temperature which is far below the Curie point, i.e.,
𝑀𝑖 = 𝑀s 𝑚𝑖 . (1)
The components of the unit vector describe the spatial direction of magnetization. Another non-conservative phase-field order
parameter 𝑐 (𝐱, 𝑡) ∈ [0, 1] is defined to describe the crack state of the material. The perfect state of the material is characterized
for 𝑐 = 0 while 𝑐 = 1 indicates that the material is broken. The numerical approximation of the topology of the crack in the
phase-field model is based on the second order fracture energy density, i.e.,
1 ( 2 )
𝛾= 𝑐 + 𝑙2 ‖∇𝑐‖2 (2)
2𝑙
in which 𝑙 is the length scale that governs the width of the diffusive crack. At first it was a numerical parameter used to regularize
the sharp crack in Fig. 1(a). However, it is now widely assumed that 𝑙 is a material parameter (Amor et al., 2009). It was
experimentally determined by Nguyen et al. (2016) using the critical stress in uniformly stretched samples. The surface energy
for crack is approximated by
𝜓 𝑐 (𝛤 ) = 𝐺c d𝑆 ≈ 𝐺c 𝛾d𝑉 (3)
∫𝛤 ∫𝛺
where 𝐺c is the critical fracture energy density.
3
2.1. Balance law
• Balance of linear momentum

The momentum of a system is constant if no external forces are acting on the system. For the body 𝛺 with a boundary 𝜕𝛺,
the quasi-static balance equation can be given as
𝜎𝑖𝑗,𝑗 + 𝜌𝑏𝑖 = 0 in 𝛺 , (4)
𝑢𝑖 = 𝑢̂ 𝑖 on 𝜕𝛺𝑢 , 𝜎𝑗𝑖 𝑛𝑗 = 𝑡̂𝑖 on 𝜕𝛺𝜎 (5)
where 𝜎𝑖𝑗 is the second-order stress tensor and 𝑏𝑖 is the body force of unit mass 𝜌. If the body force is not considered, it is
assumed that 𝑏𝑖 = 0. The Latin indices (𝑖, 𝑗) run over the range of 1–3. 𝑢̂ 𝑖 is the displacement prescribed on the boundary 𝜕𝛺𝑢 .
𝑛𝑗 is the outward unit normal vector and 𝑡̂𝑖 is the traction on the surface 𝜕𝛺𝜎 .
• Balance of angular momentum
The stress tensor is symmetric, i.e.,
𝜎𝑖𝑗 = 𝜎𝑗𝑖 . (6)
• Gauss’s law for magnetic flux

As for the ferromagnetic materials, electromagnetic motions are governed by the Maxwell’s equations. Here come the magnetic
field, 𝐻𝑖 , magnetic induction, 𝐵𝑖 , volume current density, 𝐽𝑖 , surface current density, 𝐾𝑖 , and the magnetic vector potential,
𝐴𝑖 in this work. The Gauss’s law states that no magnetic monopole exists and that the total magnetic flux 𝐵𝑖 through a closed
surface must be zero. Specifically, in any arbitrary volume
𝐵𝑖,𝑖 = 0. (7)
Faraday’s law describes how a time-varying magnetic field creates or induces an electric field. In integral form, it states that
the work per unit charge required to move a charge around a closed loop equals the rate of change of the magnetic flux
through the enclosed surface, i.e.,
𝜖𝑖𝑗𝑘 𝐸𝑘,𝑗 = −𝐵̇ 𝑖 (8)
where 𝜖𝑖𝑗𝑘 is the permutation tensor. It comes that 𝜖𝑖𝑗𝑘 = 1 if 𝑖𝑗𝑘 = 123, 231, 312, 𝜖𝑖𝑗𝑘 = −1 if 𝑖𝑗𝑘 = 321, 213, 132 and 𝜖𝑖𝑗𝑘 = 0
under all other circumstances.
Ampėre’s circuital law relates the integrated magnetic field around a closed loop to the electric current passing through the
loop. Considering a case of magnetostatics, the current is steady. The line integral of the magnetic induction around any closed
curve is proportional to the conduction current through any surface enclosed by that curve. The formulation is
𝜖𝑖𝑗𝑘 𝐻𝑘,𝑗 = 𝐽𝑖 . (9)
The magnetic boundary conditions are described as
𝐵𝑖 𝑛𝑖 = 𝐵̂ on 𝜕𝛺𝐵 , 𝜙 = 𝜙̂ on 𝜕𝛺𝜙 , 𝜖𝑖𝑗𝑘 𝐻𝑗 𝑛𝑘 = 𝐾̂ 𝑖 on 𝜕𝛺𝐾 (10)
in which 𝐵̂ 𝑖 is the specified value on the boundary part 𝜕𝛺𝐵 , 𝜙 is the scalar magnetic potential, and 𝜙̂ is the prescribed potential
on the boundary part 𝜕𝛺𝜙 . 𝐾̂ is the specified current density on the boundary part 𝜕𝛺𝐾 . The magnetic field 𝐻𝑖 is calculated
by the negative gradient of 𝜙, i.e.,
𝐻𝑖 = −𝜙,𝑖 . (11)
• Balance of microforce associated with magnetization vector 𝑀𝑖

Since magnetization vector 𝑀𝑖 is critical for describing the micromagnetic structures, the free energy will be required to depend
on 𝑀𝑖 that could be taken as a configurational quantity. Thus, following Gurtin’s microforce theory (Gurtin, 1996) and Landis’s
work (Landis, 2008), there exist a set of forces that are work conjugate to 𝑀𝑖 . These forces are named as microforces because
they are relevant to the local micromagnetic structure rather than the macroscopic movements. Here, the microforces are
associated with order parameters 𝑀𝑖 and are defined as follows.
𝜉𝑗𝑖 : the microforce stress tensor,
𝜉𝑗𝑖 𝑛𝑗 : the surface microforce with 𝑛𝑗 as the unit outer normal to 𝜕𝛺,
𝜋𝑖 : the internal microforce vector,
𝑓𝑖 : the external microforce vector.
Referring to the rotation nature of the magnetization (Gilbert, 2004), the angular momentum balance in this work can be
given as
( )
1 1 ̇
𝜖𝑖𝑗𝑘 𝑀𝑗 𝜉𝑙𝑘 𝑛𝑙 d𝑆 + 𝜖𝑖𝑗𝑘 𝑀𝑗 𝜋𝑘 d𝑉 + 𝜖𝑖𝑗𝑘 𝑀𝑗 𝑓𝑘 d𝑉 = 𝑀 d𝑉 (12)
𝜇0 ∫𝜕𝛺 ∫𝛺 ∫𝛺 ∫𝛺 𝛾0 𝑖
4
where 𝛾0 is the gyromagnetic ratio for an electron spin. The balance law must hold for any arbitrary volume. With the use of
Gauss’s law, one yields
( ) 𝜇
𝜖𝑖𝑗𝑘 𝑀𝑗 𝜉𝑙𝑘,𝑙 + 𝜋𝑘 + 𝑓𝑘 + 𝜖𝑖𝑗𝑘 𝑀𝑗,𝑙 𝜉𝑙𝑘 = 0 𝑀̇ 𝑖 . (13)
𝛾0
It will be shown Eq. (13) leads to the Landau–Lifshitz–Gilbert (LLG) equation of micromagnetics after performing the analysis
of the second law of thermodynamics. The magnetization only changes its direction, so the microforce as the driving force on
the change of magnetization should lie perpendicular to the direction of the magnetization. Therefore, the microforces along
the direction of the magnetization should vanish, i.e.,
( )
𝑀𝑘 𝜉𝑙𝑘,𝑙 + 𝜋𝑘 + 𝑓𝑘 d𝑉 = 0. (14)
∫𝛺
For the volume arbitrariness, it implies
( )
𝑀𝑘 𝜉𝑙𝑘,𝑙 + 𝜋𝑘 + 𝑓𝑘 = 0. (15)
From Eqs. (13) and (15), one yields

𝜇0 𝜖𝑟𝑖𝑠 𝜖𝑖𝑗𝑘
𝜉𝑙𝑟,𝑙 + 𝜋𝑟 + 𝑓𝑟 = 𝜖𝑟𝑖𝑠 𝑀̇ 𝑖 𝑀𝑠 − 𝑀𝑗,𝑙 𝜉𝑙𝑘 𝑀𝑠 . (16)
𝛾0 𝑀 2 𝑀2
• Balance of microforce associated with crack phase field 𝑐

The crack phase field is to describe the diffusive fracture topology. The definition for the microforce balance of the crack phase
field is expressed as:
𝜉𝑖𝑐 : the microforce stress,
𝜉𝑖𝑐 𝑛𝑖 : the surface microforce with 𝑛𝑖 as the unit outer normal to 𝜕𝛺,
𝜋 𝑐 : the internal microforce,
𝑓 𝑐 : the external microforce.
Then the balance equation is written as
𝜉𝑖𝑐 𝑛𝑖 d𝑆 + 𝜋 𝑐 d𝑉 + 𝑓 𝑐 d𝑉 = 0. (17)
∫𝜕𝛺 ∫𝛺 ∫𝛺
By means of the Gauss law, it comes
𝑐
𝜉𝑖,𝑖 + 𝜋 𝑐 + 𝑓 𝑐 = 0. (18)
• Balance of energy
When it comes to the balance law of energy, only the quasistatic and isothermal problem here is considered for simplicity.
The materials time derivative of the total energy is equal to the sum of the power due to external forces and the heat input to
the system. The system is in contact with an outside thermal reservoir and the change in the system is slow enough such that
the kinetic energy can be ignored and the temperature is always same as that of the thermal reservoir. For magnetostrictive
alloys, there is electromagnetic energy flux represented by the Poynting vector 𝐒 = 𝐄 × 𝐇 = 𝜖𝑖𝑗𝑘 𝐸𝑗 𝐻𝑘 . The electric field also
performs work on the electric current inside the body. The first law of thermodynamics can be presented as
d ( ) ( )
𝑒d𝑉 = 𝑒d𝑉
̇ = 𝑢̇ 𝜎 − 𝜖𝑖𝑗𝑘 𝐸𝑗 𝐻𝑘 + 𝜉𝑖𝑐 𝑐̇ + 𝜉𝑖𝑗 𝑀̇ 𝑗 𝑛𝑖 d𝑆 + −𝐸𝑖 𝐽𝑖 + 𝑓 𝑐 𝑐̇ + 𝑓𝑖 𝑀̇ 𝑖 d𝑉 . (19)
d𝑡 ∫𝛺 ∫𝛺 ∫𝜕𝛺 𝑗 𝑖𝑗 ∫𝛺
Using Gauss’ flux theorem to simplify the area integral, the equation can be transferred into
( )
𝑐
𝑒d𝑉
̇ = 𝑢̇ 𝜎 + 𝑢̇ 𝑗 𝜎𝑖𝑗,𝑖 − 𝜖𝑖𝑗𝑘 𝐸𝑗,𝑖 𝐻𝑘 − 𝜖𝑖𝑗𝑘 𝐸𝑗 𝐻𝑘,𝑖 + 𝜉𝑖,𝑖 𝑐̇ + 𝜉𝑖𝑐 𝑐̇ ,𝑖 + 𝜉𝑖𝑗,𝑖 𝑀̇ 𝑗 + 𝜉𝑖𝑗 𝑀̇ 𝑗,𝑖 − 𝐸𝑖 𝐽𝑖 + 𝑓 𝑐 𝑐̇ + 𝑓𝑖 𝑀̇ 𝑖 d𝑉 . (20)
∫𝛺 ∫𝛺 𝑗,𝑖 𝑖𝑗
Considering the balance law in Eqs. (4), (8), (9) and (18), one yields
𝑒̇ = 𝑢̇ 𝑗,𝑖 𝜎𝑖𝑗 + 𝐵̇ 𝑖 𝐻𝑖 − 𝜋 𝑐 𝑐̇ + 𝜉𝑖𝑐 𝑐̇ ,𝑖 + 𝜉𝑖𝑗,𝑖 𝑀̇ 𝑗 + 𝜉𝑖𝑗 𝑀̇ 𝑗,𝑖 + 𝑓𝑖 𝑀̇ 𝑖 . (21)
2.2. Constitutive relations and evolution equations
• Second law of thermodynamics

For the isothermal process, the second law of thermodynamics is simplified as the non-negative increase of the entropy 𝑠. The
Clausius–Duhem (dissipation) inequality is
𝑠d𝑉
̇ ≥ 0. (22)
∫𝛺
• Free energy imbalance

The Helmholtz free energy density 𝜓 is defined as
𝜓 = 𝑒 − 𝑇𝑠 (23)
5
where 𝑇 is the temperature. Taking the time derivatives at both sides, it comes
𝜕𝜓 𝜕𝑒 𝜕𝑠 𝜕𝑇
= −𝑇 − 𝑠 (24)
𝜕𝑡 𝜕𝑡 𝜕𝑡 𝜕𝑡
Due to the isothermal condition we consider here, the temperature remains constant. Taking into account Eq. (21) and the
thermodynamics inequality (22), one has
𝜕𝜓
≤ 𝑢̇ 𝑗,𝑖 𝜎𝑖𝑗 + 𝐵̇ 𝑖 𝐻𝑖 − 𝜋 𝑐 𝑐̇ + 𝜉𝑖𝑐 𝑐̇ ,𝑖 + 𝜉𝑖𝑗,𝑖 𝑀̇ 𝑗 + 𝜉𝑖𝑗 𝑀̇ 𝑗,𝑖 + 𝑓𝑖 𝑀̇ 𝑖 . (25)
𝜕𝑡
• Coleman–Noll type analysis
To get the thermodynamically consistent constitutive relations, the Coleman–Noll argument (Noll and Coleman, 1974) is
applied. Invoking Truesdell’s principle of equipresence (Truesdell and Noll, 2004), it is reasonable to assume that 𝜓, 𝑠, 𝑒,
𝜎𝑖𝑗 , 𝜉𝑗𝑖 and 𝜉𝑖𝑐 are functions depending on 𝜀𝑖𝑗 , 𝐵𝑖 , 𝑀𝑖 , 𝑀𝑖,𝑗 , 𝑀̇ 𝑖 , 𝑐, 𝑐,𝑖 , 𝑐.̇ In this model, the Helmholtz free energy takes the
following form
( )
𝜓 = 𝜓 𝜀𝑖𝑗 , 𝐵𝑖 , 𝑀𝑖 , 𝑀𝑖,𝑗 , 𝑀̇ 𝑖 , 𝑐, 𝑐,𝑖 , 𝑐̇ . (26)
The chain rule yields the time derivative of 𝜓 as
𝜕𝜓 𝜕𝜓 𝜕𝜀𝑖𝑗 𝜕𝜓 𝜕𝐵𝑖 𝜕𝜓 𝜕𝑀𝑖 𝜕𝜓 𝜕𝑀𝑖,𝑗 𝜕𝜓 𝜕 𝑀̇ 𝑖 𝜕𝜓 𝜕𝑐 𝜕𝜓 𝜕𝑐,𝑖 𝜕𝜓 𝜕 𝑐̇
= + + + + + + + . (27)
𝜕𝑡 𝜕𝜀𝑖𝑗 𝜕𝑡 𝜕𝐵𝑖 𝜕𝑡 𝜕𝑀𝑖 𝜕𝑡 𝜕𝑀𝑖,𝑗 𝜕𝑡 ̇
𝜕 𝑀𝑖 𝜕𝑡 𝜕𝑐 𝜕𝑡 𝜕𝑐,𝑖 𝜕𝑡 𝜕 𝑐̇ 𝜕𝑡
The magnetostrictive alloys are usually elastic, brittle and of small deformation. Therefore, the small strain theory is adopted,
i.e.,
1( )
𝜀𝑖𝑗 = 𝑢 + 𝑢𝑗,𝑖 . (28)
2 𝑖,𝑗
Substituting the chain rule, Eqs. (6) and (28) into the inequality (25), analysis gives
( ) ( ) ( ) ( )
𝜕𝜓 𝜕𝜓 𝜕𝜓 𝜕𝜓
− 𝜎𝑖𝑗 𝜀̇ 𝑖𝑗 + − 𝐻𝑖 𝐵̇ 𝑖 + − 𝜉𝑗𝑖 𝑀̇ 𝑖,𝑗 + + 𝜋 𝑐 𝑐̇
𝜕𝜀𝑖𝑗 𝜕𝐵𝑖 𝜕𝑀𝑖,𝑗 𝜕𝑐
( ) ( ) (29)
𝜕𝜓 𝜕𝜓 𝜕𝜓 𝜕𝜓
+ − 𝑓𝑖 − 𝜉𝑗𝑖,𝑗 𝑀̇ 𝑖 + − 𝜉𝑖𝑐 𝑐̇ ,𝑖 + 𝑀̈𝑖+ 𝑐̈ ≤ 0.
𝜕𝑀𝑖 𝜕𝑐,𝑖 𝜕 𝑀̇ 𝑖 𝜕 𝑐̇
Using Eq. (16), one further has
( ) ( ) ( ) ( )
𝜕𝜓 𝜕𝜓 𝜕𝜓 𝜕𝜓
− 𝜎𝑖𝑗 𝜀̇ 𝑖𝑗 + − 𝐻𝑖 𝐵̇ 𝑖 + − 𝜉𝑗𝑖 𝑀̇ 𝑖,𝑗 + + 𝜋 𝑐 𝑐̇
𝜕𝜀𝑖𝑗 𝜕𝐵𝑖 𝜕𝑀𝑖,𝑗 𝜕𝑐
( ) ( ) (30)
𝜕𝜓 𝜖 𝜖 𝜕𝜓 𝜕𝜓 ̈
+ + 𝜋𝑖 + 𝑖𝑟𝑠 𝑟𝑙𝑘 𝑀𝑙,𝑗 𝜉𝑗𝑘 𝑀𝑠 𝑀̇ 𝑖 + − 𝜉𝑖𝑐 𝑐̇ ,𝑖 + 𝑀𝑖 ≤ 0.
𝜕𝑀𝑖 𝑀2 𝜕𝑐,𝑖 𝜕 𝑀̇ 𝑖
Due to the principle of equipresence (Noll and Coleman, 1974), the coefficients of linear terms in Eq. (30) must vanish
considering the arbitrariness of the thermodynamics process. Then one has
𝜕𝜓 𝜕𝜓 𝜕𝜓 𝜕𝜓 𝜕𝜓 𝜕𝜓
= 𝜎𝑖𝑗 , = 𝐻𝑖 , = 𝜉𝑗𝑖 , = 𝜉𝑖𝑐 , = 0, = 0. (31)
𝜕𝜀𝑖𝑗 𝜕𝐵𝑖 𝜕𝑀𝑖,𝑗 𝜕𝑐,𝑖 𝜕 𝑀̇ 𝑖 𝜕 𝑐̇
Therefore, it can be concluded that the free energy density 𝜓 is independent of 𝑀̇ 𝑖 and 𝑐.̇ Then 𝜓 is simplified as
( )
𝜓 = 𝜓 𝜀𝑖𝑗 , 𝐵𝑖 , 𝑀𝑖 , 𝑀𝑖,𝑗 , 𝑐, 𝑐,𝑖 . (32)
The nonlinear terms are reduced to
( ) ( )
𝜕𝜓 𝜕𝜓 𝜖𝑖𝑗𝑘 𝜖𝑗𝑟𝑠
+ 𝜋 𝑐 𝑐̇ + + 𝜋𝑖 + 𝑀𝑟,𝑙 𝜉𝑙𝑠 𝑀𝑘 𝑀̇ 𝑖 ≤ 0. (33)
𝜕𝑐 𝜕𝑀𝑖 𝑀2
• Evolution equations
The inequality (33) can be satisfied by the following assumptions, i.e.,
𝜕𝜓
𝜋𝑐 = − − 𝜂 𝑐,̇ (34)
𝜕𝑐
𝜕𝜓 𝜖𝑖𝑗𝑘 𝜖𝑗𝑟𝑠
𝜋𝑖 = − − 𝑀𝑟,𝑙 𝜉𝑙𝑠 𝑀𝑘 − 𝛽𝑖𝑗 𝑀̇ 𝑗 , (35)
𝜕𝑀𝑖 𝑀2
with a non-negative mobility constant 𝜂 and a viscosity tensor 𝛽𝑖𝑗 which control the evolution rate of order parameters.
Combining Eqs. (18), (31) and (34), the governing equation for the fracture order parameter 𝑐 can be derived as
( )
𝜕𝜓 𝜕𝜓
− − 𝜂 𝑐̇ + 𝑓 𝑐 = 0. (36)
𝜕𝑐,𝑖 ,𝑖 𝜕𝑐
Similarly, the magnetic part of the above theory is equivalent to the Landau–Lifshitz–Gilbert equation. The viscosity tensor
must satisfy
𝛽𝑖𝑗 𝑀̇ 𝑗 𝑀̇ 𝑖 ≥ 0. (37)
6
A governing equation for the evolution of the magnetization vector reads

( )
𝜕𝜓 𝜕𝜓 𝜇0
− + 𝑓𝑖 = 𝛽𝑖𝑗 𝑀̇ 𝑗 + 𝜖𝑖𝑗𝑘 𝑀̇ 𝑗 𝑀𝑘 . (38)
𝜕𝑀𝑖,𝑗 ,𝑗 𝜕𝑀𝑖 𝛾0 𝑀 2
The viscosity parameter 𝛽𝑖𝑗 is related to damping parameter 𝛼, i.e.,
𝜇0 𝛼
𝛽𝑖𝑗 = 𝛿 . (39)
𝛾0 𝑀s 𝑖𝑗
Then Eq. (38) can be written as
( )
𝜕𝜓 𝜕𝜓 𝜇 𝛼 𝜇0
− + 𝑓𝑖 = 0 𝑀̇ 𝑖 + 𝜖𝑖𝑗𝑘 𝑀̇ 𝑗 𝑀𝑘 . (40)
𝜕𝑀𝑖,𝑗 ,𝑗 𝜕𝑀𝑖 𝛾0 𝑀s 𝛾0 𝑀 2
2.3. Free energy density
After Legendre transformation, the magnetic enthalpy ℎ = 𝜓 − 𝐵𝑖 𝐻𝑖 is used, which is a functional of 𝜀𝑖𝑗 , 𝐻𝑖 , 𝑀𝑖 , 𝑀𝑖,𝑗 , 𝑐 and 𝑐,𝑖 .
Accordingly, the constitutive relations in Eq. (31) and the governing equations in Eqs. (36) and (40) should be revised. The total
magnetic enthalpy can be written as
ℎ = ℎexc + ℎani + ℎent1 + ℎent2 + ℎe + ℎc . (41)
in which ℎexc is the magnetic exchange energy density, ℎani is the magnetocrystalline anisotropy energy density, ℎent1 is the
magnetization-dependent magnetostatic energy density, ℎent2 is the magnetic Maxwell energy describing the background magnetic
energy due to the non-zero magnetic permeability of the vacuum, ℎe is the elastic energy density, and ℎc is the crack energy density.
It is supposed that the magnetic energy may degrade in the crack phase, i.e., ℎ in the cracked phase is different from ℎ in
the intact phase. This is realized by a degradation function 𝑔 (𝑐, 𝜁), in which 𝜁 is an additional parameter to describe the extent
how the magnetic energy is degraded (i.e., the ratio of magnetization dependent magnetic energy in the crack phase to that in the
intact phase is 𝜁). The choice of 𝑔 (𝑐, 𝜁) will be discussed in the following. In detail, according to the micromagnetic theory (Brown,
1963; Hubert and Schäfer, 1998; Kronmüller and Fähnle, 2003), ℎexc can be expressed in terms of 𝑔 (𝑐, 𝜁) and the derivatives of
magnetization unit vector 𝐦 with respect to the coordinates, i.e.,
ℎexc = 𝑔 (𝑐, 𝜁) 𝐴e ‖∇𝐦‖2 = 𝑔 (𝑐, 𝜁) 𝐴e 𝑚𝑖,𝑗 𝑚𝑖,𝑗 (42)
where 𝐴e is the exchange stiffness constant. Similarly, ℎani for a cubic crystal can be written as
[ ( ) ]
ℎani = 𝑔 (𝑐, 𝜁) 𝐾1 𝑚21 𝑚22 + 𝑚22 𝑚23 + 𝑚23 𝑚21 + 𝐾2 𝑚21 𝑚22 𝑚23 (43)
where 𝐾1 and 𝐾2 are the magnetocrystalline anisotropy constants. For the uniaxial crystal, one has
( )
ℎani = 𝑔 (𝑐, 𝜁) 𝐾1 1 − ‖𝐦 ⋅ 𝐞a ‖2 (44)
where 𝐞a is the easy-axis vector. ℎent1 is expressed as
ℎent1 = −𝑔 (𝑐, 𝜁) 𝜇0 𝑀s 𝐇 ⋅ 𝐦 = 𝑔 (𝑐, 𝜁) 𝜇0 𝑀s 𝜙,𝑖 𝑚𝑖 , (45)
and ℎent2 is expressed as

1
ent2 1
ℎ = − 𝜇0 𝐇 ⋅ 𝐇 = − 𝜇0 𝜙,𝑖 𝜙,𝑖 . (46)
2 2
Then the magnetic induction can be derived as
[ ]
𝐵𝑖 = 𝜇0 𝐻𝑖 + 𝑔 (𝑐, 𝜁) 𝑀𝑖 . (47)
The elastic energy density ℎe

will be given and discussed in Section 2.4. According to Eq. (3), the crack surface energy density ℎc
is
𝐺 ( )
ℎc = c 𝑐 2 + 𝑙2 ‖∇𝑐‖2 . (48)
2𝑙
The motivation and formulation of the degradation function 𝑔 (𝑐, 𝜁) deserves to be clarified. The physics behind is that we
introduce 𝑔 (𝑐, 𝜁) to consider the magnetic energy could degrade due to material damage from 𝑐 = 0 (intact) to 𝑐 = 1 (cracked).
An additional parameter 𝜁 ∈ [0, 1] is introduced to control the magnetic degradation behavior of the crack. Following the typical
quadratic degradation function (1 − 𝑐)2 in phase-field fracture model, we propose 𝑔 (𝑐, 𝜁) as
𝑔 (𝑐, 𝜁) = 1 − 2 (1 − 𝜁) 𝑐 + (1 − 𝜁) 𝑐 2 . (49)
2
in which 𝑔 (𝑐, 𝜁 = 0) = (1 − 𝑐) represents that the magnetization dependent magnetic energy is degraded by (1 − 𝑐)2 and is zero in
the crack phase (𝑐 = 1), while 𝑔 (𝑐, 𝜁 = 1) = 1 represents that the magnetic energy is not degraded. For the crack phase (𝑐 = 1),
𝑔 (𝑐 = 1, 𝜁) = 𝜁 reflects the magnetization state of crack phase. In this way, 𝜁 = 1 represents crack with perfect magnetization,
𝜁 = 0 represents crack with no magnetization, and 0 < 𝜁 < 1 represents crack with magnetization multiplied by 𝜁. The curves
7
of degradation functions with different 𝜁 are plotted in Fig. 2(a). It can be found that the degradation function 𝑔 (𝑐, 𝜁) satisfies
certain physical conditions (Kuhn et al., 2015): (1) 𝑔 (𝑐, 𝜁) is monotonically decreasing and continuously differentiable; (2) There is
no degradation for a perfect material 𝑔 (𝑐 = 0, 𝜁) = 1 and a certain magnetization degradation for a completely damaged material
𝑔 (𝑐 = 1, 𝜁) = 𝜁; (3) The driving force of a completely damaged material is zero, i.e., 𝜕𝑔(𝑐,𝜁)
𝜕𝑐
|𝑐=1 = 0. As for the elastic energy density
ℎe , the degradation function is set as 𝑔pm (𝑐) = (1 − 𝑐)2 that is widely used in phase-field fracture model.
A simulation test is conducted to illustrate the influence of parameter 𝜁 on the distributions of magnetic variables (𝐌, 𝐇 and
𝐁) across the crack interface. The simulated square plate is 1 μm × 1 μm with an elliptical crack located in the center, as shown in
Fig. 2(b). The saturation magnetization magnitude 𝑀s is 7.76 × 105 A∕m. An external magnetic field of 106 A∕m is applied along the
𝑦 direction. The distributions of magnetic variables along three vertical lines at 𝑥 = 0.5, 0.7 and 0.9 μm are plotted. Fig. 2(c)–(e),
(f)–(h), (i)–(k) represent the distributions of 𝐌, 𝐇 and 𝐁 along the 𝑦 direction under different
[ 𝜁, respectively.
] 𝐇 field is the sum of
demagnetizing field the external field, i.e., 𝐻𝑖 = −𝜙,𝑖 + 𝐻𝑖ex . 𝐁 is calculated by 𝐵𝑖 = 𝜇0 𝐻𝑖 + 𝑔 (𝑐, 𝜁) 𝑀𝑖 .
As shown in Fig. 2(c)–(k), for the case of 𝜁 = 1, the magnetization in the crack phase is equal to that in the perfect phase.
The obtained distributions of 𝐌, 𝐇 and 𝐁 agree with the analytical solution by Grönefeld and Kronmüller (1989) and Kovacs et al.
(2022). For 𝜁 = 0.5, the magnetization is partially damaged in the crack. The magnetic field 𝐻𝑦 in the crack phase is higher
than that in the perfect phase, owing to the positive and negative magnetic surface charges on the bottom and top surface of
the elliptic crack, respectively. The magnetic induction field 𝐵𝑦 decreases in the transition region from the perfect phase to the
crack phase and then holds nearly constant in the crack phase. Once 𝜁 is 0, the magnetization in the crack phase drops to zero
and the magnetic field 𝐻𝑦 in the crack is even larger. The change of magnetic induction field 𝐵𝑦 has a similar tendency as that
for 𝜁 = 0.5. It can be found that when the magnetization in crack phase is degraded (𝑔(𝑐, 𝜁) < 1), 𝐌 smoothly and continuously
changes across the diffusive crack interface from intact phase (𝑐 = 0) to fully cracked phase (𝑐 = 1), as shown in Fig. 2(c) and (d).
The calculated demagnetization/stray field 𝐇 in the intact and fully cracked phase is largely different (Fig. 2(f) and (g)), agreeing
with the micromagnetic theory (Kronmüller and Fähnle, 2003). The calculated magnetic induction field 𝐁 is continuous across the
interface (Fig. 2(i) and (j)), also agreeing with 𝐁 field continuity in micromagnetism (Kronmüller and Fähnle, 2003). The only
difference is that in our phase-field model here, the interface is diffusive and thus the additional information on 𝐇 and 𝐁 gradually
and smoothly changing in the interface with 0 < 𝑐 < 1 is presented. Especially, 𝐇 becomes continuous across the phase-field diffusive
interface, unlike 𝐇 discontinuity across a sharp material-vacuum interface in micromagnetism (Kronmüller and Fähnle, 2003). If
the interface goes to infinitesimal (i.e., the shadow region in Fig. 2(c), (f) and (i) vanishes), our simulations could restore the 𝐇 and
𝐁 field distributions in traditional micromagnetism.
2.4. Split of the elastic energy density
The total strain consists of the elastic component 𝜀e𝑖𝑗 and the magnetostrictive component 𝜀0𝑖𝑗 , i.e.,
𝜀𝑖𝑗 = 𝜀e𝑖𝑗 + 𝜀0𝑖𝑗 . (50)
𝜀0𝑖𝑗 of a cubic crystal is expressed as
⎧3 ( )
1
⎪ 2 𝜆100 𝑚2𝑖 − 3 , if 𝑖 = 𝑗
𝜀0𝑖𝑗 =⎨ (51)
3
⎪ 2 𝜆111 𝑚𝑖 𝑚𝑗 , if 𝑖 ≠ 𝑗
⎩
where 𝜆100 and 𝜆111 are the magnetostrictive constants on the direction of ⟨100⟩ and ⟨111⟩, respectively. For polycrystalline and
amorphous alloys, the magnetostrictive constants are isotropic, i.e.,
𝜆100 = 𝜆111 = 𝜆s . (52)

Due to the loss of material stiffiness in the smeared crack zone, the elastic strain energy depends not only on the displacement
field 𝑢𝑖 but also on the crack phase field 𝑐. Nevertheless, not the whole elastic energy drives the evolution of the order parameter
𝑐. The spectral split on the elastic strain tensor 𝜺e that was first raised by Miehe et al. (2010a) is adapted here to distinguish the
different role of tension and compression in driving crack prorogation. The elastic energy is expressed as
( ) ( ) ( )
ℎe 𝑐, 𝜺e = 𝑔pm (𝑐) ℎe+ 𝜺e + ℎe− 𝜺e . (53)
The elastic energy density is decomposed into a positive part ℎe+ (𝜺e ) due to tension and a negative part ℎe− (𝜺e ) due to compression,
which are given as
( ) 𝜆⟨ ⟩2 ( )2
ℎe± 𝜺e = tr[𝜺e ] ± + 𝜇tr[ 𝜺e± ] (54)
2
in which 𝜆 and 𝜇 are Lamé constants. The decomposition is based on the positive and negative components of the elastic strain
tensor. Let
𝜺e = 𝐏Λ𝐏T (55)
( )
where 𝐏 consists of the orthonormal eigenvectors of 𝜺e and Λ = diag 𝜆1 , 𝜆2 , 𝜆3 is a diagonal matrix of principal elastic strains. It
can be defined that
𝜺e+ = 𝐏Λ+ 𝐏T (56)
8
Fig. 2. (a) Degradation curves under different parameter 𝜁. (b) Geometry of the center-cracked plate. Distribution of magnetization-vector component 𝑀𝑦 under
(c) 𝜁 = 0, (d) 𝜁 = 0.5 and (e) 𝜁 = 1. Distribution of magnetic-field component 𝐻𝑦 under (f) 𝜁 = 0, (g) 𝜁 = 0.5 and (h) 𝜁 = 1. Distribution of magnetic induction
field 𝐵𝑦 under (i) 𝜁 = 0, (j) 𝜁 = 0.5 and (k) 𝜁 = 1.
𝜺e− = 𝐏Λ− 𝐏T (57)

where
( )
Λ+ = diag ⟨𝜆1 ⟩+ , ⟨𝜆2 ⟩+ , ⟨𝜆3 ⟩+
(58)
Λ− = Λ − Λ+
and
{ {
𝑥 𝑥>0 𝑥 𝑥<0
⟨𝑥⟩+ = ⟨𝑥⟩− = (59)
0 𝑥 ≤ 0, 0 𝑥 ≥ 0.
9
2.5. Irreversible condition for crack
The external force for the crack phase field 𝑓 𝑐 is supposed to be zero. Combining Eq. (36) and the free energy density function,
the evolution of 𝑐 for the uniaxial crystal can be written as
𝜕𝑐 𝐺 { [ ( ) ]}
𝜂 = − c 𝑐 + 2 (1 − 𝑐) ℎe+ + (1 − 𝜁) 𝐴e 𝑚𝑖,𝑗 𝑚𝑖,𝑗 + 𝐾1 1 − ‖𝐦 ⋅ 𝐞a ‖2 + 𝜇0 𝑀s 𝜙,𝑖 𝑚𝑖 + 𝐺c 𝑙∇2 𝑐. (60)
𝜕𝑡 𝑙
Since the load of the specimen can be complicated, the driving force in the energy density ℎ may not increase all the time. The
healing behavior of the crack is not considered in this work. It is assumed that the fracture order parameter 𝑐 will not decrease. Miehe
et al. (2010a) introduced a history variable  (𝑡) to represent the maximum positive strain energy, i.e.,
 (𝑡) = max ℎe+ (𝑡) . (61)

𝑠∈[0,𝑡]
Learned from the work of Miehe et al. (2010a), the history variable is modified. Owing to the additional energy contribution from
micromagnetic structure, the history variable has to be revised as
{ [ ( ) ]}
 (𝑡) = max ℎe+ (𝑡) + (1 − 𝜁) 𝐴e 𝑚𝑖,𝑗 (𝑡) 𝑚𝑖,𝑗 (𝑡) + 𝐾1 1 − ‖𝐦 (𝑡) ⋅ 𝐞a ‖2 + 𝜇0 𝑀s 𝜙,𝑖 (𝑡) 𝑚𝑖 (𝑡) . (62)
𝑠∈[0,𝑡]
Then the governing equation for 𝑐 is derived as

𝜕𝑐 𝐺
𝜂 = − c 𝑐 + 2 (1 − 𝑐)  + 𝐺c 𝑙∇2 𝑐. (63)
𝜕𝑡 𝑙
3. Extension to phase-field fatigue
In the phase field model of fatigue, a degradation function 𝑓 (𝛼)

̄ is introduced to describe the decrease of the fracture toughness.
As Alessi et al. (2018b) and Carrara et al. (2020) considered the influence of the fatigue history, the phase field equation in Eq. (63)
is extended to include fatigue effect, i.e.,
𝜕𝑐 𝑓 (𝛼)
̄ 𝐺c
𝜂 =− ̄ 𝐺c 𝑙∇2 𝑐
𝑐 + 2 (1 − 𝑐)  + 𝑓 (𝛼) (64)
𝜕𝑡 𝑙
where 𝛼̄ is a cumulation describing the accumulated state of fatigue history. It can be derived as
𝑡
𝛼̄ (𝐱, 𝑡) = ̇ |𝛼|d𝜏
𝐻 (𝛼 𝛼) ̇ (65)
∫0
where 𝐻 (𝛼 𝛼) ̇ = 1 if 𝛼 𝛼̇ ≥ 0 (loading) and 𝐻 (𝛼 𝛼)
̇ is the Heaviside function, defined as 𝐻 (𝛼 𝛼) ̇ = 0 otherwise (unloading). The fatigue
variable 𝛼 is supposed to depend on the history variable in Eq. (62), i.e.,
{ [ ( ) ]}
𝛼 = (1 − 𝑐)2 ℎe+ + (1 − 𝜁) 𝐴e 𝑚𝑖,𝑗 𝑚𝑖,𝑗 + 𝐾1 1 − ‖𝐦 ⋅ 𝐞a ‖2 + 𝜇0 𝑀s 𝜙,𝑖 𝑚𝑖 . (66)
Carrara et al. (2020) considered several types of degradation function, from which a concise one is accepted in this work, i.e.,
⎧
⎪1, if 𝛼̄ (𝑡) ≤ 𝛼T
̄ = ⎨( 2𝛼 )2
𝑓 (𝛼) (67)
T
⎪ 𝛼(𝑡)+𝛼 , if 𝛼̄ (𝑡) > 𝛼T
⎩ ̄ T
where 𝛼T is a fatigue threshold under which the fracture toughness of the material experiences no fatigue effects. For precision, the
parameter should be determined on an experimental basis. Borden et al. (2012) analytically studied the homogeneous 1-D solution
of monotonous tensile fracture. By ignoring the spatial derivative of the order parameter and taking 𝜁 = 1, the homogeneous stress
and the homogeneous phase field value can be expressed as
⎧ 𝑙𝐸𝜀2hom ⎧ 𝑙𝐸𝜀2hom
⎪(1 − )2 𝐸𝜀hom , if ℎe+ =  ⎪ 𝐺 +𝑙𝐸𝜀2 , if ℎe+ = 
𝐺c +𝑙𝐸𝜀2hom
𝜎hom =⎨ 𝑐hom =⎨ c hom (68)
⎪(1 − 2𝑙
)2 𝐸𝜀hom , if ℎe+ < , ⎪ 2𝑙 , if ℎe+ < .
⎩ 2𝑙+𝐺c ⎩ 2𝑙+𝐺c
The critical stress and the corresponding strain are found as

√ √
9 𝐸𝐺c 𝐺c
𝜎c = , 𝜀c = . (69)
16 3𝑙 3𝑙𝐸
Carrara et al. (2020) and Simoes and Martínez-Pañeda (2021) tended to assume the threshold as
1
𝛼T = 𝜀 𝐸𝜀 (70)
2 c c
which is equal to the elastic strain energy at the corresponding critical strain under the linear elastic assumption. In this work, we
propose that the threshold parameter should be the critical energy calculated by
𝜀c 𝐺c
𝛼T = 𝜎d𝜀 = . (71)
∫0 8𝑙
10
4. Finite-element implementation
The phase-field model is solved by finite element method (FEM). The governing equations are discretized spatially and
temporally. In the 3D simulation, the eight-node linear elements are applied to mesh the body. Method of weighted residuals is
used to solve the equations. Galerkin method sets possible test functions as the weighted functions. Then the solution( approximate
)
to the realistic solution can be calculated under the equilibrium of the whole
( integral )domain. The displacement 𝐮 𝑢1 , 𝑢2 , 𝑢3 , the
scalar magnetic potential 𝜙, the magnetization-vector order parameters 𝐌 𝑀1 , 𝑀2 , 𝑀3 and the crack phase-field order parameter
𝑐 are taken as the degrees of freedom.
Lagrangian multiplier method is applied to satisfy the constraint ‖𝐌‖ = 𝑀s at an constant temperature. In this way, an additional
functional is given as
( ) 𝜆2𝑀
𝛱𝑐 = 𝜆𝑀 ‖𝐌‖2 − 𝑀s d𝑉 − d𝑉 (72)
∫𝛺 ∫𝛺 2𝑘a
where 𝜆𝑀 is the Lagrangian multiplier and 𝑘a is a constraint parameter. The associated functional variation is
( )
𝜆
𝛿𝛱𝑐 = ‖𝐌‖2 − 𝑀s − 𝑀 𝛿𝜆d𝑉 + 2𝜆 𝐌 ⋅ 𝛿𝐌d𝑉 . (73)
∫𝛺 𝑘a ∫𝛺 𝑀
The governing equations are converted into the weak form on the integral domain, i.e.,
𝑢
0= 𝜎𝑖𝑗 𝛷𝑖,𝑗 d𝑉 − 𝜎𝑖𝑗 𝑛𝑗 𝛷𝑖𝑢 d𝑆, (74)
∫𝛺 ∫𝜕𝛺
0= 𝐵𝑖 𝛷,𝑖𝜙 d𝑉 − 𝐵𝑖 𝑛𝑖 𝛷𝜙 d𝑆, (75)

∫𝛺 ∫𝜕𝛺
[( ) ]
𝜇0 𝛼 𝛼 𝜕ℎ 𝜕ℎ 𝜕ℎ
0= −𝑀̇ − 𝜖 𝑀̇ 𝑀 − 𝛷𝑖𝑀 − 𝛷𝑀 d𝑉 + 𝛷𝑀 𝑛 d𝑆, (76)
∫𝛺 𝛾0 𝑀s 𝑖 𝛾0 𝑀s2 𝑖𝑗𝑘 𝑗 𝑘 𝜕𝑀𝑖 𝜕𝑀𝑖,𝑗 𝑖,𝑗 ∫𝜕𝛺 𝜕𝑀𝑖,𝑗 𝑖 𝑗
[ ]
𝜕𝑐 𝑓 (𝛼)̄ 𝐺c
0= 𝜂 + 𝑐 − 2 (1 − 𝑐)  𝛷𝑐 + 𝑓 (𝛼)
̄ 𝐺c 𝑙𝑐,𝑖 𝛷,𝑖𝑐 d𝑉 − ̄ 𝐺c 𝑙𝑐,𝑖 𝑛𝑖 𝛷𝑐 d𝑆
𝑓 (𝛼) (77)
∫𝛺 𝜕𝑡 𝑙 ∫𝜕𝛺
where 𝛷𝑖𝑢 , 𝛷𝜙 , 𝛷𝑖𝑀 , 𝛷𝑐 are the test functions for 𝑢𝑖 , 𝜙, 𝑀𝑖 and 𝑐, respectively. By introducing the shape functions 𝑁 𝐼 for the
independent variables and test functions, they can be discretized as
∑ ∑ ∑ ∑
𝑢𝑖 = 𝑁 𝐼 𝑢𝐼𝑖 , 𝜙 = 𝑁 𝐼 𝜙𝐼 , 𝑀𝑖 = 𝑁 𝐼 𝑀𝑖𝐼 , 𝑐 = 𝑁 𝐼 𝑐𝐼
∑ ∑ ∑ ∑ (78)
𝛷𝑖𝑢 = 𝑁 𝐼 𝛷𝑖𝑢𝐼 , 𝛷𝜙 = 𝑁 𝐼 𝛷𝜙𝐼 , 𝛷𝑖𝑀 = 𝑁 𝐼 𝛷𝑖𝑀𝐼 , 𝛷𝑐 = 𝑁 𝐼 𝛷𝑐𝐼 .
The superscript 𝐼 is the node number. The additional residual from Lagrangian multiplier method in Eq. (73) needs to be added
into Eq. (76). Meanwhile, the residual has to be set for the additional degree of freedom 𝜆𝑀 according to Eq. (73). Inserting the
discretization into the weak forms, one can obtain the elemental residuals
𝑅𝐼𝑢 = 𝜎𝑖𝑗 𝑁,𝑗𝐼 d𝑉 − 𝜎𝑖𝑗 𝑛𝑗 𝑁 𝐼 d𝑆, (79)

𝑖 ∫𝛺 ∫𝜕𝛺
( )
𝑅𝐼𝜙 = 𝐵𝑖 𝑁,𝑖𝐼 d𝑉 = 𝜇0 −𝜙,𝑖 + 𝑀s 𝑚𝑖 𝑁,𝑖𝐼 d𝑉 − 𝐵𝑖 𝑛𝑖 𝑁 𝐼 d𝑆, (80)
∫𝛺 ∫𝛺 ∫𝜕𝛺
[( ) ]
𝜇0 𝛼 𝛼 𝜕ℎ 𝜕ℎ 𝜕ℎ
𝑅𝐼𝑀 = − 𝑀̇ 𝑖 − 𝜖𝑖𝑗𝑘 𝑀̇ 𝑗 𝑀𝑘 − 𝑁𝐼 − 𝑁,𝑗𝐼 + 2𝜆𝑀 𝑀𝑖 𝑁 𝐼 d𝑉 + 𝑁 𝐼 𝑛𝑗 d𝑆, (81)
𝑖 ∫𝛺 𝛾0 𝑀s 𝛾0 𝑀s2 𝜕𝑀𝑖 𝜕𝑀𝑖,𝑗 ∫𝜕𝛺 𝜕𝑀𝑖,𝑗
[ ]
𝜕𝑐 𝑓 (𝛼) ̄ 𝐺c
𝑅𝐼𝑐 = 𝜂 + 𝑐 − 2 (1 − 𝑐)  𝑁 𝐼 + 𝑓 (𝛼) ̄ 𝐺c 𝑙𝑐,𝑖 𝑁,𝑖𝐼 d𝑉 − ̄ 𝐺c 𝑙𝑐,𝑖 𝑛𝑖 𝑁 𝐼 d𝑆,
𝑓 (𝛼) (82)
∫𝛺 𝜕𝑡 𝑙 ∫𝜕𝛺
( )
𝜆
𝑅𝐼𝜆 = 𝑀12 + 𝑀22 + 𝑀32 − 𝑀s − 𝑀 𝑁 𝐼 d𝑉 . (83)
𝑀 ∫𝛺 𝑘a
Based on Eqs. (81) and (82), the Neumann boundary conditions for the order parameters 𝑀𝑖 and 𝑐 can be given as
𝜕ℎ 2𝑔(𝑐, 𝜁)𝐴e
𝑛 = ̄𝑖
𝑀𝑖,𝑗 𝑛𝑗 = 𝑀 on 𝜕𝛺MN , (84)
𝜕𝑀𝑖,𝑗 𝑗 𝑀s2
̄ 𝐺c 𝑙𝑐,𝑖 𝑛𝑖 = 𝑐̄
𝑓 (𝛼) on 𝜕𝛺𝑐N (85)
where 𝑀 ̄ 𝑖 and 𝑐̄ are the specified values on the boundary part 𝜕𝛺MN and 𝜕𝛺𝑐N , respectively. In the micromagnetic theory, no
exchange torque acts at the free surface (Szambolics et al., 2008) and the gradient of 𝑀𝑖 along the outward normal direction is
zero (Miehe and Ethiraj, 2012), resulting in
𝑀𝑖,𝑗 𝑛𝑗 = 0 on 𝜕𝛺MN . (86)
Similarly, for the crack phase field, one has
𝑐,𝑖 𝑛𝑖 = 0 on 𝜕𝛺𝑐N . (87)
11
For the Dirichlet boundary conditions for 𝑀𝑖 and 𝑐, one has

̂𝑖
𝑀𝑖 = 𝑀 on 𝜕𝛺MD , (88)
and
𝑐 = 𝑐̂ on 𝜕𝛺𝑐D . (89)
The array of all the nodal variables is
𝐝𝐼 = [𝑢𝐼1 , 𝑢𝐼2 , 𝑢𝐼3 , 𝜙𝐼 , 𝑀1𝐼 , 𝑀2𝐼 , 𝑀3𝐼 , 𝑐 𝐼 , 𝜆𝐼𝑀 ]T . (90)
The implicit backward Euler time integration is adopted in the temporal discretization, i.e.,
𝐝𝑛+1 − 𝐝𝑛
𝐝̇ = . (91)
𝛥𝑡
𝐑𝐼𝑛 and 𝐝𝐽𝑛 are the previous variables at time step 𝑡𝑛 . The present quantity at time step 𝑡𝑛+1 is the function of the past ones, i.e.,
( )
𝐝𝐽𝑛+1 − 𝐝𝐽𝑛
𝐑𝐼𝑛+1 = 𝐑𝐼 𝐝𝐽𝑛+1 , . (92)
𝛥𝑡
The nonlinear equations are computed with the method of Newton iteration. The iteration matrix comes into the form
1 𝐼𝐽
𝐒𝐼𝐽 = 𝐊𝐼𝐽 + 𝐃 (93)
𝛥𝑡
where 𝐊𝐼𝐽 is the stiffness matrix and 𝐃𝐼𝐽 is the damping matrix. The stiffness matrix 𝐊𝐼𝐽 is the negative derivative of the residual
with respect to the degree of freedom while the damping matrix 𝐃𝐼𝐽 is the negative derivation of the residual with respect to the
rate. After completing the matrixes, the programming is performed under the open source Multiphysics Object Oriented Simulation
Environment (MOOSE).
5. Crack driving force
The proper crack driving forces for micromagnetic-mechanical fracture is critical. This section will discuss whether the magnetic
energy drives the crack propagation in magnetostrictive alloys. A one-dimensional case based on the Stoner–Wohlfarth model is
investigated. √
In the initial state, the magnetostriction strain of the one-dimensional bar is assumed to be zero with the angle
𝜃 = arccos(1∕ 3) between the 𝑥 direction and the magnetization vector, as depicted in Fig. 3(a). The strain in the bar is expressed
as 𝜀 = 𝜀0 = 𝜀e = 0 and the bar is at the stress-free state. The tensile loading state is shown in Fig. 3(b). It is supposed that
the one-dimensional bar is stretched in the 𝑥 direction which is also parallel to the direction of magnetic field. In the case of
Stoner–Wohlfarth model, the total magnetic enthalpy is simplified as
[ ( )] ( )
1 3 1 2
ℎ = (1 − 𝑐)2 𝐸 𝜀 − 𝜆s cos2 𝜃 − + 𝑔 (𝑐, 𝜁) −𝜇0 𝑀s 𝐻 cos 𝜃 + 𝐾1 sin2 𝜃 . (94)
2 2 3
The driving force 𝑓 is formulated as
𝜕ℎ
𝑓 =− . (95)
𝜕𝑐
5.1. Evolution of fracture in case of pure mechanical driving force
For the case of pure mechanical driving force when the variable 𝜁 = 1 and thus the magnetic energy has nothing to do with the
fracture, the driving force 𝑓 is expressed as
[ ( )]
𝜕ℎela 3 1 2
𝑓 =− = 𝐸 (1 − 𝑐) 𝜀 − 𝜆s cos2 𝜃 − . (96)
𝜕𝑐 2 3
As proposed by Miehe et al. (2010b), a threshold function 𝜏 for rate-independent and discontinuous evolution of the order parameter
𝑐 is postulated, i.e.,
𝐺c
𝜏=𝑓− 𝑐. (97)
𝑙
𝐺
It characterizes that there is no diffusive crack accumulation for 𝜏 < 0. The factor 𝑙c models the resistance of the material to
damage. The necessary conditions for the evolution of the crack phase field order parameter can be proposed as
𝑐̇ ≥ 0, 𝜏 ≤ 0, ̇ = 0.
𝑐𝜏 (98)
In case of the low driving force, the energy that is needed to overcome the threshold cannot be satisfied and the fracture order
parameter will not evolve. Based on the critical state 𝜏 = 0, the solution of 𝑐 can be obtained analytically. In the presence of
12
Fig. 3. Schematic representation of the one-dimensional tensile model: (a) the initial state with zero strain and stress, (b) the loading state under tension and
magnetic field along 𝑥 direction. Mechanical responses under different driving forces: (c) stress–strain curves and (d) order parameter 𝑐 evolution under uniaxial
tension and different magnetic fields in consideration of the pure mechanical part of the driving force, (e) stress–strain curves and (f) order parameter 𝑐 evolution
under uniaxial tension and different magnetic fields in consideration of the full micromagnetic-mechanical driving force.
positive magnetostrictive strain, the initial self-equilibrium state without external magnetic field is corresponding to a negative
elastic strain, i.e.,
𝜀e ≤ 0, 𝑐 = 0 at 𝑡 = 0. (99)
The critical state 𝜏 = 0 gives

[ ( )]2
𝜀e > 0, 𝑐=
N , N = 𝐸 𝜀 − 3 𝜆s cos2 𝜃 − 1 at 𝑡 = 0. (100)
𝐺c 2 3
N+ 𝑙
The fracture irreversibility must also be considered. Defining  (𝑡) as the maximum value of N obtained in the time history, one
has
 (𝑡)
𝑐 (𝑡) = 𝐺c
∈ [0, 1] (101)
 (𝑡) + 𝑙
where
[ ]
 (𝑡) ∶= max N (𝑠) , 0 ∈ [0, ∞). (102)
𝑠∈[0,𝑡]
13
Table 1
Comparison of the critical stresses for the pure mechanical driving force.
B = 0 B = 0.05 T B = 0.1 T
Analytic 72.62 MPa 72.62 MPa 72.62 MPa
FEM 72.62 MPa 72.61 MPa 72.61 MPa
Table 2
Comparison of the critical stresses for the micromagnetic-mechanical driving force.
B = 0 B = 0.05 T B = 0.1 T
Analytic 72.62 MPa 79.30 MPa 89.52 MPa
FEM 72.62 MPa 79.51 MPa 89.91 MPa
With the result of 𝑐 at hand, the stress state of the material can be calculated as
[ ( )]
3 1
𝜎 = (1 − 𝑐)2 𝐸 𝜀 − 𝜆s cos2 𝜃 − . (103)
2 3
5.2. Evolution of fracture in case of micromagnetic-mechanical driving force
In the consideration of micromagnetic-mechanical driving force, the impermeable crack behavior is followed with the variable
𝜁 = 0. The micromagnetic-mechanical driving force is written as
{ [ ( )] }
𝜕ℎ 3 1 2
𝑓 =− = (1 − 𝑐) 𝐸 𝜀 − 𝜆s cos2 𝜃 − − 2𝜇0 𝑀s 𝐻 cos 𝜃 + 2𝐾1 sin2 𝜃 . (104)
𝜕𝑐 2 3
The energy density function N is derived as
{
−2𝜇0 𝑀s 𝐻 cos 𝜃 + 2𝐾1 sin2 𝜃, if 𝜀e ≤ 0
N= (105)
𝐸[𝜀 − 32 𝜆s (cos2 𝜃 − 13 )]2 − 2𝜇0 𝑀s 𝐻 cos 𝜃 + 2𝐾1 sin2 𝜃, if 𝜀e > 0
For short, the evolution of stress and the order parameter are similar to that in the previous Section 5.1.
5.3. Comparison of the results in the two cases of driving force
In the one-dimensional simulation, the material parameters are chosen to 𝐸 = 50 GPa, 𝜆s = 0.001, 𝜇0 = 1.257 × 10−6 H/m,
𝐺
𝑀s = 7.76 × 105 A/m, 𝐾1 = 2 × 103 J/m3 and 𝑙c = 106 J/m3 . In order to reveal effect of the external magnetic field, an external
magnetic field of 0, 0.05 and 0.1 T is applied. Both the results of the two cases are obtained from analytical calculation and FEM
simulation. The strain–stress curves and the relations between strain and order parameter 𝑐 are depicted in Fig. 3(c)–(f).
It can be seen that the FEM simulation results well match the analytic results. In both cases, the existence of magnetic field makes
the bar at a compressive state when the load is low. The external magnetic field delays the increase of the order parameter 𝑐. The
curves also depict that the micromagnetic-mechanical driving force leads to the higher maximum stress under an external magnetic
field. In contrast, the maximum stresses for different magnetic fields hold the same in the case of pure mechanical driving force. The
critical stresses are listed in Tables 1 and 2. The experiment result is far few to prove which driving force is realistic. Clatterbuck
et al. (2000) and Yamaguchi et al. (2003) discovered that the fracture properties of Incoloy 908 is not significantly affected by the
magnetic fields. In contrast, Shindo et al. (2004, 2005, 2006, 2008), Narita et al. (2016) and Peron et al. (2019) conducted various
experiments and concluded that the fracture load of Terfenol-D is greater in absence of the magnetic field. The conclusions from
both aspects fail to be completely convincing. The influence of magnetic field on the fracture properties is still inconclusive and
remains to be further investigated. Similar with the conclusion from Sridhar and Keip (2019), here the pure mechanical driving
force is accepted and 𝜁 = 1 is utilized in the following investigations.
6. Simulation results and discussions
The simulations are performed for an amorphous Terfenol-D film that is a typical giant magnetostrictive material. The material
parameters are collected from the literatures and listed in Table 3. Simulations on the tension and shear of single-edge notched plate,
three-point bending, tension of a plate with an elliptical inclusion, and cyclic loading under external magnetic field are carried out
to demonstrate the ability of the proposed model.
6.1. Monotonic tension under horizontal magnetic field
We consider a squared plate with a horizontal notch which is placed at middle height from the left outer surface to the center
of the specimen. The geometry is depicted in Fig. 4(a). The plate size is set as 1000 nm × 1000 nm × 10 nm and the notch is 500 nm.
14
Table 3
Material parameters of Terfenol-D.
Parameter Name Value Unit
𝐸 Young’s modulus 50 GPa
𝜈 Poisson’s ratio 0.3
𝐾1 Anisotropy coefficient 2 × 103 J/m3
𝐾2 Anisotropy coefficient 0 J/m3
𝑀s Saturation magnetization 7.76 × 105 A/m
𝜆100 Magnetostrictive constant 1 × 10−3
𝜆111 Magnetostrictive constant 1 × 10−3
𝜇0 Vacuum permeability 1.257 × 10−6 H/m
𝐴e Exchange coefficient 10−11 J/m
𝑘a Constraint parameter 104
𝐺c Critical fracture toughness 2 × 10−2 N/m
𝑙 Crack surface width 20 nm
Fig. 4. Geometry and boundary conditions of (a) single-edge notched tension, (b) single-edge notched pure shear, (c) symmetric three-point bending, (d)
single-edge notched tension of a plate with an elliptical inclusion.
The initial magnetization is chosen to be 𝑚1 = 1 for the whole plate. The magnetic easy-axis is supposed to be (1, 0, 0). The tensile
load is applied on the top of the specimen with a constant displacement increment of 𝛥𝑢 = 1 × 10−3 nm. An external magnetic field
of 0, 0.1 and 1 T is applied along the 𝑥 direction by giving the boundary conditions of the magnetic potential 𝜙. Both the positive
and negative magnetic fields are investigated.
The load–displacement curves are shown in Fig. 5(a). The global displacement of the plate is set to zero initially, so the plate is
experiencing the tensile force in the 𝑦 direction due to the magnetostrictive strain. Before the crack propagation, the stress at the
crack tip increases with the load. Upon the critical circumstance has reached, the load decreases rapidly. The positive magnetic fields
parallel to the 𝑥 direction will decrease the maximum forces and the critical displacement, while the negative ones will increase the
loads and deformation.
Fig. 5(b) shows the absolute value of the vertical component of magnetization unit vector 𝑚2 near the crack tip as a function
of 𝑢2 . It is found that upon loading, the magnetization vector tends to switch and point to the crack tip. Due to the principle of
energy minimization, 𝑚2 changes with the increasing displacement load and decreases the elastic strain. The external magnetic
fields influence the switching behavior of magnetization vectors. The curves in Fig. 5 show apparent asymmetric load–displacement
response under external fields of ±1 T. This asymmetry is attributed to the local magnetization vector rotation and the local
magnetostrictive strain change around the crack tip. A magnetic field of 1 T keeps the initially uniform magnetization vectors
distributed along +𝑥 direction, and the local magnetization vectors around the crack tip seldomly rotate. In contrast, a magnetic
15
Fig. 5. (a) Load–displacement curves and (b) changes of 𝑚2 around the crack tip for single-edge notched tension under external magnetic fields.
Fig. 6. (a) The crack propagation path under tension with −1 T magnetic field along 𝑥 direction. Local magnetization-vector (indicated by arrows) distribution
around the crack tip under a displacement loading of (b) 𝑢2 = 0.3 nm, (c) 𝑢2 = 0.5 nm, (d) 𝑢2 = 0.7 nm and (e) 𝑢2 = 0.9 nm.
field of −1 T makes the local magnetization vectors around the crack tip drastically rotate away from the initial +𝑥 direction. This
drastic rotation can be seen from Fig. 5(b) that presents large increase of vertical component of magnetization vector (𝑚2 ). Moreover,
the arrows in Fig. 6 also show the drastic rotation of magnetization vectors around the crack tip. Fig. 6 presents the distributions
of magnetization vectors around the crack tip under an external magnetic field of −1 T. The local contour plot of crack pattern
and magnetization-vector distributions under displacements of 𝑢2 = 0.3, 0.5, 0.7, and 0.9 nm are given in Fig. 6(b), (c), (d), and
(e), respectively. Although the macroscopic magnetostrictive strain (without the consideration of local magnetization distribution
difference) under ±1 T is the same, the local magnetization vector distribution and thus the local magnetostrictive strain around
the crack tip are different, finally leading to asymmetric load–displacement response.
Under a low magnetic field, the increasing 𝑚2 will decrease the magnetostriction 𝜀022 if 𝑚22 does not exceed 13 . The single-edge
notched tension specimen mainly bears the normal stress in the 𝑦 direction. The positive elastic energy density is dominated by the
term 12 𝜎22 𝜀e22 . Qualitatively, it can be approximated as
1 1( e )
ℎe+ =̇ 𝜎22 𝜀e22 = 𝜆𝜀𝑘𝑘 + 2𝜇𝜀e22 𝜀e22
2 2
[ ( ) ] (106)
1 9 2 3 1 1
= 𝜇𝜆 𝑚4 + − 𝜆100 𝜆𝜀e𝑘𝑘 − 3𝜆2100 𝜇 − 6𝜇𝜆100 𝜀22 𝑚22 + 𝜆𝜀e𝑘𝑘 𝜀22 + 𝜆100 𝜆𝜀e𝑘𝑘 + 2𝜇𝜀222 + 𝜇𝜆2100 + 2𝜇𝜆100 𝜀22 .
2 2 100 2 2 2 2
1 𝜆𝜃+4𝜀22 𝜇
The elastic energy can be the function of 𝑚2 . When 𝑚22 tends to reach 3
+ 6𝜆100 𝜇
, the elastic energy could achieve a minimum.
The positive magnetic field will inhibit the increase of 𝑚22 while the negative one will raise the amplitude of 𝑚22 . Under 1, 0.1 and
−0.1 T, the values of 𝑚2 have little variation. So it can be observed in Fig. 5(a) that the critical loads are close. Compared to −0.1 T,
a magnetic field of −1 T is quite large and causes a significant variation of 𝑚2 . The change of 𝑚2 decreases the elastic energy and
thus makes the critical load much higher.
16
Fig. 7. Load–displacement curves for single-edge notched tension (𝐸 = 1 GPa) with vertical magnetization (magnetized case) and without magnetization
(non-magnetized case).
6.2. Monotonic tension under vertical magnetic field
Inspired by Moreno-Mateos et al. (2023), we also examine the influence of magnetization on the loading capacity. In detail,
tension simulations of the single-edge notched sample with uniform magnetization vector along the vertical direction (magnetized
case) and without magnetization (non-magnetized case) are carried out, as shown in Fig. 7. Considering that the magnetorheological
elastomers studied by Moreno-Mateos et al. (2023) is soft with extremely low Young’s modulus, here we also reduce the Young’s
modulus to 1 GPa in order to exhibit a more notable magnetization dependent phenomenon. The rest material parameters in
this subsection are the same as in Table 3. The calculated load–displacement curves in Fig. 7 indicate that both the critical load
and displacement of the magnetized case are much higher than that of the non-magnetized case. This result is consistent with
that reported by Moreno-Mateos et al. (2023). It should be noted that experiments and simulations conducted by Moreno-Mateos
et al. (2023) show that the magnetorheological elastomers containing magnetized particles had stronger fracture toughness owing
to the remanent magnetization from the magnetized particles. Since our model here is micromagnetic-mechanically coupled at
the magnetic-domain scale, magnetization vector is taken as the degree of freedom and thus we could solve the magnetization
itself induced demagnetization/stray field in both the material and vacuum. Through the demagnetization/stray field (the origin of
magnetic field generated by magnetized particles in the work by Moreno-Mateos et al. (2023)), opposed magnetized sides of a crack
will naturally interact to influence the crack behavior. In addition, our model could resolve the magnetization vector or magnetic
domain switching around the crack tip, providing more microscale information about the magnetization dependent tension behavior.
In brief, compared to the work of Moreno-Mateos et al. (2023), the novelty and improvement of our model could lie in three
aspects: (1) mesoscale model including micromagnetics and allowing for the interaction between micromagnetic domain and crack
phase field; (2) micromagnetic-mechanically coupled phase-field fracture model with an additional LLG equation for governing the
spatial/temporal evolution of magnetization vector; (3) further extension to micromagnetic-mechanically coupled fatigue fracture.
6.3. Pure shear under magnetic field
The pure shear test is investigated by the single-edge notched plate in Fig. 4(b). There is a horizontal displacement load on
the top of the specimen and the bottom is fixed. The constant monotonic displacement increment is set as 𝛥𝑢 = 1 × 10−2 nm. It is
anticipated that the crack would propagate with a angle of 45◦ to the 𝑥 direction. The simulation results with an external magnetic
field of 1 and −1 T are discussed.
The simulated load–displacement curves are shown in Fig. 8(a) and (b). The external magnetic fields have tiny effects on the shear
test. The critical load and the magnetization are the same under different magnetic fields. The crack pattern and the magnetization-
vector distribution under −1 T are presented in Fig. 8(c)–(g). The crack propagation path is similar to that in the non-magnetostrictive
materials. It can be seen that there is no obvious change in the magnetization vector during the loading. The negative magnetic
field can only make the magnetization rotate slightly. Different from the tension test, magnetization vectors in the domain around
the crack tip do not point to the crack path.
In the shear test, the crack propagation is mainly driven by the shear stress which is related to the magnetostriction 𝜀012 . The
positive elastic energy density can be approximated as
( )2
1( ) ( )2 3
ℎe+ =̇ 𝜎 𝜀e + 𝜎21 𝜀e21 = 𝜎12 𝜀e12 = 2𝜇 𝜀e12 = 2𝜇 𝜀12 − 𝜆111 𝑚1 𝑚2 . (107)
2 12 12 2
The elastic energy is dominated by 𝑚1 𝑚2 . The competition among ℎani , ℎent and ℎe makes a little increase of 𝑚2 . The results of
𝑚1 𝑚2 are all close to 0 under different conditions, and thus magnetic field has slight effect on the elastic energy in the shear test.
Therefore, the critical force are almost the same for different cases.
17
Fig. 8. Load–displacement curves for single-edge notched pure shear under a (a) positive and (b) negative magnetic field. (c) The crack path of the specimen
under the shear load. The contour plots of 𝑚2 and shear crack tip pattern under the −1 T at different displacements according to (d)–(g) marked in (b).
6.4. Symmetric three-point bending under magnetic field
Another fracture benchmark is three-point bending test. A compressive load is imposed on a simply supported notched beam.
The geometric setup along with the loading conditions are illustrated in Fig. 4(c). The beam is initially constructed as 800 nm ×
200 nm × 10 nm with a single micromagnetic domain state 𝑚1 = 1. The notch is 20 nm wide and 40 nm high. The computation is
performed in a monotonic displacement driven context with a constant displacement increment of 𝛥𝑢 = 1 × 10−2 nm. A magnetic
field of 0.1 and 1 T is applied along the 𝑥 direction.
The crack path is shown in Fig. 9(a). The crack propagates along the beam center. The global mechanical responses under
different magnetic fields are presented in Fig. 9(b). The critical load descends due to the parallel magnetic field. The magnetization-
vector distribution under no magnetic field is depicted in Fig. 9(c)–(f). In the three-point bending test, the magnetization vectors
encircle the notch clockwisely. In the tests, the strain around the crack consists of the compressive 𝜀e22 and the tensile 𝜀e11 . So the
fracture is mainly driven by the stress 𝜎11 . It is similar to the single edge tension test that the initial 𝑚1 = 1 will decrease to make the
elastic strain lower. Thus, the total energy could be at a minimum state. The external magnetic field prevents the magnetization from
rotating and leads to a smaller critical load. It should be noted that in the work of Colussi et al. (2016), both of the experimental
and FEM results show the decrease of critical loads under external magnetic fields. Therefore, our phase-field fracture simulation
results here quantitatively agree with the experimental reports (Colussi et al., 2016).
6.5. Tension of a plate with an elliptical inclusion
Besides the single-domain situation, a complex structure with an elliptical inclusion is simulated, as shown in Fig. 4(d).
In a squared magnetostrictive plate with an elliptical inclusion, a short notch is set in the left edge. The plate size is set as
1000 nm × 1000 nm × 10 nm and the notch is 200 nm. The ellipse positioned in the plate center measures 300 nm in width and
240 nm in height. Inside the inclusion, the initial magnetization is 𝑚2 = 1 and the magnetic easy-axis is (0, 1, 0). The tensile boundary
condition applied on the top has a constant displacement increment of 𝛥𝑢 = 1 × 10−3 nm in each time step. The magnetic boundary
condition is set by giving specific values to the magnetic potential 𝜙. A magnetic field of 0, −1 and 1 T along the 𝑥 direction is
considered.
18
Fig. 9. (a) The crack path of three-point bending test. (b) Load–displacement curves for three-point bending test under positive magnetic fields.
Magnetization-vector distribution and crack tip pattern without magnetic field at different displacements according to (c)–(f) marked in (b).
The load–displacement curves are shown in Fig. 10(a). The initial magnetization distribution is given in Fig. 10(b). The crack path
and the final magnetic domain structure under different magnetic fields are presented in Fig. 10(c)–(e). It is found that the external
magnetic fields decrease the critical load regardless of the direction of the magnetic field. However, −1 T makes the specimen endure
largest deformation. The crack propagates horizontally and crosses through the inclusion. The plate in Fig. 10(c) tends to form a
multi-domain structure from the initial single-domain state. The magnetization vectors under 1 T present a single-domain state and
are parallel to the magnetic field, as shown in Fig. 10(d). A magnetic field of −1 T makes the magnetization vector rotate towards
the negative 𝑥 direction, resulting in magnetostrictive strain to compensate the elastic strain and thus the endurement of larger
external displacement, as shown in Fig. 10(e).
In another example, the fracture toughness of the elliptical magnetic particle is set as 0.2 N/m, which is 10 times as that in the
surrounding. The magnetic boundary condition is set to be open circuited, i.e., 𝐵𝑖 𝑛𝑖 = 0. Fig. 11(a) shows the load–displacement
relation and Fig. 11(b)–(e) present the crack pattern and the magnetic domain structure. Due to the open-circuited boundary
condition, the specimen exhibits a multi-domain state. There are two drops in the load–displacement curve in Fig. 11(a). At first, the
crack propagates horizontally from the notch and the load decreases slightly. When the crack reaches the left edge of the elliptical
inclusion, the larger fracture toughness of the inclusion prevents the crack from crossing through the elliptical region. Then the
external load continues to rise while the crack tends to propagate around the inclusion. In the tri-junction of the magnetic domains,
the elastic strain 𝜀e22 is remarkably large, as shown in Fig. 11(f). Therefore, a new crack nucleates in that tri-junction point and then
propagates horizontally towards both sides, as shown in Fig. 11(d) and (e).
6.6. Fatigue cracking behavior under magnetic field
A cyclic tensile test of the single-edge notched plate is performed to examine the fatigue cracking behavior. The geometry and
the micromagnetic domain structure are the same as that in Section 6.1. A cyclic tensile load is applied at the top of the plate with
𝑢 = 3 × 10−2 × (− cos (𝑡) + 1) nm. A magnetic field of 1 and 10 T is parallel to the 𝑥 axis.
The crack length vs cycle curves under positive magnetic fields are presented in Fig. 12(a). The crack path and the magnetization-
vector distribution are given in Fig. 12(b)–(g). The positive magnetic fields slightly decrease the fatigue life. Similar to the monotonic
single-edge notched tension test, the magnetization vector tends to point to the crack. 𝑚2 around the crack tip slightly increases
with the displacement loading. The positive magnetic fields would force magnetization vectors parallel to the positive 𝑥 direction,
i.e., 𝑚1 is close to 1 and 𝑚2 is almost zero.
The case of negative magnetic fields is different, as shown in Fig. 13. Unlike the positive magnetic fields, the negative ones could
improve the fatigue resistance and fatigue life. It is found that a external field of −1 T increases the fatigue life by 50%, while −10 T
field only slightly enhance the fatigue life. This different role of −1 and −10 T field could be attributed to the distinct magnetization
19
Fig. 10. (a) Load–displacement curves for a plate with an elliptical inclusion that possesses a magnetic easy-axis along the vertical direction. (b) Magnetization-
vector distribution at the initial state. Magnetization-vector distribution and crack path after the crack goes throughout the specimen under the magnetic field
of (c) 0 T, (d) 1 T, (e) −1 T.
Fig. 11. (a) Load–displacement curve for a plate with an elliptical inclusion that possesses a higher 𝐺c (0.2 N/m). Magnetization-vector distribution and crack
pattern at different displacements according to (b)–(e) marked in (a). (f) Distribution of elastic strain 𝜀e22 before the second crack propagates.
switching behaviors around the crack. At a high magnetic field of −10 T, the magnetic driving force is strong enough to completely
switch the micromagnetic single domain with magnetization along the positive 𝑥 direction (𝑚1 = 1) to along the negative 𝑥 direction
(𝑚1 = −1) in a quite short time (several cycles), as shown in Fig. 13(b)–(d). In contrast, at a low magnetic field of −1 T, the magnetic
20
Fig. 12. (a) Fatigue cracking behavior of single-edge notched tension under positive magnetic fields. Magnetization-vector distribution and crack tip pattern
without magnetic field at (b) 95 cycles, (c) 190 cycles, (d) 238 cycles. Magnetization-vector distribution and crack tip pattern under 1 T at (e) 95 cycles, (f)
190 cycles, (g) 238 cycles.
Fig. 13. (a) Fatigue cracking behavior of single-edge notched tension under negative magnetic fields. Magnetization-vector distribution and crack tip pattern
under −10 T at (b) 14 cycles, (c) 190 cycles, (d) 238 cycles. Magnetization-vector distribution and crack tip pattern under −1 T at (e) 175 cycles, (f) 190 cycles,
(g) 238 cycles, (h) 334 cycles. The distribution of fracture toughness degradation function along the crack path under −1 T at (i) 175 cycles, (j) 190 cycles, (k)
238 cycles, (l) 334 cycles.
driving force cannot induce 180◦ magnetization switching suddenly and thus the magnetization vectors gradually rotate during the
whole fatigue cracking process, as shown in Fig. 13(e)–(h). The distribution of fracture toughness degradation function 𝑓 (𝛼) ̄ under
−1 T is shown in Fig. 13(i)–(l). There are four obvious stages under −1 T field (Fig. 13(a)). In the first stage (before point (e)), the
micromagnetic domain around the crack tip gently rotates and deviates from the horizontal state to balance the elastic strain. In the
most of the region of the plate, the absolute value of 𝑚2 slowly increases. During these cycles, the crack gradually propagates. Then
there is a steep ascend in the length-cycle curve (point (e) to (f)) when the micromagnetic domain in the plate begins to rotate in
a large area. In this situation, the external magnetic field makes the magnetization vectors around the crack perpendicular to the
initial direction (𝑚1 = 1). The increase of |𝑚2 | holds the elastic strain energy. The degradation function 𝑓 (𝛼)
̄ around the crack tip
21
at point (e) is small, nearly 0, as is shown in Fig. 13(i). Thus the fracture toughness is damaged greatly and the crack propagates
easily. In the next stage (point (f) to (g)), the magnetization vectors in the crack path are gradually rotated to the direction of the
external field. Under this situation, 𝑚2 decreases to 0. The magnetostrictive strain makes elastic energy low. From Fig. 13(j), the
value of the degradation function 𝑓 (𝛼) ̄ is approximately 0.2 and thus the fatigue resistance is relatively high, leading to a extremely
slow crack growth and a plateau from point (f) to (g) in Fig. 13(a). The accumulation of positive elastic energy along the crack tip
needs more cycles and then 𝑓 (𝛼) ̄ will degrade closely to 0, shown in Fig. 13(k). Finally (after point (g)), 𝑚1 in the majority region of
the plate comes to −1 and the magnetization vector is parallel to the magnetic field. The crack continues to propagate until failure.
7. Conclusion
In summary, a thermodynamically consistent and micromagnetic-mechanically coupled phase-field model is developed for the
crack propagation of fracture and fatigue fracture in magnetostrictive alloys, in which the micromagnetic structure evolution
and crack topology are simultaneously involved. The microforce theory, thermodynamic laws, and the Coleman–Noll procedure
are applied to derive Allen–Cahn equation for the crack phase field order parameter and Landau–Lifshitz–Gilbert equation for
the magnetization-vector order parameter. The phase-field model is further extended to fatigue crack propagation. In the one-
dimensional case, the pure mechanical and micromagnetic-mechanical driving forces for crack have been analyzed. Under the pure
mechanical driving force, the maximum stresses under the different magnetic fields are the same. Nevertheless, the micromagnetic-
mechanical driving force under the magnetic field will enhance the critical stress. With the motivation of that the fracture behavior is
directly related to mechanical incentives, the pure mechanical driving force is chosen. In an amorphous Terfenol-D film as a typical
giant magnetostrictive material, it is shown that the mechanical behavior is affected by the external magnetic fields in the type-I and
three-point bending fracture. The magnetization vector is switched under the coupled magnetic-mechanical effects and the critical
load under the negative magnetic fields decreases. In the type-II fracture, the shear load is not sensitive to the magnetic fields. The
elliptical inclusion in the notched plate influences the magnetic domain evolution and magnetic fields decrease the critical load. The
larger fracture toughness of the elliptical inclusion makes the crack propagate around the ellipse and the multi-domain structure
drives the crack to nucleate in the tri-junction of magnetic domains. The fatigue cracking behavior depicts that the −1 T magnetic
field could increase the fatigue life cycles by 50%. Under this situation, the magnetization slowly changes through all the cycles
and the magnetostrictive strain at crack front delays the fatigue degradation of fracture toughness, thus slowing down the fatigue
crack growth. In sum, the phase-field model has been demonstrated capable of simulating the coupled behavior of micromagnetic
structure and static/fatigue crack propagation under external mechanical and magnetic fields. The simulation results on fatigue
behavior also imply a feasible route to improve the fatigue life of magnetostrictive alloys by regulating the micromagnetic structure
via external magnetic fields.
CRediT authorship contribution statement
Shen Sun: Writing – original draft, Visualization, Investigation, Formal analysis, Data curation, Conceptualization. Qihua Gong:
Writing – review & editing, Writing – original draft, Visualization, Supervision, Resources, Investigation, Funding acquisition, Formal
analysis, Data curation, Conceptualization. Yong Ni: Writing – review & editing, Writing – original draft, Supervision, Investigation,
Funding acquisition, Conceptualization. Min Yi: Writing – review & editing, Writing – original draft, Supervision, Resources, Project
administration, Investigation, Funding acquisition, Conceptualization.
Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared
to influence the work reported in this paper.
Data availability
Data will be made available on request.
Acknowledgments
The authors acknowledge the support from National Natural Science Foundation of China (12302134, 12272173, 11902150,
12025206), the Fundamental Research Funds for the Central Universities, China (NE2024001, NS2023054), National Youth Talents
Program of China, the Research Fund of State Key Laboratory of Mechanical and Control for Aerospace Structure, China, the Strategic
Priority Research Program of the Chinese Academy of Sciences, China (Grant No. XDB0620101), and a project Funded by the
Priority Academic Program Development of Jiangsu Higher Education Institutions, China. This work is partially supported by High
Performance Computing Platform of Nanjing University of Aeronautics and Astronautics, China.
22
References
Abdollahi, A., Arias, I., 2011. Phase-field modeling of the coupled microstructure and fracture evolution in ferroelectric single crystals. Acta Mater. 59 (12),
4733–4746. http://dx.doi.org/10.1016/j.actamat.2011.03.030.
Aguilar-Arteaga, K., Rodriguez, J., Barrado, E., 2010. Magnetic solids in analytical chemistry: a review. Anal. Chim. Acta 674 (2), 157–165. http://dx.doi.org/
10.1002/chin.201048276.
Alessi, R., Marigo, J.J., Maurini, C., Vidoli, S., 2018a. Coupling damage and plasticity for a phase-field regularisation of brittle, cohesive and ductile fracture:
One-dimensional examples. Int. J. Mech. Sci. 149, 559–576. http://dx.doi.org/10.1016/j.ijmecsci.2017.05.047.
Alessi, R., Vidoli, S., De Lorenzis, L., 2018b. A phenomenological approach to fatigue with a variational phase-field model: The one-dimensional case. Eng. Fract.
Mech. 190, 53–73. http://dx.doi.org/10.1016/j.engfracmech.2017.11.036.
Ambati, M., Gerasimov, T., De Lorenzis, L., 2015. Phase-field modeling of ductile fracture. Comput. Mech. 55 (5), 1017–1040. http://dx.doi.org/10.1007/s00466-
015-1151-4.
Amor, H., Marigo, J.J., Maurini, C., 2009. Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments. J. Mech.
Phys. Solids 57 (8), 1209–1229. http://dx.doi.org/10.1016/j.jmps.2009.04.011.
Beveridge, J.S., Stephens, J.R., Williams, M.E., 2011. The use of magnetic nanoparticles in analytical chemistry. Annu. Rev. Anal. Chem. 4 (1), 251–273.
http://dx.doi.org/10.1146/annurev-anchem-061010-114041.
Borden, M.J., Hughes, T.J., Landis, C.M., Anvari, A., Lee, I.J., 2016. A phase-field formulation for fracture in ductile materials: Finite deformation balance law
derivation, plastic degradation, and stress triaxiality effects. Comput. Methods Appl. Mech. Engrg. 312, 130–166. http://dx.doi.org/10.1016/j.cma.2016.09.
005.
Borden, M.J., Verhoosel, C.V., Scott, M.A., Hughes, T.J., Landis, C.M., 2012. A phase-field description of dynamic brittle fracture. Comput. Methods Appl. Mech.
Engrg. 217–220, 77–95. http://dx.doi.org/10.1016/j.cma.2012.01.008.
Bourdin, B., Marigo, J.J., Maurini, C., Sicsic, P., 2014. Morphogenesis and propagation of complex cracks induced by thermal shocks. Phys. Rev. Lett. 112 (1),
014301. http://dx.doi.org/10.1103/PhysRevLett.112.014301.
Brown, W., 1963. Micromagnetics. Wiley.
Carrara, P., Ambati, M., Alessi, R., De Lorenzis, L., 2020. A framework to model the fatigue behavior of brittle materials based on a variational phase-field
approach. Comput. Methods Appl. Mech. Engrg. 361, 112731. http://dx.doi.org/10.1016/j.cma.2019.112731.
Clark, A., 1980. Magnetostrictive rare earth-Fe2 compounds. In: Handbook of Ferromagnetic Materials. Elsevier, pp. 531–589. http://dx.doi.org/10.1016/s1574-
9304(05)80122-1, chapter 7.
Clatterbuck, D.M., Chan, J.W., John William Morris, J., 2000. The influence of a magnetic field on the fracture toughness of ferromagnetic steel. Mater. Trans.,
JIM 41 (8), 888–892. http://dx.doi.org/10.2320/matertrans1989.41.888.
Colussi, M., Berto, F., Mori, K., Narita, F., 2016. Fracture behavior of cracked giant magnetostrictive materials in three-point bending under magnetic fields:
strain energy density criterion. Adv. Eng. Mater. 18 (12), 2063–2069. http://dx.doi.org/10.1002/adem.201500533.
Dang, L., He, X., Tang, D., Li, Y., Wang, T., 2022. A fatigue life prediction approach for laser-directed energy deposition titanium alloys by using support vector
regression based on pore-induced failures. Int. J. Fatigue 159, 106748. http://dx.doi.org/10.1016/j.ijfatigue.2022.106748.
Eerenstein, W., Mathur, N., Scott, J.F., 2006. Multiferroic and magnetoelectric materials. Nature 442 (7104), 759–765. http://dx.doi.org/10.1038/NATURE05023.
Fang, D.N., Wan, Y.P., Soh, A.K., 2004. Magnetoelastic fracture of soft ferromagnetic materials. Theor. Appl. Fract. Mech. 42 (3), 317–334. http://dx.doi.org/
10.1016/j.tafmec.2004.09.006.
Fidler, J., Schrefl, T., 2000. Micromagnetic modelling – the current state of the art. J. Phys. D: Appl. Phys. 33 (15), R135. http://dx.doi.org/10.1088/0022-
3727/33/15/201.
Fischbacher, T., Franchin, M., Bordignon, G., Fangohr, H., 2007. A systematic approach to multiphysics extensions of finite-element-based micromagnetic
simulations: Nmag. IEEE Trans. Magn. 43 (6), 2896–2898. http://dx.doi.org/10.1109/tmag.2007.893843.
Francfort, G.A., Marigo, J.J., 1998. Revisiting brittle fracture as an energy minimization problem. J. Mech. Phys. Solids 46 (8), 1319–1342. http://dx.doi.org/
10.1016/S0022-5096(98)00034-9.
Gilbert, T.L., 2004. A phenomenological theory of damping in ferromagnetic materials. IEEE Trans. Magn. 40 (6), 3443–3449. http://dx.doi.org/10.1109/TMAG.
2004.836740.
Griffith, A., 1920. The phenomena of rupture and flow in solids. Philos. Trans. R. Soc. Lond. A221, 163–197. http://dx.doi.org/10.1098/rsta.1921.0006.
Grönefeld, M., Kronmüller, H., 1989. Calculation of strayfields near grain edges in permanent magnet material. J. Magn. Magn. Mater. 80 (2–3), 223–228.
http://dx.doi.org/10.1016/0304-8853(89)90122-4.
Gurtin, M.E., 1996. Generalized Ginzburg-Landau and Cahn-Hilliard equations based on a microforce balance. Physica D 92 (3–4), 178–192. http://dx.doi.org/
10.1016/0167-2789(95)00173-5.
Hubert, A., Schäfer, R., 1998. Magnetic Domains: The Analysis of Magnetic Microstructures. Springer, Berlin-Heidelberg-New York.
Jiles, D., 2003. Recent advances and future directions in magnetic materials. Acta Mater. 51 (19), 5907–5939. http://dx.doi.org/10.1016/j.actamat.2003.08.011.
Kovacs, A., Exl, L., Kornell, A., Fischbacher, J., Hovorka, M., Gusenbauer, M., Breth, L., Oezelt, H., Praetorius, D., Suess, D., Schrefl, T., 2022. Magnetostatics
and micromagnetics with physics informed neural networks. J. Magn. Magn. Mater. 548, 168951. http://dx.doi.org/10.1016/j.jmmm.2021.168951.
Kronmüller, H., Fähnle, M., 2003. Micromagnetism and the Microstructure of Ferromagnetic Solids. Cambridge University Press, Cambridge, UK.
Kuhn, C., Schlüter, A., Müller, R., 2015. On degradation functions in phase field fracture models. Comput. Mater. Sci. 108, 374–384. http://dx.doi.org/10.1016/
j.commatsci.2015.05.034.
Landau, L., Lifshitz, E., 1935. On the theory of the dispersion of magnetic permeability in ferromagnetic bodies. Phys. Z. Sowjetunion 8 (1), 153–169.
http://dx.doi.org/10.1063/1.4931784.
Landis, C.M., 2008. A continuum thermodynamics formulation for micro-magnetomechanics with applications to ferromagnetic shape memory alloys. J. Mech.
Phys. Solids 56 (10), 3059–3076. http://dx.doi.org/10.1016/j.jmps.2008.05.004.
Lei, N., Devolder, T., Agnus, G., Aubert, P., Daniel, L., Kim, J.-V., Zhao, W., Trypiniotis, T., Cowburn, R.P., Chappert, C., Ravelosona, D., Lecoeur, P., 2013.
Strain-controlled magnetic domain wall propagation in hybrid piezoelectric/ferromagnetic structures. Nature Commun. 4 (1), 1378. http://dx.doi.org/10.
1038/ncomms2386.
Liang, W., Shen, Y., Zhao, M., 2000. Magnetoelastic formulation of soft ferromagnetic elastic problems with collinear cracks: energy density fracture criterion.
Theor. Appl. Fract. Mech. 34 (1), 49–60. http://dx.doi.org/10.1016/S0167-8442(00)00023-9.
Lo, Y.S., Borden, M.J., Ravi-Chandar, K., Landis, C.M., 2019. A phase-field model for fatigue crack growth. J. Mech. Phys. Solids 132, 103684. http:
//dx.doi.org/10.1016/j.jmps.2019.103684.
Maugin, G.A., 1995. Material forces: concepts and applications. Appl. Mech. Rev. 48 (5), 213–245. http://dx.doi.org/10.1115/1.3005101.
Maugin, G.A., Epstein, M., Trimarco, C., 1992. Theory of elastic inhomogeneities in electromagnetic materials. Internat. J. Engrg. Sci. 30 (10), 1441–1449.
http://dx.doi.org/10.1016/0020-7225(92)90154-9.
Mei, J., Xing, S., Vasu, A., Chung, J., Desai, R., Dong, P., 2020. The fatigue limit prediction of notched components – A critical review and modified stress
gradient based approach. Int. J. Fatigue 135, 105531. http://dx.doi.org/10.1016/j.ijfatigue.2020.105531.
23
Miehe, C., Ethiraj, G., 2012. A geometrically consistent incremental variational formulation for phase field models in micromagnetics. Comput. Methods Appl.
Mech. Engrg. 245–246, 331–347. http://dx.doi.org/10.1016/j.cma.2012.03.021.
Miehe, C., Hofacker, M., Schänzel, L.M., Aldakheel, F., 2015a. Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-to-ductile
failure criteria and crack propagation in thermo-elastic-plastic solids. Comput. Methods Appl. Mech. Engrg. 294, 486–522. http://dx.doi.org/10.1016/j.cma.
2014.11.017.
Miehe, C., Hofacker, M., Welschinger, F., 2010a. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on
operator splits. Comput. Methods Appl. Mech. Engrg. 199 (45), 2765–2778. http://dx.doi.org/10.1016/j.cma.2010.04.011.
Miehe, C., Schänzel, L.M., Ulmer, H., 2015b. Phase field modeling of fracture in multi-physics problems. Part I. Balance of crack surface and failure criteria for
brittle crack propagation in thermo-elastic solids. Comput. Methods Appl. Mech. Engrg. 294, 449–485. http://dx.doi.org/10.1016/j.cma.2014.11.016.
Miehe, C., Welschinger, F., Hofacker, M., 2010b. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE
implementations. Internat. J. Numer. Methods Engrg. 83 (10), 1273–1311. http://dx.doi.org/10.1002/nme.2861.
Moreno-Mateos, M.A., Hossain, M., Steinmann, P., Garcia-Gonzalez, D., 2023. Hard magnetics in ultra-soft magnetorheological elastomers enhance fracture
toughness and delay crack propagation. J. Mech. Phys. Solids 173, 105232. http://dx.doi.org/10.1016/j.jmps.2023.105232.
Moreno-Mateos, M.A., Mehnert, M., Steinmann, P., 2024. Electro-mechanical actuation modulates fracture performance of soft dielectric elastomers. Internat. J.
Engrg. Sci. 195, 104008. http://dx.doi.org/10.1016/j.ijengsci.2023.104008.
Moreno-Mateos, M.A., Steinmann, P., 2024. Configurational force method enables fracture assessment in soft materials. J. Mech. Phys. Solids 186, 105602.
http://dx.doi.org/10.1016/j.jmps.2024.105602.
Narita, F., Shikanai, K., Shindo, Y., Mori, K., 2016. Three-point bending fracture behavior of cracked giant magnetostrictive materials under magnetic fields. J.
Test. Eval. 44 (4), 1454–1460. http://dx.doi.org/10.1520/jte20140439.
Nguyen, T.T., Yvonnet, J., Bornert, M., Chateau, C., Sab, K., Romani, R., Le Roy, R., 2016. On the choice of parameters in the phase field method for simulating
crack initiation with experimental validation. Int. J. Fract. 197 (2), 213–226. http://dx.doi.org/10.1007/s10704-016-0082-1.
Noll, W., Coleman, B.D., 1974. The thermodynamics of elastic materials with heat conduction and viscosity. In: The Foundations of Mechanics and
Thermodynamics. Springer, pp. 145–156. http://dx.doi.org/10.1007/BF01262690.
Pao, Y.-H., Yeh, C.-S., 1973. A linear theory for soft ferromagnetic elastic solids. Internat. J. Engrg. Sci. 11 (4), 415–436. http://dx.doi.org/10.1016/0020-
7225(73)90059-1.
Peron, M., Katabira, K., Viespoli, L.M., Narita, F., Berto, F., 2019. Mixed mode fracture behavior of notched giant magnetostrictive: mechanical characterization
and comparison among failure criteria. Theor. Appl. Fract. Mech. 99, 194–204. http://dx.doi.org/10.1016/j.tafmec.2018.12.005.
Schreiber, C., Kuhn, C., Müller, R., Zohdi, T., 2020. A phase field modeling approach of cyclic fatigue crack growth. Int. J. Fract. 225 (1), 89–100.
http://dx.doi.org/10.1007/s10704-020-00468-w.
Seiler, M., Linse, T., Hantschke, P., Kästner, M., 2020. An efficient phase-field model for fatigue fracture in ductile materials. Eng. Fract. Mech. 224, 106807.
http://dx.doi.org/10.1016/j.engfracmech.2019.106807.
Shi, Y., Yu, H., Shimada, T., Wang, J., Kitamura, T., 2017. Phase field simulations on domain switching-induced toughening in ferromagnetic materials. Eur. J.
Mech. A Solids 65, 205–211. http://dx.doi.org/10.1016/j.euromechsol.2017.04.007.
Shindo, Y., 1977. The linear magnetoelastic problem for a soft ferromagnetic elastic solid with a finite crack. J. Appl. Mech. 44 (1), 47–50. http://dx.doi.org/
10.1115/1.3424012.
Shindo, Y., 1978. Magnetoelastic interaction of a soft ferromagnetic elastic solid with a penny-shaped crack in a constant axial magnetic field. J. Appl. Mech.
45 (2), 291–296. http://dx.doi.org/10.1115/1.3424290.
Shindo, Y., 1980. Singular stresses in a soft ferromagnetic elastic solid with two coplanar Griffith cracks. Int. J. Solids Struct. 16 (6), 537–543. http:
//dx.doi.org/10.1016/0020-7683(80)90004-9.
Shindo, Y., Horiguchi, K., Shindo, T., 1999. Magneto-elastic analysis of a soft ferromagnetic plate with a through crack under bending. Internat. J. Engrg. Sci.
37 (6), 687–702. http://dx.doi.org/10.1016/S0020-7225(98)00095-0.
Shindo, Y., Komatsu, T., Narita, F., Horiguchi, K., 2006. Magnetic stress intensity factor for an edge crack in a soft ferromagnetic elastic half-plane under tension.
Acta Mech. 182 (3–4), 183–193. http://dx.doi.org/10.1007/s00707-005-295-2.
Shindo, Y., Ohnishi, I., Tohyama, S., 1997. Flexural wave scattering at a through crack in a conducting plate under a uniform magnetic field. J. Appl. Mech. 66
(4), 828–834. http://dx.doi.org/10.1115/1.2788988.
Shindo, Y., Sekiya, D., Narita, F., Hohiguchi, K., 2004. Tensile testing and analysis of ferromagnetic elastic strip with a central crack in a uniform magnetic
field. Acta Mater. 52 (15), 4677–4684. http://dx.doi.org/10.1016/j.actamat.2004.06.029.
Shindo, Y., Sekiya, D., Narita, F., Horiguchi, K., 2005. Magnetic stress intensity factors for two symmetric edge cracks in a soft ferromagnetic elastic strip under
tension. JSME Int. J. Ser. A 48 (1), 7–13. http://dx.doi.org/10.1299/jsmea.48.7.
Shindo, Y., Shindo, I., Narita, F., 2008. The influence of magnetic field on the fracture toughness of soft ferromagnetic materials. Eng. Fract. Mech. 75 (10),
3010–3019. http://dx.doi.org/10.1016/j.engfracmech.2007.12.015.
Shindo, Y., Tohyama, S., 1998. Scattering of oblique flexural waves by a through crack in a conducting mindlin plate in a uniform magnetic field. Int. J. Solids
Struct. 35 (17), 2183–2203. http://dx.doi.org/10.1016/S0020-7683(97)00183-2.
Simoes, M., Martínez-Pañeda, E., 2021. Phase field modelling of fracture and fatigue in Shape Memory Alloys. Comput. Methods Appl. Mech. Engrg. 373, 113504.
http://dx.doi.org/10.1016/j.cma.2020.113504.
Sridhar, A., Keip, M.A., 2019. A phase-field model for anisotropic brittle fracturing of piezoelectric ceramics. Int. J. Fract. 220 (2), 221–242. http://dx.doi.org/
10.1007/s10704-019-00391-9.
Steinke, C., Kaliske, M., 2019. A phase-field crack model based on directional stress decomposition. Comput. Mech. 63 (5), 1019–1046. http://dx.doi.org/10.
1007/s00466-018-1635-0.
Szambolics, H., Buda-Prejbeanu, L., Toussaint, J., Fruchart, O., 2008. A constrained finite element formulation for the Landau–Lifshitz–Gilbert equations. Comput.
Mater. Sci. 44 (2), 253–258. http://dx.doi.org/10.1016/j.commatsci.2008.03.019.
Tang, W., Yi, M., Chen, L.-Q., Guo, W., 2024. Classical fatigue theory informed phase-field model for high-cycle fatigue life and fatigue crack growth. Eng. Fract.
Mech. 306, 110212. http://dx.doi.org/10.1016/j.engfracmech.2024.110212.
Truesdell, C., Noll, W., 2004. The Non-Linear Field Theories of Mechanics. Springer, http://dx.doi.org/10.1007/978-3-662-10388-3_1.
Wang, J., Zhang, J., 2013. A real-space phase field model for the domain evolution of ferromagnetic materials. Int. J. Solids Struct. 50 (22–23), 3597–3609.
http://dx.doi.org/10.1016/j.ijsolstr.2013.07.001.
Xu, B.-X., Schrade, D., Gross, D., Mueller, R., 2010. Fracture simulation of ferroelectrics based on the phase field continuum and a damage variable. Int. J. Fract.
166 (1), 163–172. http://dx.doi.org/10.1007/s10704-010-9520-7.
Yamaguchi, Y., Horiguchi, K., Shindo, Y., Sekiya, D., Kumagai, S., 2003. Fracture and deformation properties of Ni–Fe superalloy in cryogenic high magnetic
field environments. Cryogenics 43 (8), 469–475. http://dx.doi.org/10.1016/s0011-2275(03)00123-1.
Yi, M., Xu, B.-X., 2014. A constraint-free phase field model for ferromagnetic domain evolution. Proc. R. Soc. A 470 (2171), 20140517. http://dx.doi.org/10.
1098/rspa.2014.0517.
Yi, M., Xue, M., Cong, P., Song, Y., Zhang, H., Wang, L., Zhou, L., Li, Y., Guo, W., 2024. Machine learning for predicting fatigue properties of additively
manufactured materials. Chin. J. Aeronaut. 37 (4), 1–22. http://dx.doi.org/10.1016/j.cja.2023.11.001.
24
Zheng, X.-J., Liu, X.-E., 2005. A nonlinear constitutive model for Terfenol-D rods. J. Appl. Phys. 97 (5), 053901. http://dx.doi.org/10.1063/1.1850618.
Zhou, Y.-H., Zheng, X.-J., 1997. A general expression of magnetic force for soft ferromagnetic plates in complex magnetic fields. Internat. J. Engrg. Sci. 35 (15),
1405–1417. http://dx.doi.org/10.1016/S0020-7225(97)00051-7.
Zhou, H.-M., Zhou, Y.-H., Zheng, X.-J., 2008. A general theoretical model of magnetostrictive constitutive relationships for soft ferromagnetic material rods. J.
Appl. Phys. 104 (2), 023907. http://dx.doi.org/10.1063/1.2957075.
Zhou, H.-M., Zhou, Y.-H., Zheng, X.-J., Ye, Q., Wei, J., 2009. A general 3-D nonlinear magnetostrictive constitutive model for soft ferromagnetic materials. J.
Magn. Magn. Mater. 321 (4), 281–290. http://dx.doi.org/10.1016/j.jmmm.2008.09.012.
Zhu, B., Lo, C.C.H., Lee, S.J., Jiles, D., 2001. Micromagnetic modeling of the effects of stress on magnetic properties. J. Appl. Phys. 89 (11), 7009–7011.
http://dx.doi.org/10.1063/1.1363604.
25

1-s2.0-S0022509624002333-main

Uploaded by

Copyright:

Available Formats

1-s2.0-S0022509624002333-main

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

1-s2.0-S0022509624002333-main

Uploaded by

Copyright:

Available Formats

J. Mech. Phys.

Solids 191 (2024) 105767

Contents lists available at ScienceDirect

Journal of the Mechanics and Physics of Solids

A micromagnetic-mechanically coupled phase-field model for

ARTICLE INFO ABSTRACT

2. Micromagnetic-mechanically coupled phase-field model for fracture

2.1. Balance law

• Balance of linear momentum

𝜎𝑖𝑗,𝑗 + 𝜌𝑏𝑖 = 0 in 𝛺 , (4)

𝑢𝑖 = 𝑢̂ 𝑖 on 𝜕𝛺𝑢 , 𝜎𝑗𝑖 𝑛𝑗 = 𝑡̂𝑖 on 𝜕𝛺𝜎 (5)

𝜎𝑖𝑗 = 𝜎𝑗𝑖 . (6)

• Gauss’s law for magnetic flux

𝜖𝑖𝑗𝑘 𝐸𝑘,𝑗 = −𝐵̇ 𝑖 (8)

𝜖𝑖𝑗𝑘 𝐻𝑘,𝑗 = 𝐽𝑖 . (9)

The magnetic boundary conditions are described as

𝐵𝑖 𝑛𝑖 = 𝐵̂ on 𝜕𝛺𝐵 , 𝜙 = 𝜙̂ on 𝜕𝛺𝜙 , 𝜖𝑖𝑗𝑘 𝐻𝑗 𝑛𝑘 = 𝐾̂ 𝑖 on 𝜕𝛺𝐾 (10)

• Balance of microforce associated with magnetization vector 𝑀𝑖

From Eqs. (13) and (15), one yields

• Balance of microforce associated with crack phase field 𝑐

𝑒̇ = 𝑢̇ 𝑗,𝑖 𝜎𝑖𝑗 + 𝐵̇ 𝑖 𝐻𝑖 − 𝜋 𝑐 𝑐̇ + 𝜉𝑖𝑐 𝑐̇ ,𝑖 + 𝜉𝑖𝑗,𝑖 𝑀̇ 𝑗 + 𝜉𝑖𝑗 𝑀̇ 𝑗,𝑖 + 𝑓𝑖 𝑀̇ 𝑖 . (21)

2.2. Constitutive relations and evolution equations

• Second law of thermodynamics

• Free energy imbalance

A governing equation for the evolution of the magnetization vector reads

2.3. Free energy density

ℎ = ℎexc + ℎani + ℎent1 + ℎent2 + ℎe + ℎc . (41)

ℎexc = 𝑔 (𝑐, 𝜁) 𝐴e ‖∇𝐦‖2 = 𝑔 (𝑐, 𝜁) 𝐴e 𝑚𝑖,𝑗 𝑚𝑖,𝑗 (42)

where 𝐞a is the easy-axis vector. ℎent1 is expressed as

ℎent1 = −𝑔 (𝑐, 𝜁) 𝜇0 𝑀s 𝐇 ⋅ 𝐦 = 𝑔 (𝑐, 𝜁) 𝜇0 𝑀s 𝜙,𝑖 𝑚𝑖 , (45)

and ℎent2 is expressed as

The elastic energy density ℎe

2.4. Split of the elastic energy density

𝜀𝑖𝑗 = 𝜀e𝑖𝑗 + 𝜀0𝑖𝑗 . (50)

𝜀0𝑖𝑗 of a cubic crystal is expressed as

𝜆100 = 𝜆111 = 𝜆s . (52)

𝜺e− = 𝐏Λ− 𝐏T (57)

2.5. Irreversible condition for crack

 (𝑡) = max ℎe+ (𝑡) . (61)

Then the governing equation for 𝑐 is derived as

3. Extension to phase-field fatigue

In the phase field model of fatigue, a degradation function 𝑓 (𝛼)

The critical stress and the corresponding strain are found as

0= 𝐵𝑖 𝛷,𝑖𝜙 d𝑉 − 𝐵𝑖 𝑛𝑖 𝛷𝜙 d𝑆, (75)

𝑅𝐼𝑢 = 𝜎𝑖𝑗 𝑁,𝑗𝐼 d𝑉 − 𝜎𝑖𝑗 𝑛𝑗 𝑁 𝐼 d𝑆, (79)

𝑀𝑖,𝑗 𝑛𝑗 = 0 on 𝜕𝛺MN . (86)

Similarly, for the crack phase field, one has

𝑐,𝑖 𝑛𝑖 = 0 on 𝜕𝛺𝑐N . (87)

For the Dirichlet boundary conditions for 𝑀𝑖 and 𝑐, one has

The array of all the nodal variables is

𝐝𝐼 = [𝑢𝐼1 , 𝑢𝐼2 , 𝑢𝐼3 , 𝜙𝐼 , 𝑀1𝐼 , 𝑀2𝐼 , 𝑀3𝐼 , 𝑐 𝐼 , 𝜆𝐼𝑀 ]T . (90)

5. Crack driving force

5.1. Evolution of fracture in case of pure mechanical driving force

The critical state 𝜏 = 0 gives

5.2. Evolution of fracture in case of micromagnetic-mechanical driving force

5.3. Comparison of the results in the two cases of driving force

6. Simulation results and discussions

6.1. Monotonic tension under horizontal magnetic field

6.2. Monotonic tension under vertical magnetic field

6.3. Pure shear under magnetic field