Phase-Field Modeling of Fracture for Ferromagnetic

Materials through Maxwell’s Equation

Nima Noiia,b,111Corresponding author (Nima Noii).
E-mail addresses<EMAIL_ADDRESS>(N. Noii<EMAIL_ADDRESS>freiberg.de (M. Ghasabeh<EMAIL_ADDRESS>(P. Wriggers). ,
Mehran Ghasabehc, Peter Wriggersb,d

a Deutsches Institut für Kautschuktechnologie (DIK e.V.)

Eupener Straße 33, 30519 Hannover, Germany

b Institute of Continuum Mechanics

Leibniz Universität Hannover, An der Universität 1, 30823 Garbsen, Germany

c Chair of Soil Mechanics and Foundation Engineering

Technische Universität Bergakademie Freiberg, 09599 Freiberg, Germany

d Cluster of Excellence PhoenixD (Photonics, Optics, and Engineering -
Innovation

Across Disciplines), Leibniz Universität Hannover, Germany

Abstract

Electro-active materials are classified as electrostrictive and piezoelectric
materials. They deform under the action of an external electric field.
Piezoelectric material, as a special class of active materials, can produce an
internal electric field when subjected to mechanical stress or strain. In
return, there is the converse piezoelectric response, which expresses the
induction of the mechanical deformation in the material when it is subjected
to the application of the electric field. This work presents a variational-
based computational modeling approach for failure prediction of ferromagnetic
materials. In order to solve this problem, a coupling between magnetostriction
and mechanics is modeled, then the fracture mechanism in ferromagnetic
materials is investigated. Furthermore, the failure mechanics of ferromagnetic
materials under the magnetostrictive effects is studied based on a variational
phase-field model of fracture. Phase-field fracture is numerically challenging
since the energy functional may admit several local minima, imposing the
global irreversibility of the fracture field and dependency of regularization
parameters related discretization size. Here, the failure behavior of a
magnetoelastic solid body is formulated based on the Helmholtz free energy
function, in which the strain tensor, the magnetic induction vector, and the
crack phase-field are introduced as state variables. This coupled formulation
leads to a continuity equation for the magnetic vector potential through well-
known Maxwell’s equations. Hence, the energetic crack driving force is
governed by the coupled magneto-mechanical effects under the magneto-static
state. Several numerical results substantiate our developments.

Keywords: Maxwell’s equation, phase-field fracture, magnetization,
magnetostriction, Ferromagnetic, magnetic vector potential, electric field,
magnetic field, magnetomechanical.

###### Contents

  1. 1 Introduction
  2. 2 Phase-field modeling of fracture for ferromagnetic materials
    1. 2.1 Primary fields for the multi-field problem
    2. 2.2 Governing equations of the failure mechanism for magnetostrictive effects
      1. 2.2.1 Elastic contribution
      2. 2.2.2 Electromagnetic contribution
      3. 2.2.3 Fracture contribution
      4. 2.2.4 Transition structure from undamaged to fully damaged states
  3. 3 Energy quantities and variational principles
    1. 3.1 Constitutive functions
      1. 3.1.1 Energy quantities
      2. 3.1.2 Magnetization contribution
      3. 3.1.3 Magnetostrictive contribution
      4. 3.1.4 Elastic contribution
      5. 3.1.5 Fracture contribution
    2. 3.2 Variational formulations for the coupled multi-field problem
      1. 3.2.1 Electromagnetic contribution
      2. 3.2.2 Elastic contribution
      3. 3.2.3 Fracture contribution
        1. 3.3 Space finite element discretization
          1. 4 Numerical examples
            1. 4.1 Example 1: Magnetostatic problem for transient magnetic vector potential
            2. 4.2 Example 2: Magneto-mechanically induced cracking in an iron beam surrounded by vacuum
            3. 4.3 Example 3: Magneto-mechanically induced cracking in a ferromagnetic material containing predefined notches and wires
              1. 4.3.1 Sub-example 1: Transversely wired plate with three winding wires
              2. 4.3.2 Sub-example 2: Transversely wired plate with nine winding wires
              3. 4.3.3 Sub-example 3: Longitudinally wired plate with four winding wires
              4. 4.3.4 Sub-example 4: Longitudinally wired plate with five winding wires
                1. 5 Conclusion

##  1 .Introduction

Active materials are a particular type of material that undergoes mechanical
deformation in response to external effects. The sources of these effects can
be pressure, thermal, electric, and magnetic fields [1]. In the literature,
these materials are known as smart materials. There are various types of
active materials, such as shape memory alloys, electrostrictive elastomers,
piezoelectric materials, ferroelectric materials, and electro- and magneto-
active polymers. [2, 3]. The basic mechanical properties of active materials,
such as strain generation capability, stiffness, strain, Hysteresis and
electrical impedance vary widely. The energy output of active materials used
in actuators is principally dominated by stiffness and the amount of strain
energy generated by the material [1]. The materials can be classified
according to their response when they are subjected to external stimuli. To
this end, active materials can be generally categorized as electro-active and
magneto-active materials.

Electro-active materials are classified as electrostrictive, and piezoelectric
materials. These materials deform mechanically under an external electric
field. Piezoelectric materials can produce an internal electric field when
subjected to mechanical stress or strain (direct piezoelectricity) [4].
Piezoelectric materials also represent the reverse piezoelectric effect
(converse piezoelectricity), where stress or strain is generated when they are
subjected to an electric field [5]. The direct piezoelectric response of
material is determined by the conversion of the mechanical into the electrical
energy; however, the converse piezoelectric response expresses the induction
of mechanical deformation as a consequence of the applied electric field [6].

Figure 1: Schematic representation of Joule effect (direct magnetostrictive
effect) in the first raw and Villari effect (reversible magnetostrictive
effect) in the second raw.

Piezoelectric materials are applied within a great range, including
piezoelectric motors, actuators in industrial sector, sensors in the medical
sector, actuators in consumer electronics (printer, speakers), piezoelectric
buzzers, piezoelectric igniters, microphones, nanopositioning in atomic force
microscope (AFM) and the scanning tunneling microscope (STM), and micro-
robotics (defense) [7].

Electrostrictive materials generate stress or strain when subjected to an
external electric field [1]. However, there is a principal difference between
piezoelectric and electrostrictive materials. The piezoelectric effect is
possible only in non-centrosymmetric materials. However, the electrostrictive
effect is not limited by symmetry and is present in all materials, even those
that are amorphous. Therefore, the electrostrictive effect exhibits a
nonlinear (second-order or quadratic) dependency of the strain on the applied
electric field [6].

Magnetostriction, a key feature of magnetic materials, is defined as the
alteration in shape and dimension of the material during the magnetization
process. It continues till the magnetic saturation of the material is
attained. The magnetostriction phenomenon was first introduced in the work of
Joule [8]. It is also known as the Joule effect. In 1846, Villari discovered
the reverse effect of magnetization (see Figure 1)[9]. Ferromagnetic material
and ferri-magnetic are known for their magnetostrictive characteristics. A
magnetostrictive material consists of tiny fragments. These fragments, usually
iron, nickel, or cobalt, have small magnetic moments as a result of their ”3d”
shells that are not completely filled with electrons. Ferromagnetics
fundamentally acts like tiny permanent bar magnets. Magnetic materials have a
property of magnetostriction that causes them to change their shape or
dimensions during the magnetization process. Such materials can convert
magnetic energy into kinetic or reverse and thus are mainly used to construct
actuators and sensors. The applicability of the materials is quantified by the
magnetostrictive coefficient. It is defined as the fractional change in length
as the magnetization of the material increases from zero to the saturation
value. This coefficient can be positive and negative. The magnetostriction
characteristic induces energy loss due to frictional heating in susceptible
ferromagnetic cores. It also produces the low-pitched humming sound that can
be heard coming from transformers, where an existing magnetic field is changed
by oscillating AC currents.

In the literature, computational approaches to modeling the electromagnetic
response of the electro- and magneto-active materials are developed within a
framework of thermodynamics. For some active-material, it is sufficient to
take into account only a quasi-static state for modeling the magnetic response
[10]. The work proposes a hybrid SBFEM-FEM approach for efficiently
calculating stray magnetic fields in unbounded domains. Most of the works aim
at constructing a coupling between electromagnetics and mechanics. The key
feature of these approaches is to investigate the stress response of the
materials under the electro- and magneto-striction effects. The coupled
formulation within the framework of magnetoelasticity is found in the work of
Brigadnov and Dorfmann [11] where a mathematical model is developed to
investigate the mechanical characteristics of the magneto-sensitive elastomer
under the application of a magnetic field. Further studies related to magneto-
sensitive elastomers are conducted in the work of Dorfmann, and his colleagues
[12, 13]. The investigation of an electromagnetic forming process is conducted
in [14] by formulating a fully coupled electromagnetic-thermomechanical model.
A computational approach incorporating the energy-based magneto-mechanical
model is developed in [15, 16, 17] for the electric electrical steel sheets.
An energy-density functional method is proposed in [18] based on an isotropic
spline-based thermodynamic approach to model magneto-mechanical behavior in
ferromagnetic material.

Electro-magneto-mechanical coupling effects are formulated, based on
incremental variational principles, within a general continuum mechanics
framework and numerical implementation of dissipative functional materials in
[19]. Micromagnetics is another concept where a geometrically consistent
incremental variational formulation is extended in a micro-magneto-elastic
model that accounts for the micro-structural evolution of both magnetically-
and mechanically-driven magnetic domains in ferromagnetic materials [20]. The
variational-based computational modeling of ferromagnetics and
magnetorheological elastomers, which present a mutually-coupled magneto-
mechanical response, is proposed in [21]. Further study on the behavior of
magnetosensitive polymers based on a multiplicative magneto-elastic model
within a micromechanically-based framework is found in [22]. Moreover, models
regarding the thermomechanical effects are developed to investigate the
coupling between electromagnetics and thermomechanics in [23] for electric
motors. The multi-physics problem regarding magnetic-thermal-mechanical
modeling of the electromagnetic rail launching is investigated in [24]. The
principles of designing an in situ fatigue testing system, including
structural resonances, transient response, grip design, and thermal insulation
performance of the device, are studied in detail in [25] by developing a
coupled thermo-mechanical model. The Cell-based smoothed finite element method
(CS-FEM) and the asymptotic homogenization method (AHM) is incorporated in
[26] to simulate the magneto-electro-elastic response of a structure under
dynamic load accurately.

In engineering applications, predicting crack initiation and propagation in
structures under mechanical loading and environmental conditions is greatly
important. Besides classical approaches for fracture mechanics, non-local
damage models and variational approaches have been developed in the
literature. The non-local damage models have been formulated to overcome the
ill-posedness due to spatial discretization in finite element simulations. The
variational approaches are introduced based on energy minimization [27, 28,
29], and their regularization is obtained by $\Gamma$-convergence, which is
fundamentally inspired by the work of image segmentation conducted by Mumford
and Shah [30]. The model is then improved by formulating a Ginzburg-Landau-
type evolution equation of the fracture phase-field [31]. In recent years, the
variationally-based phase-field approach to fracture within the
thermodynamically consistent framework is proposed in the work of Miehe et al.
[32]. In the latest work, the robust algorithmic formulation of the evolution
of diffuse crack topology in time is proposed by introducing a local history
field that determines the maximum tensile strain obtained in history. It then
acts as a crack driving force in the evolution of the crack phase field. In
the literature, few works investigate the material’s cracking response; for
example, in [33], the failure mechanism of piezoelectric is studied by
extending electromechanics to the phase-field approach. Ferroelastic ceramics
is another class of active material in which the crack initiation and
propagation are examined in [34], which proposes a computational framework
that regards the electric displacement saturation and its effect on the
hysteretic behavior of ferroelectric ceramics and the resulting cracking
response.

This work applies the phase-field approach to the magneto-mechanically driven
fracture in ferromagnetic materials. In the formulations for examining the
response of the electromagnetic materials, finite element method (FEM) is
applied to discretize the time-dependent Maxwell’s equations on a bounded
domain in three-dimensional space. At first, the formulation is developed by
deriving a weak formulation for the electric and magnetic fields with
approximate Neumann and Dirichlet boundary conditions, and the problem is
discretized both in time and space. In general, the electric and the magnetic
fields are discretized by adapting Nédélec curl-conforming and Raviart-Thomas
div-conforming finite element approaches, respectively.

The main objective of this contribution is to introduce:

  * •

Extension of Maxwell’s equations towards a coupled magneto-mechanical model to
investigate the stress response of the ferromagnetic materials;

  * •

Transition rule for the electromagnetic material properties from undamaged to
fully damaged states through the fracture phase-field, which acts as a
geometric interpolation variable;

  * •

Investigation of the magnetostrictive-induced cracking in the ferromagnetic
materials by developing a magneto-mechanical model coupled with the phase-
field model;

  * •

Representation of numerical examples to substantiate our developments in
predicting the fracture response of the ferromagnetic material.

The paper is structured as follows: For a better insight into phase-field
modeling of fracture for ferromagnetic materials, primary fields for the
multi-field problem are provided in Section 2. Through a regularized fracture
phase-field formulation, we further outline the theoretical framework for
magnetostrictive-induced fracture in ferromagnetic materials. The transition
rule for the electromagnetic material properties from undamaged to fully
damaged states is also discussed. Next, Section 3 outlines the constitutive
energy density for the magneto-mechanical formulation of fracturing solids,
which suffices to define the variational formulation setting. In Section 4,
three numerical simulations are performed to demonstrate the correctness of
our algorithmic developments. Finally, the paper is concluded with some
remarks.

##  2 .Phase-field modeling of fracture for ferromagnetic materials

This section outlines a mathematical model for magnetostrictive-induced
cracking in ferromagnetic materials by developing a magneto-mechanical model,
considering small deformations. The fracture process is modeled by employing
the well-developed phase-field formulation to resolve the sharp crack surface
topology in the regularized concept. We specifically elaborate three governing
equations to characterize the constitutive formulations for the mechanical
deformation, electromagnetic as well as the fracture phase-field. This coupled
formulation is derived by obtaining a continuity equation for the magnetic
vector potential through well-known Maxwell’s equations. Thus, the
comprehensive objective of the following section is to advance a continuum
theory of cracking response in ferromagnetic materials within a framework of
thermodynamics.

###  2.1 .Primary fields for the multi-field problem

Let $\Omega\subset{\mathbb{R}}^{\delta}$ be an arbitrary vacuum free space box
with dimension $\delta=\\{2,3\\}$ and
${\mathcal{B}}\subset{\mathbb{R}}^{\delta}$ denote the domain occupied by the
material solid, as demonstrated in Figure 2. Vacuum space $\Omega$ is
considered to be large enough such that the magnetic field induced by the
magnetization of the body ${\mathcal{B}}$ is decayed at its surface
$\partial{\mathcal{B}}\subset{\mathbb{R}}^{\delta}$. Also, in the following,
let a Dirichlet boundaries conditions of vacuum defined as
$\partial_{D}\Omega=\partial^{0}_{D}\Omega\cup\partial^{1}_{D}\Omega$ which
contains inner surface (it has an intersect with the solid body
${\mathcal{B}}$) and outer boundary surface, respectively, as shown in Figure
2b . We further assume Dirichlet boundary conditions on
$\partial_{D}{\mathcal{B}}$ and complete Neumann boundary conditions on
$\Gamma_{N}:=\partial_{N}{\mathcal{B}}\cup\mathcal{C}$, where
$\partial_{N}{\mathcal{B}}$ denotes the outer domain boundary (where the
traction imposed) and ${\mathcal{C}}\in\mathbb{R}^{\delta-1}$ is the crack
boundary, as illustrated in Figure 2c. Therefore, we have three different
domains required for the electromagnetically induced fracture to be used as
follows:

$\texttt{Vacuum region:}\;\Omega\subset{\mathbb{R}}^{\delta},\;\;\texttt{Solid
region:}\;\;{\mathcal{B}}\subset{\mathbb{R}}^{\delta}\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad\texttt{Damage
region:}\;\;{\mathcal{C}}\in\mathbb{R}^{\delta-1}.$

The response of the fracturing solid at material points
${\bm{\mathnormal{x}}}\in{\mathcal{B}}$ at time $t\in{\mathcal{T}}=[0,T]$ is
described by the displacement field
${\bm{\mathnormal{u}}}({\bm{\mathnormal{x}}},t)$ and the crack phase-field
$d({\bm{\mathnormal{x}}},t)$ as

${\bm{\mathnormal{u}}}:\left\\{\begin{array}[]{ll}{\mathcal{B}}\times{\mathcal{T}}\rightarrow\mathbb{R}^{\delta}\\\\[5.69054pt]
({\bm{\mathnormal{x}}},t)\mapsto{\bm{\mathnormal{u}}}({\bm{\mathnormal{x}}},t)\end{array}\right.\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad
d:\left\\{\begin{array}[]{ll}{\mathcal{B}}\times{\mathcal{T}}\rightarrow[0,1]\\\\[5.69054pt]
({\bm{\mathnormal{x}}},t)\mapsto
d({\bm{\mathnormal{x}}},t)\end{array}\right.\;.$ (1)

(a) (b) (c) Figure 2: Primary variable fields in magneto-elasticity, for a
solid body ${\mathcal{B}}\subset{\mathcal{R}}^{\delta}$ with dimension
$\delta\in[2,3]$. $(a)$ The displacement field ${\bm{\mathnormal{u}}}$ defined
on ${\mathcal{B}}$ and Neumann-type boundary condition for traction
$\bar{{\bm{\mathnormal{t}}}}={\bm{\tau}}\cdot{\bm{\mathnormal{n}}}$ on
$\partial_{N}{\mathcal{B}}$. $(b)$ The magnetic vector potential
${\bm{\mathnormal{A}}}$ in $\Omega$ that is continuous across the interface
$\partial{{\mathcal{B}}}$, i.e.
$\llbracket{\bm{\mathnormal{A}}}\rrbracket:={\bm{\mathnormal{A}}}^{+}-{\bm{\mathnormal{A}}}^{-}=\textbf{0}$
on $\partial{{\mathcal{B}}}$, and zero on the boundary $\partial{\Omega}$ free
space box (vacuum). $(c)$ The crack phase-field is determined by Dirichlet-
type boundary condition $d=1$ on $\mathcal{C}$ and Neumann-type boundary
condition $\nabla d\cdot{\bm{\mathnormal{n}}}=0$ on the full surface
$\partial{{\mathcal{B}}}$.

Intact and fully fractured states of the material are characterized by
$d({\bm{\mathnormal{x}}},t)=0$ and $d({\bm{\mathnormal{x}}},t)=1$,
respectively. In order to derive the variational formulation, the following
space is first defined. For an arbitrary $A\subset\mathbb{R}^{\delta}$, we set

$\displaystyle\mathrm{H}^{1}({\mathcal{B}},A):=\Big{\\{}v:{\mathcal{B}}\times{\mathcal{T}}\rightarrow
A\;\colon\;\;v\in\mathrm{H}^{1}(A)\Big{\\}}.$ (2)

We also denote the vector valued space
$\mathbf{H}^{1}({\mathcal{B}},A):=\left[\mathrm{H}^{1}({\mathcal{B}},A)\right]^{\delta}$
and define

${\mathcal{W}}^{{\bm{\mathnormal{u}}}}(\text{div};{\mathcal{B}}):=\Big{\\{}{\bm{\mathnormal{u}}}\in\mathbf{H}^{1}({\mathcal{B}},\mathbb{R}^{\delta})\;\
\colon\;\ {\bm{\mathnormal{u}}}=\overline{{\bm{\mathnormal{u}}}}\ \text{on}\
\partial_{D}{\mathcal{B}}\Big{\\}},$ (3)

and correspondingly for displacement test function, we have

${\mathcal{W}}_{0}^{{\bm{\mathnormal{u}}}}(\text{div};{\mathcal{B}}):=\Big{\\{}{\bm{\mathnormal{u}}}\in\mathbf{H}^{1}({\mathcal{B}},\mathbb{R}^{\delta})\\}.$
(4)

Concerning the crack phase-field, we set

${\mathcal{W}}^{d}:=\mathrm{H}^{1}({\mathcal{B}})\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad{\mathcal{W}}^{d}_{d^{n}}(\text{div};{\mathcal{B}}):=\Big{\\{}d\in\mathrm{H}^{1}({\mathcal{B}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}[0,1]})\;\colon\;\
d\geq d^{n}\Big{\\}},$ (5)

where $d^{n}$ is the damage value in a previous time instant. Note that
${\mathcal{W}}^{d}_{d^{n}}(\text{div};{\mathcal{B}})$ is a non-empty, closed,
and convex subset of ${\mathcal{W}}^{d}$ and introduces the evolutionary
character of the phase-field, incorporating an irreversibility condition in
incremental form. Correspondingly the phase-field test function reads

${\mathcal{W}}^{d}_{0}(\text{div};{\mathcal{B}}):=\Big{\\{}d\in\mathrm{H}^{1}({\mathcal{B}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}[0,1]})\;\colon\;\
d\geq 0\Big{\\}}.$ (6)

In order to formulate a wide variety of the failure mechanism of the electro-
and magneto-active materials into the variational equations, two different
sets of fields need to be introduced. To this end, the electric and magnetic
primary fields are introduced:

${\bm{\mathnormal{E}}}:\left\\{\begin{array}[]{ll}\Omega\times{\mathcal{T}}\rightarrow\mathbb{R}^{\delta}\\\
({\bm{\mathnormal{x}}},t)\mapsto{\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)\end{array}\right.\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad{\bm{\mathnormal{H}}}:\left\\{\begin{array}[]{ll}\Omega\times{\mathcal{T}}\rightarrow\mathbb{R}^{\delta}\\\
({\bm{\mathnormal{x}}},t)\mapsto{\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t)\end{array}\right.\;,$
(7)

such that, the intensity fields
$({\bm{\mathnormal{E}}},{\bm{\mathnormal{H}}})$ belong to the space:

${\mathcal{S}}^{{\bm{\mathnormal{v}}}}(\text{curl};\Omega):=\Big{\\{}{\bm{\mathnormal{v}}}\in\mathbf{L}^{2}(\Omega,\mathbb{R}^{\delta})\;\colon\;\nabla\times{\bm{\mathnormal{v}}}\in\mathbf{L}^{2}(\Omega,\mathbb{R}^{\delta})\Big{\\}},$
(8)

along with their test function space as

${\mathcal{S}}_{0}^{{\bm{\mathnormal{v}}}}(\text{curl};\Omega):=\Big{\\{}{\bm{\mathnormal{v}}}\in{\mathcal{S}}^{{\bm{\mathnormal{v}}}}(\text{curl};\Omega)\;\colon\;{\bm{\mathnormal{v}}}\times{\bm{\mathnormal{n}}}=0\in\mathbf{L}^{2}(\Omega,\mathbb{R}^{\delta})\Big{\\}}.$
(9)

The generalized fluxes conjugate to (7) for the electric and magnetic primary
fields read:

${\bm{\mathnormal{D}}}:\left\\{\begin{array}[]{ll}\Omega\times{\mathcal{T}}\rightarrow\mathbb{R}^{\delta}\\\
({\bm{\mathnormal{x}}},t)\mapsto{\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t)\end{array}\right.\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad{\bm{\mathnormal{B}}}:\left\\{\begin{array}[]{ll}\Omega\times{\mathcal{T}}\rightarrow\mathbb{R}^{\delta}\\\
({\bm{\mathnormal{x}}},t)\mapsto{\bm{\mathnormal{B}}}({\bm{\mathnormal{x}}},t)\end{array}\right.\;.$
(10)

Correspondingly, the flux fields
$({\bm{\mathnormal{D}}},{\bm{\mathnormal{B}}})$ are defined in the space:

${\mathcal{S}}^{{\bm{\mathnormal{v}}}}(\text{div};\Omega):=\Big{\\{}{\bm{\mathnormal{v}}}\in\mathbf{L}^{2}(\Omega,\mathbb{R}^{\delta})\;\colon\;\nabla.{\bm{\mathnormal{v}}}\in\mathbf{L}^{2}(\Omega,\mathbb{R}^{\delta})\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad{\bm{\mathnormal{v}}}=\overline{{\bm{\mathnormal{v}}}}\
\text{on}\ \partial^{0}_{D}\Omega\Big{\\}},$ (11)

and augmented with their test function space, as

${\mathcal{S}}_{0}^{{\bm{\mathnormal{v}}}}(\text{div};\Omega):=\Big{\\{}{\bm{\mathnormal{v}}}\in{\mathcal{S}}^{{\bm{\mathnormal{v}}}}(\text{div};\Omega)\;\colon\;{\bm{\mathnormal{v}}}.{\bm{\mathnormal{N}}}=0\;\;\text{on}\;\
\partial_{D}\Omega:=\partial^{0}_{D}\Omega\cup\partial^{1}_{D}\Omega\Big{\\}}.$
(12)

Here, ${\bm{\mathnormal{N}}}$ is the outward unit normal vector on the surface
$\partial\Omega$, as illustrated in Figure 2.

These electric and magnetic primary fields and their fluxes describe the
fundamental macroscopic fields governing all electromagnetic phenomena denoted
as Maxwell’s equations. These are described in detail in Section 2.2.2. To
complete this section, we need to introduce the vector potential formulation
of Maxwell’s equations denoted as
${\bm{\mathnormal{A}}}({\bm{\mathnormal{x}}},t)$. Since, the $div$ and $curl$
operator are surjective, following [35, 36] the de Rham complex reads

$\mathbb{R}\xrightarrow{id}{H^{1}(\Omega)}\xrightarrow{{}^{\text{grad}}}H(\text{curl},\Omega)\xrightarrow{{}^{\text{curl}}}{H(\text{div},\Omega)}\xrightarrow{{}^{\text{div}}}{L_{2}(\Omega)}\xrightarrow{{}^{\text{0}}}0.$
(13)

For a magnetic primary field
${\bm{\mathnormal{B}}}\in{\mathcal{S}}^{{\bm{\mathnormal{v}}}}(\text{div};\Omega)$
it holds

$\nabla.{\bm{\mathnormal{B}}}=0\quad\xrightarrow{\text{de
Rham}}\quad\exists\;{\bm{\mathnormal{A}}}\in{\mathcal{S}}^{{\bm{\mathnormal{v}}}}(\text{curl};\Omega)\;\;\colon\nabla\times{\bm{\mathnormal{A}}}=:{\bm{\mathnormal{B}}},$
(14)

in which, ${\bm{\mathnormal{A}}}$ denotes the magnetic vector potential. To
summarize, we have following sets:

$\displaystyle({\bm{\mathnormal{u}}},d)\in\Big{(}{\mathcal{W}}^{{\bm{\mathnormal{u}}}}(\text{div};{\mathcal{B}}),$
$\displaystyle{\mathcal{W}}^{d}(\text{div};{\mathcal{B}})\Big{)},\qquad({\bm{\mathnormal{E}}},{\bm{\mathnormal{H}}},{\bm{\mathnormal{A}}})\in{\mathcal{S}}^{{\bm{\mathnormal{v}}}}(\text{curl};\Omega)$
(15)
$\displaystyle({\bm{\mathnormal{D}}},{\bm{\mathnormal{B}}})\in{\mathcal{S}}^{{\bm{\mathnormal{v}}}}(\text{div};\Omega),$

and correspondingly it follows for their test functions:

$\displaystyle(\delta{\bm{\mathnormal{u}}},\delta
d)\in\Big{(}{\mathcal{W}}_{0}^{{\bm{\mathnormal{u}}}}(\text{div};{\mathcal{B}}),$
$\displaystyle{\mathcal{W}}_{0}^{d}(\text{div};{\mathcal{B}})\Big{)},\qquad(\delta{\bm{\mathnormal{E}}},\delta{\bm{\mathnormal{H}}},\delta{\bm{\mathnormal{A}}})\in{\mathcal{S}}_{0}^{{\bm{\mathnormal{v}}}}(\text{curl};\Omega),$
(16)
$\displaystyle(\delta{\bm{\mathnormal{D}}},\delta{\bm{\mathnormal{B}}})\in{\mathcal{S}}_{0}^{{\bm{\mathnormal{v}}}}(\text{div};\Omega).$

###  2.2 .Governing equations of the failure mechanism for magnetostrictive
effects

In this section, we present a theoretical development associated with the
failure mechanism of the ferromagnetic material based on the variational
phase-field model by taking into account their fully coupled magneto-
mechanical characteristics. This includes mechanical, electromagnetic, and
fracture contributions to the model. Additionally, the transition rule from
undamaged to fully damaged states is presented. This section is used as a
departure point for Section 3 in which to formulate the variational framework
for magnetostrictive-induced cracking is formulated.

####  2.2.1 .Elastic contribution.

The standard elastic energy density, so-called the effective strain energy
density [37, 38] is expressed in terms of the scalar-valued function
$\psi_{e}$. For isotropic materials, $\psi_{e}$ is defined in terms of the
strain tensor ${\bm{\varepsilon}}$. In our formulation, in order to preclude
fracture in compression, a decomposition of the effective strain energy
density into damageable and undamageable parts are employed. Thus, we perform
additive decomposition of the strain tensor into volume-changing (volumetric)
and volume-preserving (deviatoric) counterparts:

$\bm{\varepsilon}(\bm{u})=\bm{\varepsilon}^{vol}(\bm{u})+\bm{\varepsilon}^{dev}(\bm{u})\quad\mathchoice{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}\quad\bm{\varepsilon}(\bm{u}):=\frac{1}{2}(\nabla{{\bm{\mathnormal{u}}}}+\nabla^{T}{{\bm{\mathnormal{u}}}}),$
(17)

where the volumetric strain is denoted as
$\bm{\varepsilon}^{vol}(\bm{u}):=\frac{1}{3}(\bm{\varepsilon}(\bm{u}):\text{\bm{\mathnormal{I}}})\text{\bm{\mathnormal{I}}}$
and the deviatoric strain is denoted as
$\bm{\varepsilon}^{dev}(\bm{u}):=\mathbb{P}:\bm{\varepsilon}$. The fourth-
order projection tensor
$\mathbb{P}:=\mathbb{I}-\frac{1}{3}\text{\text{\bm{\mathnormal{I}}}}\otimes\text{\bm{\mathnormal{I}}}$
is introduced to map the full strain tensor onto its deviatoric component.
Therein,
$\mathbb{I}_{ijkl}:=\frac{1}{2}\big{(}\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\big{)}$
is the fourth-order symmetric identity tensor.

Next we construct the mechanical BVP. The solid geometry ${\mathcal{B}}$ is
loaded by prescribed deformations $\overline{{\bm{\mathnormal{u}}}}$ and an
external traction vector $\overline{\bm{t}}$ on the boundary, defined by time-
dependent Dirichlet conditions and Neumann conditions as

${\bm{\mathnormal{u}}}=\overline{{\bm{\mathnormal{u}}}}({\bm{\mathnormal{x}}},t)\
\textrm{on}\
\partial^{{\bm{\mathnormal{u}}}}_{D}{\mathcal{B}}\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad{\bm{\tau}}\cdot{\bm{\mathnormal{n}}}=\overline{{\bm{\mathnormal{t}}}}({\bm{\mathnormal{x}}},t)\
\textrm{on}\ \partial^{{\bm{\mathnormal{u}}}}_{N}{\mathcal{B}},$ (18)

where ${\bm{\mathnormal{n}}}$ is the outward unit normal vector on the surface
$\partial_{N}{\mathcal{B}}$. The stress tensor ${\bm{\tau}}$ is the
thermodynamic dual to the strain ${\bm{\varepsilon}}$.

The global mechanical form of the equilibrium equation for the solid body can
be represented by a second-order PDE for the multi-field system as

$\framebox{$\displaystyle\mbox{Div}\,{\bm{\tau}}({\bm{\mathnormal{u}}},{\bm{\mathnormal{B}}},d)+\overline{}{\bm{\mathnormal{b}}}=\textbf{0}\,,$}$
(19)

which is valid for a quasi-static response and where denote
$\overline{}{\bm{\mathnormal{b}}}$ is the prescribed body force.

####  2.2.2 .Electromagnetic contribution.

Maxwell’s equations can principally express the macroscopic electromagnetic
phenomena in a ferromagnetic material. These equations are related to the
calculations of electromagnetic fields, including the electric and the
magnetic fields with eddy currents (their corresponding fluxes)

$\begin{split}\text{Maxwell-Faraday
law}:&&\nabla\times{\bm{\mathnormal{E}}}&=-\frac{\partial{{\bm{\mathnormal{B}}}}}{\partial{t}}\hskip
42.67912pt&\text{in}\;\;\Omega\times(0,\text{T}),\\\\[7.11317pt]
\text{Maxwell-Amp\\`{e}re
law}:&&\nabla\times{\bm{\mathnormal{H}}}&=\frac{\partial{{\bm{\mathnormal{D}}}}}{\partial{t}}+{\bm{\mathnormal{J}}}\hskip
42.67912pt&\text{in}\;\;\Omega\times(0,\text{T}),\\\\[7.11317pt]
\text{Maxwell-Gauss law}:&&\nabla\cdot{\bm{\mathnormal{D}}}&=\rho_{0}(d)\hskip
42.67912pt&\text{in}\;\;\Omega\times(0,\text{T}),\\\\[7.11317pt] \text{No
magnetic monopole law}:&&\nabla\cdot{\bm{\mathnormal{B}}}&=0\hskip
42.67912pt&\text{in}\;\;\Omega\times(0,\text{T}),\end{split}$ (20)

where ${\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)$ and
${\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t)$ are the electric and magnetic
fields, ${\bm{\mathnormal{D}}}({\bm{\mathnormal{x}}},t)$ and
${\bm{\mathnormal{B}}}({\bm{\mathnormal{x}}},t)$ are the corresponding
electric and magnetic flux densities, respectively, along with
${\bm{\mathnormal{J}}}({\bm{\mathnormal{x}}},t)$ that is the electric current
density, and $\rho$ is the electric charge density. In addition to these
equations, there are constitutive equations that describe the macroscopic
characteristics of a electromagnetic materials. These equations construct
relations between the magnetic and electric flux densities and the magnetic
and the electric fields

${\bm{\mathnormal{D}}}=\epsilon_{0}(d){\bm{\mathnormal{E}}},\quad\quad{\bm{\mathnormal{B}}}={\mu_{0}}(d){\bm{\mathnormal{H}}}\quad\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad\quad{\bm{\mathnormal{J}}}=\sigma_{0}(d){\bm{\mathnormal{E}}}+{\bm{\mathnormal{J}}}_{s}\;.$
(21)

######  Remark 2.1.

We note that polarization is to be considered in the case of electromagnetic
wave propagation that contains the electric and magnetic fields whose
oscillations are perpendicular to the wave direction. The polarization is a
feature of an electromagnetic wave that determines the geometrical orientation
of the oscillations [39]. When an electromagnetic wave propagates in a
homogeneous isotropic medium, the oscillations of electric and magnetic fields
are perpendicular to each other and perpendicular to wave propagation’s
direction. The polarization expresses the direction of the electric field.
Depending on the oscillation direction of the fields, there are two types of
polarization, linear and circular the latter can be specified as the left
circular and right circular polarization [40]. The polarization effect is
regarded in Maxwell’s equation by modifying (21) in terms of the electric
magnetic field ${\bm{\mathnormal{D}}}$ as [12]

${\bm{\mathnormal{D}}}=\epsilon_{0}(d){\bm{\mathnormal{E}}}+{\bm{\mathnormal{P}}},$
(22)

where ${\bm{\mathnormal{P}}}$ represents the polarization of the
electromagnetic wave propagation. In our work we further assume
${\bm{\mathnormal{P}}}$ has negligible effect, thus
${\bm{\mathnormal{P}}}=\textbf{0}$. Therefore, this effect is not anymore take
into account in our study.

Here, ${\bm{\mathnormal{J}}}_{s}$ denotes the specific current density. Also,
material constants ${\epsilon_{0}}$, ${\mu_{0}}$, and $\sigma_{0}$ represent
the permittivity, the permeability, and the conductivity of the medium,
respectively. We further note that ${\epsilon_{0}}$, ${\mu_{0}}$, and
$\sigma_{0}$ are piecewise smooth, real, bounded, and positive and may vary in
different Cartesian coordinates. We need to mention that dependency of these
variables to the crack phase-field is further elaborated in Section 2.2.4.
Thus, we have for every ${\bm{\mathnormal{x}}}\in\Omega$:

$\displaystyle 0<{\epsilon_{0,min}}<{\epsilon_{0}}<{\epsilon_{0,max}}<\infty,$
(23) $\displaystyle 0<{\mu_{0,min}}<{\mu_{0}}<{\mu_{0,max}}<\infty,$
$\displaystyle 0<{\sigma_{0,max}}<\infty.$

Additional to (20) and (21), the Dirichlet boundary condition are for the
electric field and magnetic vector potential:

${\bm{\mathnormal{E}}}=\bar{{\bm{\mathnormal{E}}}}({\bm{\mathnormal{x}}},t)\
\textrm{on}\
\partial^{1}_{D}\Omega~{}~{}\text{and}~{}~{}{\bm{\mathnormal{A}}}=\bar{{\bm{\mathnormal{A}}}}({\bm{\mathnormal{x}}},t)\
\textrm{on}\ \partial^{1}_{D}\Omega.$ (24)

On the boundary of the free space, outer Dirichlet and Neumann boundary
conditions have to be satisfied as:

$\bar{{\bm{\mathnormal{A}}}}=\textbf{0}\quad\text{on}\;\;\partial^{0}_{D}{\Omega}\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad\nabla\bar{{\bm{\mathnormal{A}}}}\cdot{\bm{\mathnormal{N}}}=\textbf{0}\quad\text{on}\;\;\partial^{0}_{D}{\Omega},$
(25)

respectively, such that:

$\displaystyle\lim\limits_{{\bm{\mathnormal{x}}}\to\infty}\bar{{\bm{\mathnormal{A}}}}({\bm{\mathnormal{x}}},t)\rightarrow\textbf{0}.$

In the case of the magnetized material, the constitutive relation between
magnetic flux ${\bm{\mathnormal{H}}}$ and magnetic field
${\bm{\mathnormal{B}}}$ needs to be further modified. This adjustment will
help to accommodate our proposed formulation of a magneto-mechanical problem
of fracturing solids. To do so, we describe three common types of macroscopic
magnetization responses as follows:

(a) (b) (c) Figure 3: Partial alignment of the magnetic moments in the
presence of the external magnetic field for classified based on relative
magnetic permeability of materials, (a) diamagnetic with $-1\leq 1+\chi<0$,
(b) for paramagnetic with $0\leq 1+\chi<\text{eps}$, and (c) for ferromagnetic
$1\ll\chi$ in parallel orientation of magnetic moment in the presence of the
external magnetic field.

  * •

Diamagnetic. These types of materials are those which become weakly magnetized
when they are affected by an external magnetic field. The tendency of
magnetization for paramagnetic materials is to move in the direction of a
strong magnetic field towards weak parts of the external magnetic field [41].
So, diamagnetic materials are magnetized in the opposite direction of the
external magnetic field. Moreover, these material present no magnetic
hysteresis.

  * •

Paramagnetic. Similar to the diamagnetic response, this type of material
becomes weakly magnetized. Unlike diamagnetic materials, the tendency of
magnetization for paramagnetic materials is to move in the direction of weak
magnetic fields to strong ones [41]. So, paramagnetic materials are magnetized
in the direction of the external magnetic field. Here, small magnetic
hysteresis could be observed .

  * •

Ferromagnetic. This class of magnetic materials is described through a
permanent magnetization effect and mainly has a profound response to magnetic
fields. These materials are characterized by a highly nonlinear response in
which magnetic hysteresis typically occurs [42].

To formulate a different class of magnetic material, the constitutive relation
between ${\bm{\mathnormal{B}}}$ and ${\bm{\mathnormal{H}}}$ is described
through:

${\bm{\mathnormal{B}}}=\bm{\mu}\cdot{\bm{\mathnormal{H}}}\quad\mathchoice{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}\quad\bm{\mu}=\mu_{0}(\text{\bm{\mathnormal{I}}}+\bm{\chi}),$
(26)

in which magnetic susceptibility $\bm{\chi}$ (as a dimensionless quantity),
measures the degree of magnetization response of material if there is an
external magnetic field leading to:

${\bm{\mathnormal{M}}}=\bm{\chi}\cdot{\bm{\mathnormal{H}}}\quad\mathchoice{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}\quad\bm{\chi}=\frac{\partial{\bm{\mathnormal{M}}}}{\partial{\bm{\mathnormal{H}}}}.$
(27)

The term $\text{\bm{\mathnormal{I}}}+\bm{\chi}$ in (26) typically refers to
the relative magnetic permeability of materials. Thus, the coupling of (26)
and (27) results in:

${\bm{\mathnormal{B}}}=\mu_{0}(\text{\bm{\mathnormal{I}}}+\bm{\chi}):{\bm{\mathnormal{H}}}={\mu_{0}}({\bm{\mathnormal{H}}}+{\bm{\mathnormal{M}}}).$
(28)

On the basis of the magnitude of magnetic susceptibility $||\bm{\chi}||_{2}$
in (28), magnetized material is classified into the above mentioned:

$\displaystyle\texttt{Diamagnetic: }-1\leq||\bm{\chi}||_{2}<0,\quad$ (29)
$\displaystyle\texttt{Paramagnetic: }0<||\bm{\chi}||_{2}<\text{eps}\approx
0,\quad$ $\displaystyle\texttt{Ferromagnetic: }||\bm{\chi}||_{2}\gg 1,$

see Figure 3. In our case, we use the magnetic susceptibility as an isotropic
tensor for the ferromagnetic response:
$\bm{\chi}=\chi_{0}\text{\bm{\mathnormal{I}}}$ with $\chi_{0}\gg 1$, which is
valid for a variety of amorphous solids or materials with a uniform crystal
structure [43]. Since ${\bm{\mathnormal{B}}}$ is solenoidal (divergence-free)
according to Maxwell’s equations, we know through (14) that
${\bm{\mathnormal{B}}}$ must be the curl of the magnetic vector potential
field ${\bm{\mathnormal{A}}}$, so
${\bm{\mathnormal{B}}}=\nabla\times{\bm{\mathnormal{A}}}$. In a common way,
the electromagnetic fields, ${\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)$,
${\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t)$,
${\bm{\mathnormal{D}}}({\bm{\mathnormal{x}}},t)$ and
${\bm{\mathnormal{B}}}({\bm{\mathnormal{x}}},t)$ are presumed to be time-
harmonic. It means that the fields can harmonically oscillate with a single
frequency $\omega$. In such cases, they can be written as

$\displaystyle{\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)$
$\displaystyle=\exp(-i\omega
t){\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},\omega),$ (30)
$\displaystyle{\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t)$
$\displaystyle=\exp(-i\omega
t){\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},\omega),$
$\displaystyle{\bm{\mathnormal{J}}}_{s}({\bm{\mathnormal{x}}},t)$
$\displaystyle=\exp(-i\omega
t){\bm{\mathnormal{J}}}_{s}({\bm{\mathnormal{x}}},\omega).$

The first and the second-order derivatives are

$\displaystyle\frac{\partial{\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)}{\partial
t}=-i\omega{\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t),$
$\displaystyle\frac{\partial^{2}{\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)}{\partial
t^{2}}=-\omega^{2}{\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t),$ (31)
$\displaystyle\frac{\partial{\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t)}{\partial
t}=-i\omega{\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t),$
$\displaystyle\frac{\partial^{2}{\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t)}{\partial
t^{2}}=-\omega^{2}{\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t),$
$\displaystyle\frac{\partial{\bm{\mathnormal{J}}}_{s}({\bm{\mathnormal{x}}},t)}{\partial
t}=-i\omega{\bm{\mathnormal{J}}}_{s}({\bm{\mathnormal{x}}},t),$
$\displaystyle\frac{\partial^{2}{\bm{\mathnormal{J}}}_{s}({\bm{\mathnormal{x}}},t)}{\partial
t^{2}}=-\omega^{2}{\bm{\mathnormal{J}}}_{s}({\bm{\mathnormal{x}}},t).$

Therefore, the time-harmonic Maxwell’s equations are obtained by substituting
(30) and (21) into (20)

$\left\\{\,\vbox{\halign{\strut@$\m@th\displaystyle#$\hfil&\quad\strut@{#}\hfil\cr\hbox{\vrule
height=0.0pt,depth=0.0pt,width=0.0pt}$\displaystyle\nabla\times{\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)=i\omega\mu_{0}{\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t)&\hbox{\vrule
height=0.0pt,depth=0.0pt,width=0.0pt}{\quad\text{in}\;
$\Omega\times(0,\text{T})$\\\\\\\
\nabla\times{\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t) =-i
\omega\left(\epsilon_0 + i
\frac{\sigma_{0}}{\omega}\right){\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)
+ {\bm{\mathnormal{J}}}_s({\bm{\mathnormal{x}}},t) &\quad\text{in}\;
$\Omega\times(0,\text{T})$\crcr}}\right..}$ (32)

We now take a curl function, i.e., $\nabla\times$, from the first equation of
(32), which results in:

$\nabla\times{\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t)=\nabla\times\left(\frac{1}{i\omega\mu_{0}}\nabla\times{\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)\right).$
(33)

Next, we use the curl function from the second expression of (32) and with
(33), we obtain the equation for electric field
${\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)$ as

$\nabla\times\left(\frac{1}{\mu_{0}}\nabla\times{\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)\right)-\omega^{2}\left(\epsilon_{0}+i\frac{\sigma_{0}}{\omega}\right){\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)=i\omega{\bm{\mathnormal{J}}}_{s}({\bm{\mathnormal{x}}},t).$
(34)

In a similar manner, we can derive the equation for the magnetic field
${\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t)$ by taking the curl from the
second expression of (32) as

$\nabla\times(\nabla\times{\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t))=-i\omega\left(\epsilon_{0}+i\frac{\sigma_{0}}{\omega}\right)\nabla\times{\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)+\nabla\times{\bm{\mathnormal{J}}}_{s}({\bm{\mathnormal{x}}},t),$
(35)

where
$\nabla\times{\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)=i\omega\mu_{0}{\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t)$.
Therefore, we obtain

$\nabla\times(\nabla\times{\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t))-\omega^{2}\mu_{0}\left(\epsilon_{0}+i\frac{\sigma_{0}}{\omega}\right)\
{\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t)=\nabla\times{\bm{\mathnormal{J}}}_{s}({\bm{\mathnormal{x}}},t).$
(36)

The time-domain vector wave equations in terms of the electric field
${\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)$ and the magnetic field
${\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},t)$ are further simplified by
substituting (31) into (34) and (36). This yields:

$\left\\{\,\vbox{\halign{\strut@$\m@th\displaystyle#$\hfil&\quad\strut@{#}\hfil\cr\hbox{\vrule
height=0.0pt,depth=0.0pt,width=0.0pt}$\displaystyle\epsilon_{0}\frac{\partial^{2}{\bm{\mathnormal{E}}}}{\partial
t^{2}}+\sigma_{0}\frac{\partial{\bm{\mathnormal{E}}}}{\partial
t}+\nabla\times\left(\frac{1}{\mu_{0}}\nabla\times{\bm{\mathnormal{E}}}\right)=-\frac{\partial{\bm{\mathnormal{J}}}_{s}}{\partial
t}&\hbox{\vrule height=0.0pt,depth=0.0pt,width=0.0pt}{\quad\text{in}
\quad$\Omega\times(0,\text{T})$ \\\\\\\ \epsilon_0 \frac{\partial^{2}
{\bm{\mathnormal{H}}}}{\partial t^{2}} + \sigma_0
\frac{\partial{\bm{\mathnormal{H}}}}{\partial t} +
\nabla\times\left(\frac{1}{\mu_{0}} \nabla\times{\bm{\mathnormal{H}}}\right) =
\frac{1}{\mu_{0}}\nabla\times{\bm{\mathnormal{J}}}_s &\quad\text{in}
\quad$\Omega\times(0,\text{T})$\\\\\crcr}}\right..}$ (37)

In order to demonstrate the modeling capacity of (37), a representative
simulation is presented in Figure 4. In this boundary-value problem, the
evolution of the electric field ${\bm{\mathnormal{E}}}$ is represented in the
$x$ and $y$ directions under the variation of the specific current
${\bm{\mathnormal{J}}}_{s}$. In this example, the specific current
${\bm{\mathnormal{J}}}_{s}$ is defined as a bi-linear function of the nodal
coordinates and time. As a boundary condition, the electric field is set to
${\bm{\mathnormal{E}}}=\textbf{0}$ at the boundaries of the vacuum. The
results show that the variation of the electric field is symmetric along the
$x$ and $y$ directions as the electric specific current is uniformly applied
on the silicon beam. In fact, this is expected due to symmetric response of
electric field imposed to domain.

(a) (b) (c)  Figure 4: The geometry of BVP, and the the representation of the
$x$\- and $y$\- component of the electric field ${\bm{\mathnormal{E}}}$.

Next, we reduce (37) toward Maxwell’s equations in terms of vector potential.
The Maxwell-Faraday law in (20), based on magnetic vector potential
${\bm{\mathnormal{A}}}({\bm{\mathnormal{x}}},t)$ given in (14), is then
related to electric field through:

$\nabla\times{\bm{\mathnormal{E}}}=-\frac{\partial({\nabla\times{\bm{\mathnormal{A}}}})}{\partial{t}}\quad\rightarrow\quad\exists\;\Phi({\bm{\mathnormal{x}}},t)\;:\;\;-\nabla\Phi={\bm{\mathnormal{E}}}+\frac{\partial{\bm{\mathnormal{A}}}}{\partial
t}.$ (38)

Here, $\Phi({\bm{\mathnormal{x}}},t)$ is so-called the electric potential, see
for a detailed discussion [44, 45]. From the mathematical point of view, the
homogeneous Maxwell equations for Maxwell-Faraday law (20)1, and Maxwell-
Ampère law in (20)2 results in the existence set of
$({\bm{\mathnormal{A}}},\Phi)$, and thus implies (38). We refer interested
readers to [46] for more details. Considering Maxwell-Ampère law in (20), and
constitutive equations (21), yields:

$\frac{1}{\mu_{0}}\nabla\times{\bm{\mathnormal{B}}}=\epsilon_{0}(d)\frac{\partial{\bm{\mathnormal{E}}}}{\partial
t}+{\bm{\mathnormal{J}}}$ (39)

that is

$\nabla\times\left(\frac{1}{\mu_{0}}\nabla\times{\bm{\mathnormal{A}}}\right)=-\epsilon_{0}(d)\left[\frac{\partial\nabla\Phi}{\partial
t}+\frac{\partial^{2}{\bm{\mathnormal{A}}}}{\partial
t^{2}}\right]-\sigma_{0}(d)\left[\nabla\Phi+\frac{\partial{\bm{\mathnormal{A}}}}{\partial
t}\right]+{\bm{\mathnormal{J}}}_{s}$ (40)

where (38) and (21) are used. Next, we assume that for a time dependent
magnetic vector potential ${\bm{\mathnormal{A}}}({\bm{\mathnormal{x}}},t)$,
and the electric potential $\Phi({\bm{\mathnormal{x}}},t)$, we have:

$\mbox{Div}{\bm{\mathnormal{A}}}+\mu_{0}\epsilon_{0}\frac{\partial\Phi({\bm{\mathnormal{x}}},t)}{\partial
t}=0,$ (41)

which is the so-called Lorenz gauge condition [47, 45]. Thus, (40) together
with (41) reads:

$\framebox{$\displaystyle\frac{1}{\mu_{0}}\nabla^{2}{{\bm{\mathnormal{A}}}}=\epsilon_{0}(d)\left[\frac{\partial^{2}{\bm{\mathnormal{A}}}}{\partial
t^{2}}\right]+\sigma_{0}(d)\left[\nabla\Phi+\frac{\partial{\bm{\mathnormal{A}}}}{\partial
t}\right]-{\bm{\mathnormal{J}}}_{s}$}$ (42)

Here, the identity
$\nabla\times\nabla\times{\bm{\mathnormal{A}}}=\nabla(\mbox{Div}{\bm{\mathnormal{A}}})-\nabla^{2}{{\bm{\mathnormal{A}}}}$
is used.

We note that, the magnetic potential $A(x,t)$ as a solution of (42) is related
to crack parameter $d$ by PDE coefficients. So, the permittivity, the
permeability, and the conductivity of the medium
$\big{(}\epsilon_{0}(d),\mu_{0}(d),\sigma_{0}(d)\big{)}$ are function of $d$,
see Section 2.2.4.

The setting presented in [48] constitutes the starting point to formulate the
mechanical response within the electromagnetic analysis where the electric
scalar potential $\nabla{\Phi}$ in (42) is a given quantity. Following [48,
49], the magneto-static state $\nabla{\Phi}$ has the direction orthogonal to
the $x$-$y$ plane, i.e., $\nabla{\Phi}(t)=[0,0,\nabla\Phi_{0}(t)]$, where
$\Phi_{0}(t)$ is the pulse power supply parameter which is approximated as:

$\nabla\Phi_{0}(t)=a+b\,e^{-c\,t}.$ (43)

It is determined as a function of the constants $a$, $b$ and $c$ [48, 49] 222
Clearly, the representation on electric scalar potential given in (43) can be
derived through Maxwell’s equations, as well, and not as a given quantity. By
means of Maxwell-Gauss law in (20)3, and constitutive equation (21)1, we have:
$-\nabla^{2}{\Phi}-\nabla.\left(\frac{\partial{\bm{\mathnormal{A}}}}{\partial
t}\right)=-\frac{\rho_{0}(d)}{\epsilon_{0}(d)}$ (44) where (38) is used. Then,
by imposing Lorenz gauge condition in (44), results in the wave-like equation
by:
$c^{2}_{\Phi}\frac{\partial^{2}\Phi}{\partial^{2}t}=\nabla^{2}\Phi+f_{\Phi}\quad\text{in}\quad\Omega\times(0,\text{T}),$
(45) with
$f_{\Phi}=-\frac{\rho_{0}(d)}{\epsilon_{0}(d)}\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad
c^{2}_{\Phi}=\mu_{0}(d)\epsilon_{0}(d)$ (46) The the electric potential PDE
given (45) is a linear second-order differential equation with the variable
$\Phi\in H^{1}(\Omega,\mathbb{R})$ represents the wave-like propagates in the
free space with the constant speed of $c_{\Phi}$, and $f_{\Phi}$ states the
source of the wave propagation. Thus, set of equations (42) and (45) are
replaced by (37) to represent the Maxwell’s equations in terms of vector
potential.. In all our representative numerical examples in Section 4, these
values are set as $a=-1\times 10^{-3}$ V/m, $a=-1\times 10^{-3}$ V/m, and
$c=1\times 10^{-2}$ s-1.

####  2.2.3 .Fracture contribution.

In the smeared fracture framework, a sharp crack interface denoted by
$\mathcal{C}$ for satisfying the continuity of the crack topology is further
regularized, which is denoted as $\mathcal{C}_{l}$ as outlined [28]. To govern
a regularized fracture surface $\mathcal{C}_{l}$, it is required to
incorporate a continuous field variable – the so-called order parameter –
denoted by $d$, which differentiates between multiple physical phases within a
given system through a smooth transition. In the context of fracture, such an
order parameter (termed the crack phase-field) describes the smooth transition
between the fully broken ($d=1$) and intact material phases ($d=0$), thereby
approximating the sharp crack discontinuity. This geometrical perspective is
in agreement with the framework of [50], which was conceived as a
$\Gamma$-convergence regularization of the variational approach to Griffith
fracture [51]. A variety of research studies have recently extended the phase-
field approach to fracture toward the cohesive-frictional materials [52]
including thermal effects [53, 54, 55], ductile failure [56, 57, 58],
hydraulic fracture [59, 60, 61], stochastcic analysis [62, 63, 64, 65],
degradation of the fracture toughness [57], topology optimization [66], and
multi-scale approach [67, 38, 68, 69], and electro-mechanical approach [70,
71, 72] among others. In this manuscript, for the case of isotropic materials,
the regularized functional is given by

${\mathcal{C}}_{l}(d)=\int_{{\mathcal{B}}}\gamma(d,\nabla d)\,\text{d}{\cal
B},$ (47)

with possitivness for crack disspation as:

$\frac{d}{dt}{\mathcal{C}}_{l}(d)\geq 0\;.$ (48)

In line with standard phase-field models, a general surface density function
for the isotropic part $\gamma(d,\nabla d)$ is defined as

$\gamma(d,\nabla
d):=\frac{1}{c_{d}}\,\bigg{(}\frac{f(d)}{l_{d}}+\frac{l_{d}}{2}\nabla
d\cdot\nabla
d\bigg{)}\quad\mathchoice{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}\quad
c_{d}:=4\int_{0}^{1}\sqrt{\omega(b)}\,db,$ (49)

where $\omega(d)$ is a monotonic and continuous _local fracture energy
function_ such that $\omega(0)=0$ and $\omega(1)=1$. A variety of suitable
choices for $\omega(d)$ are available in the literature [73, 74, 75]. Here,
the widely adopted linear and quadratic formulations are considered, which
yield, models with and without an elastic stage, respectively. Specifically,
we define

$f(d):=\left\\{\,\vbox{\halign{\strut@$\m@th\displaystyle#$\hfil&\quad\strut@{#}\hfil\cr\hbox{\vrule
height=0.0pt,depth=0.0pt,width=0.0pt}$\displaystyle d\phantom{{}^{2}}\implies
c_{d}=8/3\quad&\hbox{\vrule height=0.0pt,depth=0.0pt,width=0.0pt}{\text{model
with an elastic stage}, \\\ d^2 \implies c_d=2 \quad&\text{model without an
elastic stage}.\crcr}}\right.}$ (50)

Following [76], the local evolution of the crack phase-field equation in the
domain ${\mathcal{B}}$ is:

$\framebox{$\displaystyle[\,d-l_{d}^{2}\Delta
d\,]+\eta_{d}\dot{d}+(1-\kappa)(d-1){{\mathcal{H}}}=0\,,$}$ (51)

augmented with its homogeneous Neumann boundary conditions, i.e., $\nabla
d\cdot{\bm{\mathnormal{n}}}=0$ on $\partial{\mathcal{B}}$. Following [37, 38],
the small residual scalar $0<\kappa\ll 1$ is introduced to prevent numerical
instabilities, which states the third equation in the coupled system.
Additionally, the damage viscosity material parameter denoted by $\eta_{d}\geq
0$ is used to characterize the viscosity term of the crack propagation. The
maximum absolute value for the crack driving state function denoted by
${\mathcal{D}}$ is defined by the crack driving force ${\mathcal{H}}$, which
reads

${\mathcal{H}}=\max_{s\in[0,t]}{\mathcal{D}}({\bm{\mathnormal{x}}},s)\geq 0,$
(52)

that accounts for the irreversibility of the crack phase-field evolution by
filtering out a maximum value of ${\mathcal{D}}$. We define the crack driving
state function in Section 3.

In summary, following formulation has to be solved for three-filed
$({\bm{\mathnormal{u}}},d,{\bm{\mathnormal{A}}})$ multi-physics problem.

######  Formulation 2.1 (Strong form of the Euler-Lagrange equations).

Let $\epsilon_{0}$, $\sigma_{0}$, $\mu_{0}$, $G_{c}$, $K$, and $\mu$ be given
with the initial conditions $\bm{u}^{0}=\bm{u}(\bm{x},0)$,
$d^{0}=d(\bm{x},0)$, and
${\bm{\mathnormal{A}}}^{0}={\bm{\mathnormal{A}}}(\bm{x},0)$. For the loading
increments $n=1,2,\ldots,N$, we solve a displacement equation where we seek
$\bm{u}:=\bm{u}^{n}:{\mathcal{B}}\rightarrow\mathbb{R}^{\delta}$ such that

$\displaystyle-\mathop{\operator@font
div}({\bm{\tau}})={\bm{\mathit{0}}}\quad\quad\quad\quad$ $\displaystyle
in\;{\mathcal{B}},$
$\displaystyle{\bm{u}}=\bar{}{\bm{\mathnormal{u}}}\quad\quad\quad\quad$
$\displaystyle on\;\partial_{D}{\mathcal{B}},$
$\displaystyle{\bm{\tau}}\cdot{\bm{\mathnormal{n}}}=\bar{{\bm{\mathnormal{t}}}}\quad\quad\quad\quad$
$\displaystyle on\;\partial_{N}{\mathcal{B}},$

in terms of the stress tensor ${\bm{\tau}}$ defined in (79) and the given
displacement field $\bar{}{\bm{\mathnormal{u}}}$. The phase-field system
consists of four parts: the PDE, the inequality constraint and a compatibility
condition along with the Neumann-type boundary conditions. Find
$d:=d^{n}:{\mathcal{B}}\rightarrow[0,1]$ such that

$\displaystyle(1-\kappa)(1-d){{\mathcal{H}}(I_{1},I_{2},t)}-[d-l_{d}^{2}\Delta
d\,]=\eta_{d}\dot{d}\qquad$ $\displaystyle in\;{\mathcal{B}},$
$\displaystyle\dot{d}\geqslant 0\qquad$ $\displaystyle in\;{\mathcal{B}},$

in terms of the crack driving force given (83), along with $\nabla
d\cdot{\bm{\mathnormal{n}}}=0\;\text{on}\;\partial\Omega$. along with a
second-order hyperbolic problem for the following magnetic vector potential.
Find
${\bm{\mathnormal{A}}}:={\bm{\mathnormal{A}}}^{n}:\Omega\rightarrow\mathbb{R}^{\delta}$
such that

$\displaystyle\quad\quad\quad\;\frac{1}{\mu_{0}}\nabla^{2}{{\bm{\mathnormal{A}}}}=\epsilon_{0}(d)\left[\frac{\partial^{2}{\bm{\mathnormal{A}}}}{\partial
t^{2}}\right]+\sigma_{0}(d)\left[\nabla\Phi+\frac{\partial{\bm{\mathnormal{A}}}}{\partial
t}\right]-{\bm{\mathnormal{J}}}_{s}\qquad\text{in}\;\Omega.$

where the Lorenz gauge condition is considered.

####  2.2.4 .Transition structure from undamaged to fully damaged states.

In this section, we provide a transition rule to connect the transition
between the intact (solid) region and the fracture domain. In fact, the
coupling of electromagnetic response to the crack phase field is achieved by
introducing constitutive functions, which are characterized by degraded
related material constants. Following [77], two state functions need to be
defined. The first set of the formulation is denoted as
$\texttt{P}^{s}({\bm{\mathnormal{x}}})\in{\mathcal{B}}\backslash{\mathcal{C}}_{l}$,
which exists in the solid part of materials to describe the degradation of
elastic stored contribution in energy. In order to extend
$\texttt{P}^{s}({\bm{\mathnormal{x}}})$ to the entire domain ${\mathcal{B}}$,
we define $\widehat{\texttt{P}}^{s}({\bm{\mathnormal{x}}},d)$ as:

$\texttt{P}^{s}({\bm{\mathnormal{x}}})\in{\mathcal{B}}\backslash{\mathcal{C}}_{l}\quad\xrightarrow{\quad}\quad\widehat{\texttt{P}}^{s}({\bm{\mathnormal{x}}},d)\in{\mathcal{B}}\quad\mathchoice{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}\quad\widehat{\texttt{P}}^{s}({\bm{\mathnormal{x}}},d)={\mathcal{G}}^{s}(d)\texttt{P}^{s}({\bm{\mathnormal{x}}})\;,$
(53)

with decreasing crack phase-field. Here, ${\mathcal{G}}^{s}(d)$ has a
descending order in terms of the crack phase-field. This set of functions can
be used to describe, e.g., heat permeability, fluid permeability, or in our
study, the permittivity and the permeability of the ferromagnetic medium. More
precisely, all sets of solid response are a subset of
$\widehat{\texttt{P}}^{s}({\bm{\mathnormal{x}}},d)$ such that in the limit
case this state function
$\widehat{\texttt{P}}^{s}({\bm{\mathnormal{x}}},d\rightarrow
0)\rightarrow\texttt{P}^{s}({\bm{\mathnormal{x}}})$ thus resembling classical
constitutive state variables for the solid response.

The second set of the formulation is denoted as
$\texttt{P}^{d}({\bm{\mathnormal{x}}})\in{\mathcal{C}}_{l}$, which exists in
the damaged domain to describe the fractured constitutive response. To further
extend $\texttt{P}^{f}({\bm{\mathnormal{x}}})$ to the entire domain, we
define:

$\texttt{P}^{f}({\bm{\mathnormal{x}}})\in{\mathcal{C}}_{l}\quad\xrightarrow{\quad}\quad\widehat{\texttt{P}}^{f}({\bm{\mathnormal{x}}},d)\in{\mathcal{B}}\quad\mathchoice{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}\quad\widehat{\texttt{P}}^{f}({\bm{\mathnormal{x}}},d)={\mathcal{G}}^{f}(d)\texttt{P}^{f}({\bm{\mathnormal{x}}})\;,$
(54)

with increasing crack phase-field. Here, ${\mathcal{G}}^{f}(d)$ has an
ascending order in terms of the crack phase-field, which has a inverse effect
of ${\mathcal{G}}^{s}(d)$ into the constitutive response. As a result, we have
the following properties:

$\begin{aligned} &\texttt{solid
part:}\quad\quad\;\;{\mathcal{G}}^{s}(d\rightarrow 0)\rightarrow
1,\quad{\mathcal{G}}^{s}(d\rightarrow 1)\rightarrow
0\quad\mathchoice{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}\quad\partial_{d}{\mathcal{G}}^{s}\leq
0\\\ &\texttt{fracture part:}\quad{\mathcal{G}}^{f}(d\rightarrow 0)\rightarrow
0,\quad{\mathcal{G}}^{f}(d\rightarrow 1)\rightarrow
1\quad\mathchoice{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}\quad\partial_{d}{\mathcal{G}}^{f}\geq
0\end{aligned}\;\;\;.$ (55)

Thus, ${\mathcal{G}}^{f}(d=0)$ is zero in the solid domain, while
${\mathcal{G}}^{s}(d=1)$ is zero in the fracture domain. In this contribution,
we propose the following formulation for ${\mathcal{G}}^{s}(d)$ and
${\mathcal{G}}^{f}(d)$ of order $m$ for a given set of parameters
$(c_{1},c_{2})$, for the solid contribution as:

$\displaystyle{\mathcal{G}}^{s}(d)=\Big{(}\frac{c_{2}-d}{c_{2}-c_{1}}\Big{)}^{m}\quad\mathchoice{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}\quad{\mathcal{G}}^{s}(d)=1\;\;\text{if}\;\;d<c_{1}\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad{\mathcal{G}}^{s}(d)=0\;\;\text{if}\;\;d>c_{2}\;,$
(56)

and the fracture contribution as:

$\displaystyle{\mathcal{G}}^{f}(d)=\Big{(}\frac{d-c_{1}}{c_{2}-c_{1}}\Big{)}^{m}\quad\mathchoice{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}{\hbox{{with}}}\quad{\mathcal{G}}^{f}(d)=0\;\;\text{if}\;\;d<c_{1}\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad{\mathcal{G}}^{f}(d)=1\;\;\text{if}\;\;d>c_{2},$
(57)

(a) (b) (c)  Figure 5:  (a) Degradation function for solid material denoted as
${\mathcal{G}}^{s}(d)$ (b) with representation the fracture domain if
$d>c_{2}$ and unfracture domain if $d<c_{1}$, and (c) Degradation function for
regularized crack phase-field ${\mathcal{G}}^{f}(d)$. Here, we used
$m=\\{1,2,3,4\\}$ with adjustable constants $(c_{1},c_{2})=(0.1,0.9)$.

which is in line properties (55). In the following numerical study we use the
quadratic order for both ${\mathcal{G}}^{s}(d)$ and ${\mathcal{G}}^{f}(d)$, so
$m=2$, and we further set $(c_{1},c_{2})=(0.1,0.9)$, see Figure 5.

######  Remark 2.2.

For the linear order $m=1$ our proposed degradation function is similar with
linear indicator functions given in [78], and for the quadratic order $m=2$
with constants $(c_{1},c_{2})=(0,1)$ it mimics quadratic indicator functions
given in [77]. Hence, (56) and (57) are generalized formulation versus [77]
and [78] .

To denote the effect of the magnetization-induced cracking in the material, we
propose an anisotropic electromagnetic material constants in Maxwell’s
equations. This provides a simple constitutive assumption for an
electromagnetic material to represent the transition rule. These constants can
be re-formulated through crack phase-field $d$. Thus, by means of state
functions for solid and fracture regions described in (56) and (57), the
definition of the electromagnetic material constants is then proposed as
follows:

$\begin{array}[]{ll}&\texttt{permeability:}\quad{\mu_{0}}(d)={\mathcal{G}}^{s}(d)\mu^{s}_{0}+{\mathcal{G}}^{f}(d)\mu^{f}_{0},\\\
\\\
&\texttt{permittivity:}\quad{\epsilon_{0}}(d)={\mathcal{G}}^{s}(d)\epsilon^{s}_{0}+{\mathcal{G}}^{f}(d)\epsilon^{f}_{0},\\\
\\\
&\texttt{conductivity:}\quad{\sigma_{0}}(d)={\mathcal{G}}^{s}(d)\sigma^{s}_{0}+{\mathcal{G}}^{f}(d)\sigma^{f}_{0},\\\
\\\ &\texttt{electric charge
density:}\quad{\rho_{0}}(d)={\mathcal{G}}^{s}(d)\rho^{s}_{0}+{\mathcal{G}}^{f}(d)\rho^{f}_{0}.\end{array}$
(58)

The first term in the anisotropic materials expression in (58) for a set of
$(\mu_{0},\epsilon_{0},\sigma_{0},\rho_{0})$ can be indicated as a classical
isotropic intrinsic material property for permeability, permittivity,
conductivity, and electric charge density, respectively. Also, in general, the
number of constants
$(\mu^{f}_{0},\epsilon^{f}_{0},\sigma^{f}_{0},\rho^{f}_{0})$ are much less
than material constants corresponding to solid parts and greater than the
number of constants describing vacuum counterparts. We note that in our study,
the second term in (58) is also a material property, which one can relate to
classical crack jump, e.g., in failure mechanics of a thermo-elastic solid in
[77], or in hydraulic fracturing of fluid-saturated porous media [79, 64, 65].

##  3 .Energy quantities and variational principles

To outline the variational formulation setting, it suffices to define the
constitutive energy density functions $W_{elas}$, $W_{mag}$, $W_{mos}$, and
$W_{frac}$ corresponds to elastic, magnetostrictive, magnetization, and
fracture contributions, respectively. These energy quantities lead to
establishing the multi-field evolution problem in terms of the primary fields
described in Section 2. To this end, the related constitutive relations are
provided. On this basis, a coupling between the magnetostriction and the
mechanics is formulated to investigate the mechanical deformation of materials
under the magnetic forces. Furthermore, the cracking response under the
magnetic forces is examined. The energetic crack driving force is governed by
the coupling magneto-mechanical effects under the magneto-static state. We
extend the variational formulations of the coupled multi-physics system within
the magneto-mechanical formulation of fracturing solids. We complete this
section by providing a compact algorithmic framework for magnetostrictive-
induced cracking.

###  3.1 .Constitutive functions

The coupled BVP is formulated by introducing three specific fields for
magneto-mechanically-induced cracking of the magneto-active materials

$\mbox{Global Primary Fields}:\ {\Bgothic
U}:=\\{{\bm{\mathnormal{u}}},d,{\bm{\mathnormal{A}}}\\}.$ (59)

Here, ${\bm{\mathnormal{u}}}$ is the displacement (mechanical deformation),
${\bm{\mathnormal{A}}}$ denotes the magnetic vector potential field, and $d$
is the crack phase-field ($0\leq d\leq 1$). From the numerical implementation
standpoint, to guarantee $0\leq d\leq 1$ holds, we project $d>1$ to 1 and
$d<0$ to 0 to avoid unphysical crack phase-field solution [61]. The
constitutive equations for the magnetostrictive phase-field fracture are
written in terms of the set

$\mbox{Constitutive State Variables}:\ {\Bgothic
C}:=\\{{\bm{\varepsilon}},d,\nabla d,{\bm{\mathnormal{B}}}\\},$ (60)

representing a combination of a magneto-mechanical model through the Maxwell’s
equation and a first-order gradient damage model.

####  3.1.1 .Energy quantities.

A pseudo-energy density per unit volume is then defined as $W:=W({\Bgothic
C})$, which is additively decomposed into an elastic contribution $W_{elas}$,
a magnetization contribution $W_{mag}$, a magnetostriction contribution
$W_{mos}$, and a (regularized) fracture contribution $W_{frac}\;$ through:

$\boxed{W({\Bgothic
C}):={W}_{elas}({\bm{\varepsilon}},d)+{W}_{mag}({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)+{W}_{mos}({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)+{W}_{frac}(d,\nabla
d).}$ (61)

We note that $W({\Bgothic C})$ is a state function that contains both
energetic and dissipative contributions. With this function at hand, a pseudo
potential energy functional can be constructed

$\displaystyle{\mathcal{E}}({\bm{\mathnormal{u}}},{\bm{\mathnormal{A}}},d)$
$\displaystyle:=\int_{\Omega}{W}({\Bgothic C})\,\text{d}{\bm{\mathnormal{x}}}$
(62)
$\displaystyle=\int_{{\mathcal{B}}}\Big{[}{{W}_{elas}({\bm{\varepsilon}},d)}+{W}_{frac}(d,\nabla
d)\Big{]}\,\text{d}{\bm{\mathnormal{x}}}+\int_{\Omega}\Big{[}{W}_{mag}({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)+{W}_{mos}({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)\Big{]}\,\text{d}{\bm{\mathnormal{x}}}.$

All these counterparts will be elaborated in following sections.

####  3.1.2 .Magnetization contribution.

Magnetization is defined as a vector field describing the density of permanent
or induced magnetic dipole moments in a magnetic material [80]. The magnetic
dipole moments are caused by either microscopic electric currents resulting
from the motion of electrons in atoms or the spin of the electrons or the
nuclei. In a magnetic field, a strong magnetization response is observed in
ferromagnetic materials. Any ferromagnetic material is capable of being
magnetized in the absence of an external magnetic field and becoming a
permanent magnet. It is not necessary that the magnetization distributes
uniformly within the material, so the material may present an anisotropic
response when subjected to the magnetization effect. Ferromagnetic materials
have inherent strong coupled magnetic and mechanical behavior, so the material
magnetization affects the stress response of material via magnetostriction
phenomenon, and the stress state of the material also affects magnetization
via inverse magnetostriction [81]. Following [16], the magnetization density
function is formulated as:

${W}_{mag}(I_{1},I_{4},d):=\frac{1}{2}\sum^{4}_{i=0}\frac{g_{i}(I_{1})}{i+1}\left(\frac{I_{4}}{B^{2}_{\text{ref}}}\right)^{i}I_{4},$
(63)

with the fourth invariant
$I_{4}={\bm{\mathnormal{B}}}\cdot{\bm{\mathnormal{B}}}$.

Here, $B_{\text{ref}}$ is a reference value of the magnetic flux density.
Additionally, functions $g_{i}(I_{1})$ for $i=\\{0,1,2,3,4\\}$ can be further
approximated, as it is shown in Appendix A leading to:

$\displaystyle g_{0}$
$\displaystyle=\frac{3}{4}\alpha_{0}\exp\left(\frac{3}{4}I_{1}\right)-\frac{1}{3}(\frac{1}{\mu_{0}}-\alpha_{5})$
(64) $\displaystyle g_{i}$
$\displaystyle=\frac{3(i+1)}{4}\alpha_{i}\exp\left(\frac{4(i+1)}{3}I_{1}\right)$
$\displaystyle\text{where}\quad i=\\{1,2,3,4\\}.$

The parameters $\alpha_{i}$ for $i=\\{0,1,2,3,4,5\\}$ refer to the
ferromagnetic structure due to the magnetization effect, see [16].

####  3.1.3 .Magnetostrictive contribution.

Magnetostriction is a constitutional characteristic of magnetic material
related to electron spin or orbit orientations, their interactions, and
molecular lattice geometry. The magnetostriction effect does not degrade in
time or during usage. A magnetostrictive material does not present the
magnetostriction effect when it is heated above its Curie temperature [82].
Curie temperature is defined as a transition point between ferromagnetic and
paramagnetic behaviors. In magnetic materials, three types of magnetostriction
mechanisms are determined to express the magnetostrictive feature of the
magnetic material. These mechanisms are paramagnetic magnetostriction,
spontaneous magnetostriction, and magnetostriction induced by magnetic field
[83]. Paramagnetism is a weak form of magnetism that is induced by an external
field in the direction of the applied magnetic field. It will disappear when
the magnetic field is removed. The spontaneous magnetization is followed by
paramagnetism when a ferromagnetic material is cooled below its Curie
temperature, so a transition from paramagnetism to ferromagnetism occurs, and
magnetic moments cause a spontaneous magnetization effect [84]. In a magnetic
material, the action of the external magnetic field [83] induces a
magnetostriction effect by rotating and moving the internal magnetic domain
walls.

We first employ two additional deformation-magnetization-dependent scalar-
valued function (invariants) to represent magnetostrictive effect through:

$I_{5}={\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}\cdot{\bm{\mathnormal{B}}}\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad
I_{6}={\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}^{2}\cdot{\bm{\mathnormal{B}}}.$
(65)

The fifth and sixth invariants account for the magnetostrictive response as a
function of the magnetic flux density ${\bm{\mathnormal{B}}}$ and the
mechanical strain ${\bm{\varepsilon}}$. The magnetostrictive density function
is formulated as a function of $(I_{5},I_{6})$ by:

${W}_{mos}(I_{5},I_{6},d):=\frac{1}{2}\alpha_{5}I_{5}+\frac{1}{2}\alpha_{6}I_{6}.$
(66)

Now we are able to define the magnetization vector ${\bm{\mathnormal{M}}}$ as
follows:

$\displaystyle{{\bm{\mathnormal{M}}}}({\bm{\varepsilon}},{\bm{\mathnormal{B}}};d):=-\frac{{\partial
W({\bm{\mathnormal{B}}},{\bm{\varepsilon}})}}{\partial{\bm{\mathnormal{B}}}}$
$\displaystyle=-\sum_{i}\frac{{\partial
W({\bm{\varepsilon}},{\bm{\mathnormal{B}}})}}{\partial I_{i}}\frac{{\partial
I_{i}}}{\partial{\bm{\mathnormal{B}}}}$ (67) $\displaystyle=-\frac{{\partial
W_{mag}({\bm{\varepsilon}},{\bm{\mathnormal{B}}})}}{\partial
I_{4}}\frac{{\partial
I_{4}}}{\partial{\bm{\mathnormal{B}}}}-\sum^{6}_{i=5}\frac{{\partial
W_{mos}({\bm{\varepsilon}},{\bm{\mathnormal{B}}})}}{\partial
I_{i}}\frac{{\partial I_{i}}}{\partial{\bm{\mathnormal{B}}}}$
$\displaystyle=-\sum_{i=0}^{4}\frac{g_{i}}{i+1}\partial_{I_{4}}\frac{I_{4}^{i+1}}{B^{2i}_{\text{ref}}}-\frac{1}{2}\left(\alpha_{5}{\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}-\alpha_{6}{\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}^{2}\right).$

For a detailed derivation, see Appendix A.

####  3.1.4 .Elastic contribution.

Since, the fracturing material behaves quite differently in bulk and shear
parts of the domain, we employ a consistent split for the strain energy
density function (68) into tension and compression counterparts, respectively.
Following additive split for strain tensor in (17), the effective strain
energy density $\psi_{e}$ admits the additive decomposition

$\psi_{e}\big{(}I_{1},I_{2}\big{)}:=\psi_{e}^{vol}(I_{1})+\psi_{e}^{dev}(I_{1},I_{2}),$
(68)

with

$\displaystyle\psi_{e}^{vol}\big{(}I_{1}\big{)}=\frac{K_{n}}{2}I^{2}_{1}\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad\psi_{e}^{dev}\big{(}I_{1},I_{2}\big{)}=\mu\Big{(}\frac{I_{1}^{2}}{3}-I_{2}\Big{)},\quad\text{note}\quad
K_{n},\mu>0.$ (69)

Here, $K_{n}=\lambda+\frac{2}{3}\mu$ is the bulk modulus which includes
elastic Lamé’s constant $\lambda$ and the shear modulus $\mu$. Additionally,
$I_{1}$ and $I_{2}$ denote the invariants

$I_{1}=\text{tr}[{\bm{\varepsilon}}]\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad
I_{2}=\mathop{\operator@font tr}[{\bm{\varepsilon}}^{2}].$ (70)

The isotropic strain energy density function given in (68) is additively
decomposed into tension and compression contributions:

$\psi_{e}\big{(}I_{1},I_{2}\big{)}=\psi_{e}^{+}(I_{1},I_{2})+\psi_{e}^{-}(I_{1}),$
(71)

where

$\displaystyle{\psi_{e}^{+}}(I_{1},I_{2})={H{{}^{+}}[I_{1}]}\psi_{e}^{vol}\big{(}I_{1}\big{)}+\psi_{e}^{dev}\big{(}I_{1},I_{2}\big{)}~{}\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad{\psi_{e}^{-}}(I_{1})=\big{(}1-{H{{}^{+}}[I_{1}]}\big{)}\psi_{e}^{vol}\big{(}I_{1}\big{)}~{}.$
(72)

Therein, $H{{}^{+}}[I_{1}({\bm{\varepsilon}}^{e})]$ is a positive Heaviside
function which returns one if $I_{1}({\bm{\varepsilon}}^{e})>0$ and zero if
$I_{1}({\bm{\varepsilon}}^{e})\leq 0$. The total elastic contribution to the
pseudo-energy density in (61) finally reads

${W}_{elas}({\bm{\varepsilon}},d):=g_{e}(d)\;\big{[}\psi_{e}^{+}(I_{1},I_{2})\big{]}+\psi_{e}^{-}(I_{1}),$
(73)

where $g_{e}(d)$ is a elastic degradation function. Here, the standard
monotonically decreasing quadrature degradation function, reads as
$g_{e}(d):=(1-\kappa)(1-d)^{2}+\kappa$. Here, $\kappa$ is a very small
parameter (to avoid numerical instabilities), and mathematically it is also
dependents on the discretization space, see [85].

We note that, the current resulting degraded stress tensor by the bulk energy
is rely on the $vol/dev$ decomposition. But, it is also possible to enhance
(73), for the material which are sensitive to the shear fracture as it is
investigated in rock-like materials, e.g. [86, 87]. So, to explore different
characteristic behavior of fracture which is subjected to compression and
shear modes for ferromagnetic materials, we could introduce different crack
phase-field driving force. This subject is open for further investigation.

The Cauchy stress tensor is defined as a derivative of the free energy
function with respect to the strain tensor as follows

$\displaystyle{{\bm{\sigma}}}({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d):=$
$\displaystyle\frac{{\partial
W}}{\partial{\bm{\varepsilon}}}=\sum_{i}\frac{{\partial W}}{\partial
I_{i}}\frac{{\partial I_{i}}}{\partial{\bm{\varepsilon}}}$ (74)
$\displaystyle=$
$\displaystyle\underbrace{g_{e}(d){\widetilde{{\bm{\sigma}}}}^{+}({\bm{\varepsilon}})+{\widetilde{{\bm{\sigma}}}}^{-}({\bm{\varepsilon}})}_{\text{elasticity
term}}+\underbrace{\frac{1}{2}\sum_{i=0}^{4}\frac{1}{i+1}\partial_{I_{1}}g_{i}\frac{I_{4}^{i+1}}{B^{2i}_{\text{ref}}}\textbf{I}}_{\text{magnetization
term}}$
$\displaystyle+\underbrace{\frac{1}{2}\alpha_{5}{\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}+\frac{1}{2}\alpha_{6}\left({\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}+{\bm{\varepsilon}}\cdot{\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}\right)}_{\text{magnetostrictive
term}}.$

with,

${\widetilde{{\bm{\sigma}}}}^{+}({\bm{\varepsilon}})=K_{n}H{{}^{+}}(I_{1})(\bm{\varepsilon}:\textbf{I})\textbf{I}+2\mu{\bm{\varepsilon}}^{dev},\quad\text{and}\quad{\widetilde{{\bm{\sigma}}}^{-}}({\bm{\varepsilon}})=K_{n}\big{(}1-H{{}^{+}}(I_{1})\big{)}(\bm{\varepsilon}:\textbf{I})\textbf{I}.$
(75)

In addition to the Cauchy stress tensor, we have a stress tensor induced by
the electromagnetic effects, denoted as ${\bm{\tau}}_{m}$ [88]. This
electromagnetic-induced stress tensor is expressed by the Lorentz force, which
furnishes an association between electromagnetism and the mechanics of the
material. In accordance to [16, 88, 17], the second-order electromagnetic
stress tensor defined as:

$\displaystyle{\bm{\tau}}_{m}({\bm{\mathnormal{B}}};d)=$
$\displaystyle\underbrace{\frac{E_{M}}{\mu_{0}}\Big{(}{\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}-\frac{1}{2}({\bm{\mathnormal{B}}}\cdot{\bm{\mathnormal{B}}})\textbf{I}\Big{)}+({\bm{\mathnormal{M}}}\cdot{\bm{\mathnormal{B}}})\textbf{I}-{\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{M}}}}_{\text{magnetic
response}}+$ (76)
$\displaystyle\underbrace{\epsilon_{0}\Big{(}{\bm{\mathnormal{E}}}\otimes{\bm{\mathnormal{E}}}-\frac{1}{2}({\bm{\mathnormal{E}}}\cdot{\bm{\mathnormal{E}}})\textbf{I}+{\bm{\mathnormal{E}}}\times{\bm{\mathnormal{B}}}\otimes\dot{{\bm{\mathnormal{x}}}}\Big{)}}_{\text{electric
response}}.$

We note that in the case of polarization effect, the above equation, i.e.,
(76), should also include this effect, see [88] Chapter 15. Here, $E_{M}>0$ is
a constant material property which is known as the magnetostrictive viscosity
constant [20]. Since, in the vacuum/air
(${\bm{\mathnormal{x}}}\in\Omega\backslash{\mathcal{B}}$) the magnetization
effect disappeares, thus results in ${\bm{\mathnormal{M}}}=\textbf{0}$, in the
electromagnetic stress (76). Additionally, since we only consider here for
ferromagnetic material thus, we only take into account magnetic response in
the quasi-static state, so the electric response in (76) is not considered.
This assumption is acceptable for many materials which only represent the
magnetic forces, e.g., iron; see for more detailed discussion [17, 16]. For
further information related to the derivation of the stress expression, which
contains the electromagnetically-dependent fields and parameters, you can
refer to [12, 11, 13, 23].

Finally, the reduced second-order stress tensor (by removing the effect of
electric response for the quasi-static state) is defined as

$\left\\{\,\vbox{\halign{\strut@$\m@th\displaystyle#$\hfil&\quad\strut@{#}\hfil\cr\hbox{\vrule
height=0.0pt,depth=0.0pt,width=0.0pt}$\displaystyle{\bm{\tau}}_{m}({\bm{\mathnormal{B}}};d)=\frac{E_{M}}{\mu_{0}}({\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}-\frac{1}{2}({\bm{\mathnormal{B}}}\cdot{\bm{\mathnormal{B}}})\textbf{I})+({\bm{\mathnormal{M}}}\cdot{\bm{\mathnormal{B}}})\textbf{I}-{\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{M}}}\quad\text{if}\;{\bm{\mathnormal{x}}}\in{\mathcal{B}}\\\\\\\\{\bm{\tau}}_{m}({\bm{\mathnormal{B}}};d)=\frac{E_{M}}{\mu_{0}}({\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}-\frac{1}{2}({\bm{\mathnormal{B}}}\cdot{\bm{\mathnormal{B}}})\textbf{I})\hskip
119.50148pt\quad\text{if}\;{\bm{\mathnormal{x}}}\in\Omega\backslash{\mathcal{B}}\crcr}}\right..$
(77)

Therefore, the total stress tensor is given in the solid domain
${\mathcal{B}}$ for a ferromagnetic material as:

${{\bm{\tau}}}({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)={{\bm{\sigma}}}({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)+{{\bm{\tau}}}_{m}({\bm{\mathnormal{B}}};d).$
(78)

By substituting the Cauchy stress tensor and the electromagnetic stress
considering the magnetization effect, the total stress tensor follows

$\displaystyle{{\bm{\tau}}}({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)=$
$\displaystyle
g_{e}(d){\widetilde{{\bm{\sigma}}}}^{+}({\bm{\varepsilon}})+{\widetilde{{\bm{\sigma}}}}^{-}({\bm{\varepsilon}})+\frac{1}{2}\sum_{i=0}^{4}\frac{1}{i+1}\partial_{I_{1}}g_{i}\frac{I_{4}^{i+1}}{B^{2i}_{\text{ref}}}\textbf{I}$
(79)
$\displaystyle+\left(E_{M}\mu_{0}^{-1}+\sum_{i=0}^{4}\frac{g_{i}}{i+1}\partial_{I_{4}}\frac{I_{4}^{i+1}}{B^{2i}_{\text{ref}}}+\frac{1}{2}\alpha_{5}\right){\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}$
$\displaystyle-\left(\frac{E_{M}}{2}\mu_{0}^{-1}+\sum_{i=0}^{4}\frac{g_{i}}{i+1}\partial_{I_{4}}\frac{I_{4}^{i+1}}{B^{2i}_{\text{ref}}}\right){\bm{\mathnormal{B}}}\cdot{\bm{\mathnormal{B}}}\textbf{I}$
$\displaystyle+\frac{1}{2}\alpha_{5}\Big{(}{\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}-({\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}\cdot{\bm{\mathnormal{B}}})\text{\bm{\mathnormal{I}}}\Big{)}$
$\displaystyle+\frac{1}{2}\alpha_{6}\Big{(}{\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}+{\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}^{2}-({\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}^{2}\cdot{\bm{\mathnormal{B}}})\text{\bm{\mathnormal{I}}}+({\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}+{\bm{\varepsilon}}\cdot{\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}})\Big{)}.$

Additionally the magnetic flux ${\bm{\mathnormal{H}}}$ can be specified

${\bm{\mathnormal{H}}}({\bm{\varepsilon}},{\bm{\mathnormal{B}}};d)=\frac{1}{\mu_{0}}{\bm{\mathnormal{B}}}-2\sum_{i=0}^{4}\frac{g_{i}}{i+1}\frac{\partial
I_{4}^{i+1}}{\partial
I_{4}}-\frac{1}{2}(\alpha_{5}{\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}+\alpha_{6}{\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}^{2}).$
(80)

######  Remark 3.1.

We note the traction force vector
$\bar{{\bm{\mathnormal{t}}}}={\bm{\tau}}\cdot{\bm{\mathnormal{n}}}$ is
generally decomposed through the mechanic, electric, and magnetic fields’ so-
called generalized traction vector. If we define $\bar{{\bm{\mathnormal{t}}}}$
as a mechanical traction force vector only only mechanical contribution
${\widetilde{{\bm{\sigma}}}}$ exist, while the traction force vector for the
electric and magnetic fields will be neglected, and set as zero at Neumann
boundary conditions [89]. Additionally, the traction force due to the electric
field is discussed in [90], which one can evaluate the electrostatic traction
through the so-called Maxwell stress tensor. Also, regarding the magnetic
force vector, the resultant magnetic force is divided into magnetic traction
and body forces [91]. Let’s define traction vector as the interaction between
two adjacent surface elements. Traction magnetic force is computed based on
potential magnetic force inside and outside of the magnetizable body. Thus,
magnetic traction can be seen as the jumping value of a physical quantity
between two sides of the surface [91]. So, the Neumman boundary conditions can
be set for the magnetization jump between two adjacent magnetized surfaces.
The traction force in magnetic separators is also investigated in [92].
Additionally, the electromagnetic traction due to electromagnetic stress
tensor is described in [93]. They discussed whether the electromagnetic
effects enter as a body force or as a surface stress, but not both. This is
same as our case where we have only body force.

####  3.1.5 .Fracture contribution.

The fracture contribution of the pseudo-energy density given in (61) reads:

$\displaystyle{W}_{frac}(d,\nabla d):=G_{c}\gamma_{l}(d,\nabla d),$ (81)

where ${G}_{c}>0$ is so-called a Griffith’s energy release rate. In this
study, we use a small time increment $\Delta t=t_{n+1}-t_{n}$ where $t_{n}$
denotes the previous time step. This reduces the complexity of the problem. We
further assume that:

$\displaystyle\mu_{0}(d,t)\approx\mu_{0}(d,t_{n}),\;\;\epsilon_{0}(d,t)\approx\epsilon_{0}(d,t_{n}),$
(82)
$\displaystyle\sigma_{0}(d,t)\approx\sigma_{0}(d,t_{n}),\;\;\rho_{0}(d,t)\approx\rho_{0}(d,t_{n}).$

This type of assumption has already been employed in hydraulic phase-field
fracture. See, for instance, [94, 95], where the second-order anisotropic
permeability tensor is computed of the previous time $t_{n}$. Finally, by
taking the first variational derivative $\delta_{d}W$ of (61), the positive
crack driving state function ${\mathcal{D}}({\bm{\mathnormal{x}}},t)$ follows

${\mathcal{D}}:=\frac{2l_{d}}{G_{c}}\psi_{e}^{+}(I_{1},I_{2})\ .$ (83)

###  3.2 .Variational formulations for the coupled multi-field problem

The variational formulations with respect to the four PDEs described in the
previous sections for failure mechanics of ferromagnetic materials under the
magnetostrictive effects are further discussed. The rate-dependent PDEs
models, discussed in the previous section, are defined in a temporal domain
$[t_{n},t_{n+1}]$ with $\Delta t=t_{n+1}-t_{n}>0$ holds. For the time-
dependent problem in electromagnetic-induced fracture, we require initial
conditions at $t=0$ as:

$\displaystyle{\bm{\mathnormal{u}}}^{0}={\bm{\mathnormal{u}}}({\bm{\mathnormal{x}}},0),\quad\;d^{0}=d({\bm{\mathnormal{x}}},0),$
(84)
$\displaystyle{\bm{\mathnormal{E}}}^{0}={\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},0),\quad{\bm{\mathnormal{H}}}^{0}={\bm{\mathnormal{H}}}({\bm{\mathnormal{x}}},0),$
$\displaystyle{\bm{\mathnormal{D}}}^{0}={\bm{\mathnormal{D}}}({\bm{\mathnormal{x}}},0),\quad{\bm{\mathnormal{B}}}^{0}={\bm{\mathnormal{B}}}({\bm{\mathnormal{x}}},0),$

where ${\bm{\mathnormal{H}}}^{0}$ has to satisfy the divergence -free
assumption

$\nabla.\left(\mu_{0}{\bm{\mathnormal{H}}}^{0}({\bm{\mathnormal{x}}})\right)=0\quad\text{in}\;\Omega\;\;\;\;\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad\;\;\;\;{\bm{\mathnormal{H}}}^{0}\cdot{\bm{\mathnormal{N}}}=0\quad\text{on}\;\partial\Omega\;.$
(85)

####  3.2.1 .Electromagnetic contribution.

The magnetostatic response is a special case of the electromagnetic problem
which is obtained when the electric field
${\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)$ does not change by time, and
thus ${\partial{\bm{\mathnormal{D}}}}/{\partial t}=0$. This results in

$\frac{\partial{{\bm{\mathnormal{D}}}({\bm{\mathnormal{x}}},t)}}{\partial{t}}=\epsilon_{0}\frac{\partial{{\bm{\mathnormal{E}}}({\bm{\mathnormal{x}}},t)}}{\partial{t}}=\bm{0}.$
(86)

The variation of the electric displacement in time is 0, and thus we can
assume that permittivity is zero, i.e., $\epsilon_{0}=0$. The time-discretized
variational formulation of the evolution equation of the magnetic vector
potential ${\bm{\mathnormal{A}}}({\bm{\mathnormal{x}}},t)$ in (42) is derived
by multiplying it with the test function $\delta{{\bm{\mathnormal{A}}}}$ and
then applying integrating by part

$\displaystyle\mathcal{E}_{{\bm{\mathnormal{A}}}}({\Bgothic
U},\delta{\bm{\mathnormal{A}}})$ $\displaystyle=\Delta
t\int_{\Omega}\frac{1}{\mu_{0}(d)}\nabla{\bm{\mathnormal{A}}}:\nabla(\delta{\bm{\mathnormal{A}}})\,\text{d}{\bm{\mathnormal{x}}}+\int_{\Omega}\sigma_{0}(d)\left[({\bm{\mathnormal{A}}}-{\bm{\mathnormal{A}}}^{n})+\Delta
t\nabla\Phi\right]\cdot\delta{\bm{\mathnormal{A}}}\,\text{d}{\bm{\mathnormal{x}}}$
(87) $\displaystyle-\Delta
t\int_{\Omega}{\bm{\mathnormal{J}}}_{s}(t)\cdot\delta{\bm{\mathnormal{A}}}\,\text{d}{\bm{\mathnormal{x}}}=0\qquad\qquad\forall\,\delta{\bm{\mathnormal{A}}}\in{\mathcal{S}}_{0}^{{\bm{\mathnormal{v}}}}(\text{curl};\Omega)\;,$

where (42) is imposed.

####  3.2.2 .Elastic contribution.

The variational formulations of the strong formulation of the Euler-Lagrange
equations with respect to the ${\bm{\mathnormal{u}}}({\bm{\mathnormal{x}}},t)$
in Formulation 2.a reads:

$\mathcal{E}_{\bm{\varepsilon}}({\Bgothic
U},\delta{\bm{\mathnormal{u}}})=\displaystyle\int_{\mathcal{B}}\Big{[}{\bm{\tau}}:\delta{\bm{\varepsilon}}-\bar{{\bm{\mathnormal{b}}}}\cdot\delta{\bm{\mathnormal{u}}}\Big{]}\text{d}{\bm{\mathnormal{x}}}-\int_{\partial^{{\bm{\mathnormal{u}}}}_{N}{\mathcal{B}}}\bar{{\bm{\mathnormal{t}}}}\cdot\delta{\bm{\mathnormal{u}}}\;\text{d}s=0\quad\forall\,\delta{\bm{\mathnormal{u}}}\in{\mathcal{W}}_{0}^{{\bm{\mathnormal{u}}}}(\text{div};{\mathcal{B}}).$
(88)

We note that $\mathcal{E}_{\bm{\varepsilon}}({\Bgothic
U},\delta{\bm{\mathnormal{u}}})$ unlike the magnetic vector potential given in
(87), lives in ${\mathcal{B}}$ (and not $\Omega$). Since ${\bm{\tau}}$ is a
function of ${\bm{\mathnormal{B}}}$ (and not ${\bm{\mathnormal{A}}}$), the
local constitutive equation
${\bm{\mathnormal{B}}}=\nabla\times{\bm{\mathnormal{A}}}$ has to be imposed
after solving magnetic vector potential from (87).

####  3.2.3 .Fracture contribution.

The variational formulation of the strong formulation of fracture with respect
to the crack phase-field $d({\bm{\mathnormal{x}}},t)$ in Formulation 2.a
reads:

$\displaystyle\mathcal{E}_{d}({\Bgothic U},\delta d)$
$\displaystyle=(1-\kappa)\Delta
t\int_{{\mathcal{B}}}(d-1){\mathcal{H}}\cdot\delta
d\;\text{d}{\bm{\mathnormal{x}}}+\int_{{\mathcal{B}}}\Big{[}\Delta
td\cdot\delta d+\eta_{d}(d-d^{n})\cdot\delta
d\Big{]}\;\text{d}{\bm{\mathnormal{x}}}$ (89)
$\displaystyle+\int_{{\mathcal{B}}}\Big{[}l_{d}^{2}\Delta t\nabla
d\cdot\nabla(\delta
d)\Big{]}\,\text{d}{\bm{\mathnormal{x}}}=0\quad\quad\forall\,\delta
d\in{\mathcal{W}}_{0}^{d}(\text{div};{\mathcal{B}}).$

As in the weak form $\mathcal{E}_{\bm{\varepsilon}}({\Bgothic
U},\delta{\bm{\mathnormal{u}}})$, the crack phase-field equation lives in
${\mathcal{B}}$.

Thus, by means of (87), (88), and (89), the fully coupled variational multi-
field problem to describing magnetostrictive induced fractures in
ferromagnetic material is formulated in the compact form as:

$\displaystyle\mathcal{E}_{{\Bgothic U}}({\Bgothic U},\delta{\Bgothic
U})=\mathcal{E}_{\bm{\varepsilon}}({\Bgothic U},\delta{\bm{\mathnormal{u}}})$
$\displaystyle+\mathcal{E}_{A}({\Bgothic
U},\delta{\bm{\mathnormal{A}}})+\mathcal{E}_{d}({\Bgothic U},\delta d)=0$ (90)
$\displaystyle\forall\;(\delta{\bm{\mathnormal{u}}},\delta{\bm{\mathnormal{A}}},\delta
d)\in\Big{(}{\mathcal{W}}_{0}^{{\bm{\mathnormal{u}}}}(\text{div};{\mathcal{B}}),{\mathcal{S}}_{0}^{{\bm{\mathnormal{v}}}}(\text{curl};\Omega),{\mathcal{W}}_{0}^{d}(\text{div};{\mathcal{B}})\Big{)}.$

In order to solve the magnetostrictive induced fracture system (90), we first
solve the first two equations monolithically (simultaneously obtain
$({\bm{\mathnormal{u}}},{\bm{\mathnormal{A}}})$). Then, a staggered approach
is employed to obtain the phase-field fracture $d$. To that end, we fix
alternately $(\bm{u},{\bm{\mathnormal{A}}})$ and estimate $d$ and vice versa.
The procedure is continued until convergence (using a given
$\texttt{TOL}_{\mathrm{Stag}}$). The alternate minimization scheme applied to
the (90) is summarized in Algorithm 1.

Input: loading data $\bar{\bm{t}}_{n}$ on
$\partial^{{\bm{\mathnormal{u}}}}_{N}{\mathcal{B}}$ and
${\bm{\mathnormal{J}}}_{s}(t)$ on wires [2mm] solution
$(\bm{u}^{n-1},{\bm{\mathnormal{A}}}^{n-1},d^{n-1})$ from step $n-1$. [2mm]
Initialization, $k=1$:

• set
$(\bm{u}^{0},{\bm{\mathnormal{A}}}^{0},d^{0}):=(\bm{u}^{n-1},{\bm{\mathnormal{A}}}^{n-1},d^{n-1})$.

Staggered iteration between $(\bm{u},{\bm{\mathnormal{A}}})$ and $d$:

• solve following system of equations in a monolithic manner by given
$d^{k-1}$,

$\displaystyle\left\\{\,\vbox{\halign{\strut@$\m@th\displaystyle#$\hfil&\quad\strut@{#}\hfil\cr\hbox{\vrule
height=0.0pt,depth=0.0pt,width=0.0pt}$\displaystyle\Delta
t\int_{\Omega}\frac{1}{\mu_{0}(d)}\nabla{\bm{\mathnormal{A}}}:\nabla(\delta{\bm{\mathnormal{A}}})\,\text{d}{\bm{\mathnormal{x}}}+\int_{\Omega}\sigma_{0}(d)\left[({\bm{\mathnormal{A}}}-{\bm{\mathnormal{A}}}^{n})+\Delta
t\nabla\Phi\right]\cdot\delta{\bm{\mathnormal{A}}}\,\text{d}{\bm{\mathnormal{x}}}-\Delta
t\int_{\Omega}{\bm{\mathnormal{J}}}_{s}(t)\cdot\delta{\bm{\mathnormal{A}}}\,\text{d}{\bm{\mathnormal{x}}}=0,\\\\\\\\\displaystyle\int_{\mathcal{B}}\Big{[}{\bm{\tau}}({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d):\delta{\bm{\varepsilon}}-\bar{{\bm{\mathnormal{b}}}}\cdot\delta{\bm{\mathnormal{u}}}\Big{]}\text{d}{\bm{\mathnormal{x}}}-\int_{\partial^{{\bm{\mathnormal{u}}}}_{N}{\mathcal{B}}}\bar{{\bm{\mathnormal{t}}}}\cdot\delta{\bm{\mathnormal{u}}}\;\text{d}s=0\quad
with\;\;{\bm{\mathnormal{B}}}=\nabla\times{\bm{\mathnormal{A}}},\crcr}}$

for $(\bm{u},{\bm{\mathnormal{A}}})$, set
$(\bm{u},{\bm{\mathnormal{A}}})=:(\bm{u}^{k},{\bm{\mathnormal{A}}}^{k})$,

• given $(\bm{u}^{k},{\bm{\mathnormal{B}}}^{k})$, solve

$\displaystyle(1-\kappa)\Delta
t\displaystyle\int_{{\mathcal{B}}}(d-1){\mathcal{H}}\cdot\delta
d\;\text{d}{\bm{\mathnormal{x}}}+\displaystyle\int_{{\mathcal{B}}}\Big{[}\Delta
td\cdot\delta d+\eta_{d}(d-d^{n})\cdot\delta
d\Big{]}\;\text{d}{\bm{\mathnormal{x}}}+\displaystyle\int_{{\mathcal{B}}}\Big{[}l_{d}^{2}\Delta
t\nabla d\cdot\nabla(\delta d)\Big{]}\,\text{d}{\bm{\mathnormal{x}}}=0$

for $d$, set $d=:d^{k}$,

• for the obtained pair $(\bm{u}^{k},{\bm{\mathnormal{A}}}^{k},d^{k})$, check
staggered residual by

$\displaystyle\quad\quad\;\;\mathrm{Res}_{\mathrm{Stag}}^{k}:|\mathcal{E}_{{\bm{\varepsilon}}}({\Bgothic
U}^{k};\delta{\bm{\mathnormal{u}}})|+|\mathcal{E}_{{\bm{\mathnormal{A}}}}({\Bgothic
U}^{k};\delta{\bm{\mathnormal{A}}})|\leq\texttt{TOL}_{\mathrm{Stag}}$

• if fulfilled, set ${\Bgothic
U}^{k}=(\bm{u}^{k},{\bm{\mathnormal{A}}}^{k},d^{k})=:(\bm{u}^{n},{\bm{\mathnormal{A}}}^{n},d^{n})={\Bgothic
U}^{n}$ then stop;

• else $k+1\rightarrow k$. [2mm] Output: solution
$(\bm{u}^{n},{\bm{\mathnormal{A}}}^{n},d^{n})$ at $n^{\text{th}}$ time-step.
[2mm]

######  Remark 3.2.

Since the permeability of the ferromagnetic material is computed as a function
of the crack phase-field with a degrading function, the magnetic vector
potential ${\bm{\mathnormal{A}}}$, the magnetic field ${\bm{\mathnormal{B}}}$
and the magnetic flux ${\bm{\mathnormal{H}}}$ are implicitly determined as a
function of the crack phase-field. Therefore, when a crack starts to propagate
in the material, the degradation function actively represents its degrading
effect on the magnetic response of the material. Then, the magnetic vector
potential ${\bm{\mathnormal{A}}}$, the magnetic field ${\bm{\mathnormal{B}}}$
and the magnetic flux ${\bm{\mathnormal{H}}}$ decrease along the crack paths.

Algorithm 1 The staggered iterative solution process for (90) at a fixed time-
step $n$.

###  3.3 .Space finite element discretization

In this section, we provide spatial discretizations of the variational forms
given in Algorithm 3.2.3. This results in a discretized multi-field problem to
be solved for three-field unknowns represented by
$({\bm{\mathnormal{u}}},{\bm{\mathnormal{A}}},d)$. In the case of the
magnetostatic problem, one can use either Nédélec elements or alternatively
$H^{1}(\Omega)$ conforming elements with continuous piecewise polynomials, see
[96]. We use a Galerkin finite element method to discretize the equations,
employing second order isoparametric curl-conforming triangular element
$P_{2}$ for the ansatz and test spaces of all primary fields see [97, 98, 99,
100, 101, 102]. The continuous solid domains ${\mathcal{B}}$, and the vacuum
domain $\Omega$ are approximated by ${\mathcal{B}}_{h}$, and $\Omega_{h}$ such
that ${\mathcal{B}}\approx{\mathcal{B}}_{h}$, and $\Omega\approx\Omega_{h}$.
The approximated domains ${\mathcal{B}}_{h}$, and $\Omega_{h}$ are decomposed
with non-overlapping linear triangular finite element
${\mathcal{B}}_{e}\subset{\mathcal{B}}_{h}$ such that

${\mathcal{B}}\approx{\mathcal{B}}_{h}=\bigcup_{e}^{n_{e}}{\mathcal{B}}_{e}\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad\Omega\approx\Omega_{h}=\bigcup_{e}^{n_{e}}\Omega_{e}.$
(91)

The unknown fields $({\bm{\mathnormal{u}}},{\bm{\mathnormal{A}}},d)$ are
approximated by

${\bm{\mathnormal{u}}}_{h}={\bm{\mathnormal{N}}}_{\bm{\mathnormal{u}}}\hat{\bm{u}},\quad
d_{h}=\bm{N}_{d}\hat{\bm{d}},\quad{\bm{\mathnormal{A}}}_{h}=\bm{N}_{\bm{\mathnormal{A}}}\hat{{\bm{\mathnormal{A}}}}.$
(92)

with a set of nodal solution values
$(\hat{\bm{u}},\hat{\bm{d}},\hat{{\bm{\mathnormal{A}}}})$. Their derivatives
are given by

$\bm{\varepsilon}(\bm{u}_{h})=\bm{B}_{u}\hat{\bm{u}},\quad\nabla
d_{h}=\bm{B}_{d}\hat{\bm{d}}\;,\quad\nabla{\bm{\mathnormal{A}}}_{h}=\bm{B}_{\bm{\mathnormal{A}}}\hat{\bm{A}}\;,$
(93)

which represent constitutive state variables
$({\bm{\varepsilon}}_{{\bm{\mathnormal{u}}}_{h}},\nabla
d_{h},\nabla{\bm{\mathnormal{A}}}_{h})$. Here,
${\bm{\mathnormal{B}}}_{{\bm{\mathnormal{u}}}}$, ${\bm{\mathnormal{B}}}_{d}$,
${\bm{\mathnormal{B}}}_{{\bm{\mathnormal{A}}}}$ are the matrix representation
of the shape function’s derivatives, corresponding to the global deformation,
crack phase-field, and the magnetic vector potential, respectively.
Additionally, the magnetic flux densities (which is enter to deformation
equation) are approximated by the curl operator as:

${\bm{\mathnormal{B}}}({\bm{\mathnormal{A}}}_{h})=\nabla\times{\bm{\mathnormal{A}}}_{h}=\bm{\Pi}\hat{{\bm{\mathnormal{A}}}}.$
(94)

Here, $\bm{\Pi}$ is the matrix representation of the shape function’s
derivative, see [16] (Section 3.2), and [103] for a detailed discussion. For
each primary field, the set of discretized weak forms which is based on the
residual force vector denoted by ${\bm{\mathnormal{R}}}_{\bullet}$ has to be
determined. First, the mechanical weak form leads to

${{\bf{R}}_{\bm{\mathnormal{u}}}}({\bm{\mathnormal{u}}},d,{\bm{\mathnormal{A}}})=\sum_{e}\int_{{\mathcal{B}}_{e}}\big{[}{\bm{\mathnormal{B}}}^{T}_{\bm{\mathnormal{u}}}{\bm{\tau}}({\bm{\varepsilon}}_{h},{\bm{\mathnormal{B}}}_{h},d_{h})-\overline{{\bm{\mathnormal{b}}}}\cdot{\bm{\mathnormal{N}}}_{\bm{\mathnormal{u}}}\big{]}\,\text{d}{\bm{\mathnormal{x}}}-\int_{\partial{\mathcal{B}}_{N}}\overline{{\bm{\mathnormal{t}}}}\cdot{\bm{\mathnormal{N}}}_{\bm{\mathnormal{u}}}\,ds={\bm{0}},$
(95)

next, for the magnetic vector potential yields,

$\displaystyle{{\bf{R}}_{\bm{\mathnormal{A}}}}({\bm{\mathnormal{u}}},d,{\bm{\mathnormal{A}}})=\sum_{e}\int_{\Omega_{e}}\bigg{[}\frac{\Delta
t}{\mu_{0}(d)}{\bm{\mathnormal{B}}}^{T}_{\bm{\mathnormal{A}}}\nabla{\bm{\mathnormal{A}}}$
$\displaystyle+\sigma_{0}(d)\big{[}({\bm{\mathnormal{A}}}-{\bm{\mathnormal{A}}}^{n})+\Delta
t\nabla\Phi\big{]}\cdot{\bm{\mathnormal{N}}}_{{\bm{\mathnormal{A}}}}$ (96)
$\displaystyle-\Delta
t{\bm{\mathnormal{J}}}_{s}(t)\cdot{\bm{\mathnormal{N}}}_{{\bm{\mathnormal{A}}}}\;\bigg{]}\,\text{d}{\bm{\mathnormal{x}}}={\bm{0}}\;,$

and lastly the crack phase-field follows as:

$\displaystyle{{\bf{R}}_{d}}({\bm{\mathnormal{u}}},d,{\bm{\mathnormal{A}}})=\sum_{e}\int_{{\mathcal{B}}_{e}}\bigg{[}\bigg{(}(1-\kappa)(d_{h}-1){\mathcal{H}}+{d}$
$\displaystyle+\frac{\eta_{d}}{\Delta
t}(d_{h}-d_{h}^{n})\bigg{)}\cdot{\bm{\mathnormal{N}}}_{d}$ (97)
$\displaystyle+l_{d}^{2}{\bm{\mathnormal{B}}}_{d}^{T}\cdot\nabla
d_{h}\bigg{]}\,\text{d}{\bm{\mathnormal{x}}}={\bm{0}}.$

##  4 .Numerical examples

In this section, the capability of the proposed framework in predicting the
magnetorestrictively-induced cracking in the ferromagnetic materials is
approved by demonstrating the numerical examples. At first, the evolution of
the magnetic vector potential and the magnetic field in a ferromagnetic
material is investigated. Then, cracking in the ferromagnetic materials is
investigated by adopting the magneto-mechanical model coupled with the phase-
field approach. The boundary-value problem of the first example is related to
a domain that includes copper wires and an iron beam that are surrounded by
vacuum. In the remaining examples, the domain of the boundary-value problems
contains predefined notches and copper wires where a constant electric current
source is enforced. These representative numerical examples are solved by
developing a finite element software based on the FEniCS software library, see
[104, 105]. The material parameters used in the numerical examples are listed
in Table LABEL:material-parameters. For the staggered approach, see Algorithm
1, we set $\texttt{TOL}_{\mathrm{Stag}}=10^{-4}$. We also note that the second
term of anisotropic material constants in (58) is assumed for quasi-
impermeable crack faces, hence we set
$(\mu^{f}_{0},\epsilon^{f}_{0},\sigma^{f}_{0},\rho^{f}_{0})$ identical to
vacuum counterparts.

Table 1: Material parameters used in the numerical examples based on [88] No. | Parameter | Name | Example 1-2 | Example 3 | Unit
---|---|---|---|---|---
1. | $E$ | Young’s modulus | $160$ | $16$ | $\mathrm{GPa}$
2. | $\nu$ | Poisson’s ratio | $0.33$ | $0.2$ | –
3. | $\mu_{0,\Omega}$ | Permeability of vacuum | $4\pi\times 10^{-7}$ | $--$ | $\mathrm{H/mm}$
4. | $\mu_{0,{\mathcal{W}}}$ | Permeability of wires | $1.26\times 10^{2}$ | $1.26\times 10^{3}$ | $\mathrm{H/mm}$
5. | $\mu_{0,{\mathcal{B}}}$ | Permeability in solid | $1000$ | $10^{7}$ | $\mathrm{H/mm}$
6. | $\sigma_{0}$ | Electromagnetic conductivity | $1$ | $1$ | $\mathrm{S/{mm}}$
7. | $B_{\text{ref}}$ | The magnetic flux density | $2\times 10^{6}$ | $2\times 10^{6}$ | $\mathrm{T}$
8. | $(\alpha_{0},...,\alpha_{3})$ | Magnetisation constants | $2\times 10^{-6}$ | $2\times 10^{-6}$ | –
9. | $(\alpha_{4},\alpha_{5})$ | Magnetostrictive constants | $(2.5,7.75)$ | $(2.5,7.75)$ | –
10. | $E_{M}$ | Magnetostrictive viscosity | $10^{-9}$ | $10^{-9}$ | –
11. | $\zeta$ | Permeability transition exponent | $50$ | $50$ | –
12. | $\eta_{d}$ | damage viscosity | $1\times 10^{-6}$ | $1\times 10^{-6}$ | $\mathrm{N/(mm^{2})s}$
13. | $l_{d}=2h_{e}$ | fracture regularization parameter | $2.196$ | $2.196$ | –
14. | $\kappa$ | Stabilization parameter | $10^{-8}$ | $10^{-8}$ | –
15. | $G_{c}$ | Griffith’s energy release rate | $0.0027$ | $0.0027$ |
16. | $(m,c_{1},c_{2})$ | degradation parameters | $(2,0.1,0.9)$ | $(2,0.1,0.9)$ | –

###  4.1 .Example 1: Magnetostatic problem for transient magnetic vector
potential

The first insight into the performance of the variational form of Maxwell’s
equations in a finite element context are gained by a two-dimensional
transient magnetostatic problem in a domain composed of copper wires and an
iron cylinder that are surrounded by a vacuum. The geometrical configuration
of the 2D test problem and its finite element model are depicted in Figure 6,
in which the subdomains for the iron cylinder, the copper wires, and the
surrounding vacuum with the radius of $r_{v}=5$ mm are clearly visible. The
inner radius of the cylinder is $r_{1}=1$ mm and its thickness is $n=0.2$ mm.
The domain contains ten copper wires with a radius of $r_{w}=0.1$ mm. The
discretization of the domain is performed with $64652$ triangular elements
with the minimum element size of $h_{\text{min}}\approx 0.04$ mm.

The wires are assumed to carry the electric current, which is varied linearly
with time. The electromagnetic waves are neglected, since the problem is
supposed to be magnetostatics. To examine this problem, the Poisson-type
equation (37) is solved. The specific current ${\bm{\mathnormal{J}}}_{s}(t)$
along the $z$-direction is given by

exterior wires:
$\displaystyle{\bm{\mathnormal{J}}}_{s}(t,z)=-t\;\;\text{A/mm${}^{2}$},\quad\texttt{interior
wires:}\;\;$
$\displaystyle{\bm{\mathnormal{J}}}_{s}(t,z)=t\;\;\text{A/mm${}^{2}$},$

for the interior set of the winding copper wires. The results are obtained in
terms of the magnetic vector potential in $z$-direction,
$A_{z}({\bm{\mathnormal{x}}},t)$, for the different time steps, and the
magnetic flux density ${\bm{\mathnormal{B}}}({\bm{\mathnormal{x}}},t)$ for the
last time step. These results are shown in Figure 7 and Figure 8,
respectively.

(a) (b) (c)  Figure 6: Example 1. The representation of (a) whole domain of
the magnetostatic problem including copper wires winding around the iron
cylinder and vacuum, (b) geometry and dimensions of iron cylinder, and (c)
finite element mesh generated for the test problem. The subdomains for the
iron cylinder and copper wires are clearly demonstrated.  (a) (b) (c)  Figure
7: Example 1. The representation of the $z$-component of the magnetic vector
potential $A_{z}$ at (a) $t=40$ s (b) $t=80$ s (c) $t=120$ s.  (a) (b)  Figure
8: Example 1. The representation of (a) the $x-$ and (b) $y-$ components of
the magnetic field field ${\bm{\mathnormal{B}}}$ at $t=120$ s

###  4.2 .Example 2: Magneto-mechanically induced cracking in an iron beam
surrounded by vacuum

The second example investigates the magnetostrictive-induced cracking in a
ferromagnetic notched beam. This problem is solved by coupling the Poisson-
type transient magnetic equation with the phase-field model. The domain of the
boundary value problem consists of vacuum, two copper wires, and the single
notched simply-supported iron beam. The electric current flows through the top
and the bottom wires. The geometrical configuration and finite element model
of the boundary value problem are displayed in Figure 9, respectively. The
whole domain is composed of vacuum with the radius of $r_{v}=8$ mm, two copper
wires with the radius of $r_{w}=0.5$ mm, and the notched iron beam with
$n=0.8$ mm and $m=3.6$ mm located at the middle of the vacuum with the
dimension of $A=4$ mm and $B=2$ mm hence ${\mathcal{B}}=(0,4)\times(0,2)$ mm2.
The whole domain is discretized with $28325$ triangular elements. As the
domain is divided into three subdomains, and the crack phase-field equation is
only solved in the iron beam subdomain, the discretization is refined in that
region. In this case, the minimum element size is $h_{\text{min}}\approx
0.005$ mm and the maximum element size is $h_{\text{max}}\approx 1.098$ mm,
see Figure 9(c).

(a) (b) (c)  Figure 9: Example 2. The representation of geometries and
dimensions of (a) whole domain, (b) iron beam, and (c) finite element model of
BVP. The subdomains are clearly invisible.

The wires located on the top and bottom of the iron beam contain the electric
current. The specific current ${\bm{\mathnormal{J}}}_{s}(t)$ along the
$z$-direction is given in the form of a linear equation which varies in time
as wire.

top wire:
$\displaystyle{\bm{\mathnormal{J}}}_{s}(t,z)=2.5\,t\;\;\text{A/mm${}^{2}$},\quad\texttt{bottom
wire:}\;\;$
$\displaystyle{\bm{\mathnormal{J}}}_{s}(t,z)=-1.5\,t\;\;\text{A/mm${}^{2}$}.$

At first, the evolution of magnetic vector potential
$A_{z}({\bm{\mathnormal{x}}},t)$ and magnetic field
${\bm{\mathnormal{B}}}({\bm{\mathnormal{x}}},t)$ are calculated through
Maxwell’s equations in the whole domain. The stress response regarding the
magnetostrictive effects is determined in the iron beam. The electric current
flow in the wires generates magnetization in the notched iron beam that forces
it to deform. The deformation causes stress development around the notch.
Therefore, cracking starts to initiate when the maximum principal stress
exceeds the critical value. The crack initiation and propagation in the beam
at different time steps are demonstrated in Figure 10. The magnetic vector
potential $A_{z}({\bm{\mathnormal{x}}},t)$ and the magnetic field
${\bm{\mathnormal{B}}}({\bm{\mathnormal{x}}},t)$ at the last time step
$t=0.68$ s are exhibited in Figure 11.

Figure 10: Example 2. The representation of the crack initiation and
propagation in the iron beam under coupling electromagneto-mechanical effects
at (a) $t=0.43$ s, (b) $t=0.45$ s, (c) $t=0.56$ s, (d) $t=0.59$ s, (e)
$t=0.64$ s and (f) $t=0.68$ s.  (a) (b)  Figure 11: Example 2. The
representation of (a) the magnetic vector potential $A_{z}$, and (b) the
magnitude of the magnetic field ${\bm{\mathnormal{B}}}$ at $t=0.68$ s (final
failure state).

###  4.3 .Example 3: Magneto-mechanically induced cracking in a ferromagnetic
material containing predefined notches and wires

This example comprises four sub-examples in which the magneto-mechanically
driven fracture in a ferromagnetic material is examined. In the first two
examples, the electric current is enforced through the notches, and in the
last two examples, the electric current flows through the copper wires. The
cracking response of the ferromagnetic material is examined by solving the
coupled magneto-mechanically driven crack phase-field problem. The following
sub-examples investigate a square plate of ferromagnetic material. The
geometrical configuration of the first two examples with the predefined
notches is shown in Figure 12. We set $A=40$ mm thus
${\mathcal{B}}=(0,80)^{2}$ mm2. The discretization of the domain is performed
with $51200$ triangular elements with an element size of $h\approx 0.707$ mm.
In the last two sub-examples, the square plates contain copper wires with the
radius of $r_{w}=0.3$ mm, see Figure 13. The third sub-example with four
copper wires is discretized with $31132$ triangular elements. In the domain,
the uniform element size is $h\approx 1.02$ mm. The fourth sub-example with
five copper wires is discretized with $30984$ triangular elements by the
uniform element size as $h\approx 1.02$ mm. In these BVPs, all displacements
are fixed in both directions, and the electric vector potential is set to zero
at the boundaries $\partial_{D}\cal{B}$. The notches and copper wires carry
the constant electric source current of $\bar{J}=0.03$ A/mm2.

####  4.3.1 .Sub-example 1: Transversely wired plate with three winding wires.

In this example, the domain contains three predefined notches
${\mathcal{C}}_{1}$, ${\mathcal{C}}_{2}$ and ${\mathcal{C}}_{3}$ of length
$n=8$. The notch ${\mathcal{C}}_{1}$ is located horizontally in the center of
the domain. The other vertical notches (${\mathcal{C}}_{2}$ and
${\mathcal{C}}_{3}$) are located with distances of $m=20$ mm from the right
and left boundaries, see Figure 12(a). The results in terms of the crack
phase-field, the magnetic vector potential $A_{z}$, and von Mises stress are
exhibited in Figure 14. Moreover, the magnetic field ${\bm{\mathnormal{B}}}$
developed in the plate is shown in Figure 18(a). The variation of the magnetic
vector potential average value with respect to time it includes the effect of
damage is represented in Figure 20.a.

####  4.3.2 .Sub-example 2: Transversely wired plate with nine winding wires.

In this example, the square plate considered in the BVP contains nine
predefined notches ${\mathcal{C}}_{i}$ with $=1,2,3,...,9$. The dimensions of
the notches are the same as the previous sub-example, see Figure 12(b). The
results in terms of the crack phase-field, magnetic vector potential, and von
Mises stress are presented in Figure 15 at different time steps. The magnetic
field arising in the plate under the constant electric current source is shown
at the last time step in Figure 18.b. The average value of the magnetic vector
potential over time in the square plate regarding the cracking effect is
depicted in Figure 20(b).

(a) (b)  Figure 12: Transversely wired plates. The representation of geometry,
the corresponding dimensions and the boundary conditions of the square plate
contains (a) three predefined notches (Example 3.1), and (b) nine predefined
notches (Example 3.2)

In both sub-examples, it is observed that cracks start to evolve in the square
plate when a constant electric current is enforced in the predefined notches.
As the crack propagates over time, the electric vector potential increases
along the crack. However, its value decreases when it does not reach the
cracking state. The evolution of the average value of the magnetic vector
potential over time confirms that when the plate is completely damaged, the
value decreases. The comparison of the results indicates that the average
value of the magnetic vector potential in the plate with nine notches reaches
its maximum value earlier than in the plate with three notches.

####  4.3.3 .Sub-example 3: Longitudinally wired plate with four winding
wires.

In the third case, the square plate contains four copper wires,
${\mathcal{C}}_{i}$ with $=1,2,3,4$. The radius of the wires is
$r_{\text{wire}}=3$ mm. The dimensions of the plate are the same as the
previous sub-examples, see Figure 12(a). The wires are placed of a distance of
$m=20$ mm from the edge of the plate. The crack initiation and propagation,
the evolution of the magnetic vector potential, and the von Mises stress for
the different time steps are presented in Figure 15. The magnetic field under
the constant electric current is presented in Figure 18(a) for the last time
step. The average value of the magnetic vector potential over time in the
square plate regarding the cracking effect is provided in Figure 21.a.

(a) (b)  Figure 13: Longitudinally wired plates. The representation of
geometry and the corresponding dimensions of the square plate contains (a)
four copper wires (Example 3.3), and (b) five copper wires (Example 3.4).
Figure 14: Example 3.1. Transversely wired plate with three winding wires.
Reference results of the magnetostrictively induced crack driven by constant
electric current source induction along the predefined notches. Evolution of
the crack phase-field, magnetic vector potential and von Mises stress at
different time steps.  Figure 15: Example 3.2. Transversely wired plate with
nine winding wires. Reference results of the magnetostrictively induced crack
driven by constant electric current source induction along the predefined
notches. Evolution of the crack phase-field, magnetic vector potential and von
Mises stress at different time steps.  Figure 16: Example 3.3. Longitudinally
wired plate with four winding wires. Reference results of the
magnetostrictively induced crack driven by constant electric current source
induction on the copper wires. Evolution of the crack phase-field, magnetic
vector potential and von Mises stress at different time steps.

####  4.3.4 .Sub-example 4: Longitudinally wired plate with five winding
wires.

In the last example, the square plate contains five copper wires,
${\mathcal{C}}_{i}$ with $=1,2,3,4,5$. The geometries of the plate and wires
and their dimensions are as in the previous example. In addition to Sub-
example 3, this one has an additional copper wire located at the center of the
plate with $n=20$ mm, see Figure 13(b). The results in terms of the crack
phase-field, magnetic vector potential, and von Mises stress are depicted in
Figure 17. The additional findings, including the magnetic field at the last
time step and the average value of the magnetic vector potential, are shown in
Figure 18(b) and Figure 21(b), respectively.

Figure 17: Example 3.4. Longitudinally wired plate with five winding wires.
Reference results of the magnetostrictively induced crack driven by constant
electric current source induction on the copper wires. Evolution of the crack
phase-field, magnetic vector potential and von Mises stress at different time
steps.  (a) (b)  Figure 18: Transversely wired plates. The representation of
magnetic field induced by the constant electric current source along the
predefined notches on (a) square plate with three notches (Example 3.1), and
(b) square plates with nine notches (Example 3.2).  (a) (b)  Figure 19:
Longitudinally wired plates. The representation of magnetic field induced by
the constant electric current source along the copper wires on square plate
with (a) four copper wires (Example 3.3), and (b) five copper wires (Example
3.4).  (a) (b)  Figure 20: Longitudinally wired plates. The average value of
magnetic vector potential $A_{\text{ave}}$ over time for the square plate
contains (a) three notches (Example 3.1), and (b) nine notches (Example 3.2).
(a) (b)  Figure 21: Transversely wired plates. The average value of magnetic
vector potential $A_{\text{ave}}$ over time for the square plate contains (a)
four copper wires (Example 3.3), and (b) five copper wires (Example 3.4).

The results obtained for the last two examples show that induction of a
constant electric current in the copper wires causes mechanical deformation
and corresponding stresses in the plate. When the maximum principal stress
produced under the magnetostrictive effects exceeds the critical stress value
of a ferromagnetic material, the plate starts to degrade. The evolution of the
magnetic vector potential confirms that when the material reaches the cracking
state, its average value decreases. Furthermore, the average value of the
magnetic vector potential in the plate with five wires reaches its maximum
value earlier than in the plate with four wires.

At this point it is necessary to remark, a further extension to three-
dimensional setting is computationally time demanding. This is mainly due to
the non-convexity of the energy functional to be minimized with respect to the
displacement and the phase-field and non-linearity which calls for the
resolution of small length scales. As such, an idea of using adaptive mesh
refinement [106, 37], or alternatively using global/local approach [38] seems
particularly appealing.

##  5 .Conclusion

A coupled electro-magneto-mechanical model, along with the phase-field
approach is developed for simulating crack growth in ferromagnetic material.
The proposed coupled electro-magneto-mechanical model evaluates the evolution
of the stress response and deformation of the material under the
magnetostrictive effects. The magneto-mechanically-driven cracking is examined
by applying the phase-field approach. To this end, we extended Maxwell’s
equation to a variational-based electro-magneto-elastic model, then coupled it
with the phase-field approach. In this model, the total stress is additively
decoupled into three parts, containing purely mechanical, purely-magnetic, and
the coupled magneto-mechanical stresses. The magneto-driven deformation in a
solid body is computed as a function of the magnetic field, which is
implicitly defined through the magnetic vector potential.

The cracking response of the ferromagnetic material is computationally
formulated on the basis of the degradation of the effective stress tensor. The
degradation function is formulated as a function of the damage variable. Its
initiation and propagation is developed through the thermodynamically
consistent variational form of the crack phase field. The total energy
functional includes the constitutive energy density functions that correspond
to elastic, magnetostrictive, magnetization, and fracture contributions. In
the phase-field method, the crack driving force is determined in terms of the
mechanical free energy function. In the current work, the capability of the
proposed model in predicting the response of the ferromagnetic material is
validated by a couple of representative numerical examples. At first, we used
Maxwell’s equations to investigate the variation of the electric and magnetic
fields in a solid body surrounded by a vacuum. Then, we extended the coupled
magneto-elastic model along with the crack phase-field method to examine the
cracking response under the mechanical deformation of a ferromagnetic material
induced by magnetization.

Several topics for further research emerge from the present study. First, the
assumption of monotonic increasing for specific current density
${\bm{\mathnormal{J}}}_{s}$ could be extended to time-dependent cyclic
magnetization. So, the possibility of fatigue failure in the ferromagnetic
materials with magnetostrictive effects could be elaborated. For this purpose,
the degradation of Griffith’s energy release rate can be considered. Moreover,
thermal effects, which can contribute as an additional source of failure in
the ferromagnetic material, can be taken into consideration. In that case, the
interaction between magnetic and thermal fields needs to be considered (e.g.,
in an electromagnetic rail launcher). In addition, the assumption of
magnetostatics could be further relaxed (requiring the use of Nédélec
elements), such that the variation of the electric displacement in time does
not vanish anymore and hence causes non-zero permittivity in the electric
equation.

### Acknowledgment

N. Noii founded by the Priority Program DFG-SPP 2020 within its second funding
phase. P. Wriggers were funded by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) under Germany’s Excellence Strategy within the
Cluster of Excellence PhoenixD, EXC 2122 (project number: 390833453).

### Appendix A. Approximation of the five material constants in ${W}_{mag}$

The objective is to derive the five material constants denoted as
$g_{i}({\bm{\mathnormal{x}}})$ in ${W}_{mag}$ of (63). Recall the total form
of the Helmholtz free energy function for the description of the
electromagneto-mechanical response of a ferromagnetic material is introduced
by

$W({\Bgothic
C}):={W}_{elas}({\bm{\varepsilon}},d)+{W}_{mag}({\bm{\varepsilon}},{\bm{\mathnormal{B}}};d)+{W}_{mos}({\bm{\varepsilon}},{\bm{\mathnormal{B}}};d).$
(A.1)

and bases on five invariants

$W({\Bgothic
C}):={W}_{elas}(I_{1},I_{2},d)+{W}_{mag}(I_{1},I_{4};d)+{W}_{mos}(I_{5},I_{6};d).$
(A.2)

This expression is based on invariants which express the mechanics, the
magnetization and the magnetostriction of the material. The first two of these
invariants describe the isotropic characteristics of the solid material as a
function of the total strain tensor as follows

$I_{1}({\bm{\varepsilon}})=\mathop{\operator@font
tr}({\bm{\varepsilon}})\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad
I_{2}({\bm{\varepsilon}})=\frac{1}{2}({\bm{\varepsilon}}:{\bm{\varepsilon}}).$
(A.3)

The remaining invariants $I_{4},I_{5}$ and $I_{6}$ describe the single-valued
magnetization and the magnetostrictive curves as a function of the magnetic
flux density ${\bm{\mathnormal{B}}}$. They are provided in the following
expressions

$I_{4}({\bm{\mathnormal{B}}})={\bm{\mathnormal{B}}}\cdot{\bm{\mathnormal{B}}}\,,\quad
I_{5}({\bm{\varepsilon}},{\bm{\mathnormal{B}}})={\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}\cdot{\bm{\mathnormal{B}}}\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad
I_{6}({\bm{\varepsilon}},{\bm{\mathnormal{B}}})={\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}^{2}\cdot{\bm{\mathnormal{B}}}.$
(A.4)

The Cauchy stress tensor is derived by taking the partial derivative of the
total Helmholtz free energy function with respect to the total strain tensor

${\bm{\sigma}}({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d):=\frac{{\partial
W({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)}}{\partial{\bm{\varepsilon}}}=\displaystyle\sum_{i}\frac{{\partial
W({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)}}{\partial I_{i}}\frac{{\partial
I_{i}}}{\partial{\bm{\varepsilon}}}.$ (A.5)

Also, the magnitization, see (65), is determined by taking the partial
derivative of the total Helmholtz free energy function with respect to the
magnetic flux density

${{\bm{\mathnormal{M}}}}({\bm{\varepsilon}},{\bm{\mathnormal{B}}};d):=-\frac{{\partial
W({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)}}{\partial{\bm{\mathnormal{B}}}}=-\displaystyle\sum_{i}\frac{{\partial
W({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)}}{\partial I_{i}}\frac{{\partial
I_{i}}}{\partial{\bm{\mathnormal{B}}}}.$ (A.6)

In (A.5) and (A.6), the partial derivatives of the basic invariants with
respect to the strain tensor and the magnetic flux density vector are defined
as

$\displaystyle\frac{\partial{I_{1}}}{{\bm{\varepsilon}}}=\text{\bm{\mathnormal{I}}}\,\quad\frac{\partial{I_{2}}}{{\bm{\varepsilon}}}={\bm{\varepsilon}}\,,\quad\frac{\partial{I_{5}}}{{\bm{\varepsilon}}}={\bm{\varepsilon}}\otimes{\bm{\varepsilon}},$
(A.7)
$\displaystyle\frac{\partial{I_{6}}}{{\bm{\varepsilon}}}={\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}+{\bm{\varepsilon}}\cdot{\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}},$
$\displaystyle\frac{\partial{I_{4}}}{{\bm{\mathnormal{B}}}}=2{\bm{\mathnormal{B}}}\,,\quad\frac{\partial{I_{5}}}{{\bm{\mathnormal{B}}}}={\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}\quad\mathchoice{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}{\hbox{{and}}}\quad\frac{\partial{I_{6}}}{{\bm{\mathnormal{B}}}}={\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}^{2}.$

Therefore, the Cauchy stress tensor in follows with

$\displaystyle{{\bm{\sigma}}}({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)$
$\displaystyle:=\frac{{\partial
W({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)}}{\partial{\bm{\varepsilon}}}$
(A.8)
$\displaystyle=g_{e}(d){\widetilde{{\bm{\sigma}}}}^{+}({\bm{\varepsilon}})+{\widetilde{{\bm{\sigma}}}}^{-}({\bm{\varepsilon}})+\frac{1}{2}\sum_{i=0}^{4}\frac{1}{i+1}\partial_{I_{1}}g_{i}\frac{I_{4}^{i+1}}{B^{2i}_{\text{ref}}}\textbf{I}$
$\displaystyle+\frac{g(d)}{2}\alpha_{5}{\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}+\frac{g(d)}{2}\alpha_{6}\left({\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}+{\bm{\varepsilon}}\cdot{\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}\right),$

see (72) and (73), as well. The constitutive equation for pure mechanics state
is additively split to purely tensile contribution
${\widetilde{{\bm{\sigma}}}}^{+}({\bm{\varepsilon}})$ and purely compression
contribution ${\widetilde{{\bm{\sigma}}}}^{-}({\bm{\varepsilon}})$, reads

${\bm{\sigma}}(\bm{\varepsilon},s):=\frac{\partial{W}_{elas}({\bm{\varepsilon}},d)}{\partial{\bm{\varepsilon}}}=g_{e}(d)\frac{\partial{\psi}^{+}}{\partial\bm{\varepsilon}}+\frac{\partial{\psi}^{-}}{\partial\bm{\varepsilon}}=g_{e}(d){\bm{\widetilde{\sigma}}^{+}}+{\bm{\widetilde{\sigma}}^{-}},$
(A.9)

with

${\widetilde{{\bm{\sigma}}}}^{+}({\bm{\varepsilon}})=K_{n}H{{}^{+}}(I_{1})(\bm{\varepsilon}:\textbf{I})\textbf{I}+2\mu{\bm{\varepsilon}}^{dev},\quad\text{and}\quad{\widetilde{{\bm{\sigma}}}^{-}}({\bm{\varepsilon}})=K_{n}\big{(}1-H{{}^{+}}(I_{1})\big{)}(\bm{\varepsilon}:\textbf{I})\textbf{I}.$
(A.10)

Correspondingly magnetization vector denoted as ${\bm{\mathnormal{M}}}$
yields:

$\displaystyle{{\bm{\mathnormal{M}}}}({\bm{\varepsilon}},{\bm{\mathnormal{B}}};d)$
$\displaystyle:=-\frac{{\partial
W({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)}}{\partial{\bm{\mathnormal{B}}}}$
(A.11)
$\displaystyle=-\sum_{i=0}^{4}\frac{g_{i}}{i+1}\partial_{I_{4}}\frac{I_{4}^{i+1}}{B^{2i}_{\text{ref}}}-g(d)\Big{(}\frac{1}{2}\alpha_{5}{\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}-\frac{1}{2}\alpha_{6}{\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}^{2}\Big{)},$

see (65), as well. Thus, the total stress tensor is derived as

$\displaystyle{{\bm{\tau}}}({\bm{\varepsilon}},{\bm{\mathnormal{B}}},d)=$
$\displaystyle
g_{e}(d){\widetilde{{\bm{\sigma}}}}^{+}({\bm{\varepsilon}})+{\widetilde{{\bm{\sigma}}}}^{-}({\bm{\varepsilon}})+\frac{1}{2}\sum_{i=0}^{4}\frac{1}{i+1}\partial_{I_{1}}g_{i}\frac{I_{4}^{i+1}}{B^{2i}_{\text{ref}}}\textbf{I}$
(A.12)
$\displaystyle+\left(E_{M}\mu_{0}^{-1}+\sum_{i=0}^{4}\frac{g_{i}}{i+1}\partial_{I_{4}}\frac{I_{4}^{i+1}}{B^{2i}_{\text{ref}}}+\frac{1}{2}\alpha_{5}\right){\bm{\mathnormal{B}}}\otimes{\bm{\mathnormal{B}}}$
$\displaystyle-\left(\frac{E_{M}}{2}\mu_{0}^{-1}+\sum_{i=0}^{4}\frac{g_{i}}{i+1}\partial_{I_{4}}\frac{I_{4}^{i+1}}{B^{2i}_{\text{ref}}}\right){\bm{\mathnormal{B}}}\cdot{\bm{\mathnormal{B}}}\textbf{I}$
$\displaystyle-\Big{(}\frac{1}{2}\alpha_{5}\left({\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}\cdot{\bm{\mathnormal{B}}}\right)\textbf{I}-\frac{1}{2}\alpha_{6}\left({\bm{\mathnormal{B}}}\cdot{\bm{\varepsilon}}^{2}\cdot{\bm{\mathnormal{B}}}\right)\textbf{I}$