Computational modeling of magnetic hysteresis with thermal effects

# Computational modeling of magnetic hysteresis with thermal effects

## Abstract

We study computational behavior of a mesoscopic model describing temperature/external magnetic field-driven evolution of magnetization. Due to nonconvex anisotropy energy describing magnetic properties of a body, magnetization can develop fast spatial oscillations creating complicated microstructures. These microstructures are encoded in Young measures, their first moments then identify macroscopic magnetization. Our model assumes that changes of magnetization can contribute to dissipation and, consequently, to variations of the body temperature affecting the length of magnetization vectors. In the ferromagnetic state, minima of the anisotropic energy density depend on temperature and they tend to zero as we approach the so-called Curie temperature. This brings the specimen to a paramagnetic state. Such a thermo-magnetic model is fully discretized and tested on two-dimensional examples. Computational results qualitatively agree with experimental observations. The own MATLAB code used in our simulations is available for download.

Keywords: dissipative processes, hysteresis, micromagnetics, numerical solution, Young measures

## 1 Introduction

In the isothermal situation, the configuration of a rigid ferromagnetic body occupying a bounded domain is usually described by a magnetization which denotes density of magnetic spins and which vanishes if the temperature is above the so-called Curie temperature . Brown [7] developed a theory called “micromagnetics” relying on the assumption that equilibrium states of saturated ferromagnets are minima of an energy functional. This variational theory is also capable of predictions of formation of domain microstructures. We refer e.g. to [17] for a survey on the topic. Starting from a microscopic description of the magnetic energy we will continue to a mesoscopic level which is convenient for analysis of magnetic microstructures.

On microscopic level, the magnetic Gibbs energy consists of several contributions, namely an anisotropy energy , where is the-so called anisotropy energy density describing crystallographic properties of the material, an exchange energy penalizing spatial changes of the magnetization, the non-local magnetostatic energy , work done by an external magnetic field which reads , and a calorimetric term . The anisotropic energy density depends on the material properties and defines the so-called easy axes of the material, i.e., lines along which the smallest external field is needed to magnetize fully the specimen. There are three types of anisotropy: uniaxial, triaxial, and cubic. Furthermore, is supposed to be a nonnegative function, even in its first variable, i.e., are assigned the same anisotropic energy. In the magnetostatic energy, is the magnetostatic potential related to by the Poisson problem arising from Maxwell equations. Here denotes the characteristic function of and is the permeability of vacuum.

A widely used model describing steady-state isothermal configurations is due to Landau and Lifshitz [20, 21] (see also e.g. Brown [7] or Hubert and Schäfer [13]), relying on minimization of Gibbs’ energy with as a fixed parameter, i.e.,

 Unknown environment '% (1)

where the anisotropy energy is considered in the form

 ψ(m,θ):=ϕ(m)+a0(θ−θc)|m|2−ψ0(θ), (2)

where determines the intensity of the thermo-magnetic coupling. To see a paramagnetic state above Curie temperature , one should consider . The isothermal part of the anisotropy energy density typically consists of two components , where is chosen in such a way to attain its minimum value (typically zero) precisely on lines , where each , determines an axis of easy magnetization. Typical examples are for uni-axial, for triaxial, and for cubic magnets. We can consider a uniaxial magnet with , for instance. Here, the easy axis coincides with the -th axis of the Cartesian coordinate system, i.e., . On the other hand, is used to ensure that, for , is minimized at for and that is coercive. Such energy has already been used in [27]. For , the exchange energy guarantees that the problem (1) has a solution . Zero-temperature limits of this model consider, in addition, that the minimizers to (1) are constrained to be valued on the sphere with the radius and were investigated, e.g., by Choksi and Kohn [10], DeSimone [11], James and Kinderlehrer [14], James and Müller [15], Pedregal [24, 25], Pedregal and Yan [26] and many others.

In [5], the authors first consider a mesoscopic micromagnetic energy arising for setting in (1). Moreover, it is assumed that changes of magnetization cause dissipation which is transformed into heat. Increasing temperature of the specimen influences its magnetic properties. Therefore, they analyze an evolutionary anisothermal mesoscopic model of a magnetic material. The aim of this paper is to discretize this model in space and time, and to perform numerical experiments. The plan of our work is as follows. In Section 2 we describe the stationary mesoscopic model. The evolutionary problem is introduced in Section 3. Section 4 provides us with a numerical approximation and some computational experiments. We finally conclude with a few remarks in Section 5. An appendix then briefly introduces an important tool for the analysis as well as for numerics, namely Young measures.

## 2 Mesoscopic description of magnetization

For small, minimizers of (1) typically exhibit fast spatial oscillations, usually called microstructure. Indeed, the anisotropy energy, which forces magnetization vectors to be aligned with the easy axis (axes), competes with the magnetostatic energy preferring divergence-free magnetization fields. It was shown in [11] by a scaling argument that for large domains the exchange energy contributions becomes less and less significant in comparison with other terms and thus the so-called ”no-exchange” formulation is a justified approximation. This generically leads, however, to nonexistence of a minimum for uniaxial ferromagnets as shown in [14] without an external field . Hence, various ways to extend the notion of a solution were developed. The idea is to capture the limiting behavior of minimizing sequences of as . This leads to a “relaxed problem” (3) involving possibly so-called Young measures ’s [34] which describe fast spatial changes of the magnetization and can capture limit patterns.

It can be proved [11, 24] that this limit configuration solves the following minimization problem involving temperature as a parameter and a “mesoscopic” Gibbs’ energy :

 Missing dimension or its units for \hskip (3)

where the “momentum” operator “” is defined by and similarly for which denotes the identity and . Here, the set of Young measures can be viewed as a collection of probability measures such that is a probability measure on for almost every . It means that is a positive Radon measure such that . We refer to Appendix for more details on Young measures.

In [5], the authors built and analyzed a mesoscopic model in anisothermal situations. A closely related thermodynamically consistent model on the microscopic level was previously introduced in [27] to model a ferro/para magnetic transition. Another related microscopic model with a prescribed temperature field was investigated in [4]. The goal of this contribution is to discretize the model from [5] and test it on computational examples. In order to make our exposition reasonably self-content, we closely follow the derivation of the model presented in [5]. We also point out that computationally efficient numerical implementation of isothermal models can be found in [8, 16, 18, 19], where such a model was used in the isothermal variant.

In what follows we use a standard notation for Sobolev, Lebesgue spaces and the space of continuous functions. We denote by the space of continuous functions vanishing at infinity. Further, , and .

## 3 Evolution problem and dissipation

If the external magnetic field varies during a time interval with a horizon , the energy of the system and magnetic states evolve, as well. Changes of the magnetization may cause energy dissipation. As the magnetization is the first moment of the Young measure, , we relate the dissipation on the mesoscopic level to temporal variations of some moments of and consider these moments as separate variables. This approach was already used in micromagnetics in [30, 31] and proved to be useful also in modeling of dissipation in shape memory materials, see e.g. [23]. In view of (2), we restrict ourselves to the first two moments defining giving rise to the constraint

 λ=L\thinspace∙ν , % where  L(m):=(m,|m|2) (4)

and consider the specific dissipation potential depending on a “yield set”

 ζ(\LARGE.λ):=δ∗S(\LARGE.λ)+ϵq|\LARGE.λ|q,q≥2, (5)

The set determines activation threshold for the evolution of . It is a convex compact set containing zero in its interior. The function is the Fenchel conjugate of the indicator function of . Consequently, it is convex and degree-1 positively homogeneous with . In fact, the first term describes purely hysteretic losses, which are rate-independent and which we consider dominant, and the second term models rate-dependent dissipation.

In view of (2)–(3), the specific mesoscopic Gibbs free energy, expressed in terms of , and , reads as

 g(t,ν,λ,θ):=ϕ\thinspace∙ν+(θ−θc)→a⋅λ−ψ0(θ)+12m⋅∇um−h(t)⋅m (6a) with  m={\rm id}\thinspace∙ν (6b)

where we denoted with from (2) and, of course, again from (1), which makes non-local.

As done already in [5], we relax the constraint (4) by augmenting the total Gibbs free energy (i.e., integrated over ) by the term with (presumably large) and with . Thus, ’s no longer exactly represent the “macroscopic” momenta of the magnetization but rather are in a position of a phase field or an internal parameter of the model. We define the mesoscopic Gibbs free energy as

 G(t,ν,λ,θ):=∫Ω(g(t,ν,λ,θ)+ϰ2|∇Δ−1(λ−L% \thinspace∙ν)|2)dx (7)

with meaning the inverse of the homogeneous Dirichlet boundary-value problem for the Laplacean defined as a map .

The value of the internal parameter may influence the magnetization of the system and vice versa and, on the other hand, dissipated energy influences the temperature of the system, which, in turn, may affect the internal parameters. In order to capture all these effects, we employ the concept of generalized standard materials [12] known from continuum mechanics and couple our micromagnetic model with the entropy balance with the rate of dissipation on the right-hand side; cf. (9). Then the Young measure is considered to evolve quasistatically according to the minimization principle of the Gibbs energy while the dissipative variable is governed by the flow rule:

 ∂ζ(\LARGE.λ)=∂λg(t,ν,λ,θ) (8)

with denoting the subdifferential of the convex functional and similarly is the subdifferential of the convex functional . In our specific choice, (8) takes the form . Furthermore, we define the specific entropy by the standard Gibbs relation for entropy, i.e. , and write the entropy equation

 θ\LARGE.s+divj=ξ(\LARGE.λ)= heat production rate, (9)

where is the heat flux governed by the Fourier law

 j=−K∇θ (10)

with a heat-conductivity tensor . In view of (5),

 Missing dimension or its units for \hskip (11)

Now, since , it holds . Using also , we may reformulate the entropy equation (9) as the heat equation

 Missing dimension or its units for \hskip (12)

where is the specific heat capacity.

Altogether, we can formulate our problem for unknowns , and which was first set and analyzed in [5] as

 minimize∫Ω(ϕ\thinspace∙ν+(θ−θc)→a⋅λ(t)−ψ0(θ(t))+12m⋅∇um−h(t)⋅m+ϰ2∣∣∇Δ−1(λ(t)−L\thinspace∙ν)∣∣2)dxsubject tom={\rm id}\thinspace∙ν on Ω,div(μ0∇um−χΩm)=0 on Rd,ν∈Yp(Ω;Rd), m∈Lp(Ω;Rd), um∈H1(Rd),⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭for t∈[0,T], (13a) Missing dimension or its units for \hskip (13b) cv(θ)\LARGE.% θ−div(K(λ,θ)∇θ)=δ∗S(\LARGE.λ)+ϵ|\LARGE.λ|q+→a⋅θ\LARGE.λin Q, (13c) (K(λ,θ)∇θ)⋅n+bθ=bθext on Σ:=[0,T]×Γ, (13d)

where we accompanied the heat equation (9) by the Robin-type boundary conditions with denoting the outward unit normal to the boundary , with a phenomenological heat-transfer coefficient, and with an external temperature, both assumed non-negative. Eventually, we equip this system with initial conditions

 λ(0,⋅)=λ0,θ(0,⋅)=θ0  on  Ω, (14)

Transforming (9) by the so-called enthalpy transformation, we obtain a different form of (13) simpler for the analysis. For this, let us introduce a new variable , called enthalpy, by

 w=ˆcv(θ)=∫θ0cv(r)dr. (15)

It is natural to assume positive, hence is, for increasing and thus invertible. Therefore, denote

 Θ(w):={ˆcv−1(w)if w≥00if w<0

and notice that, in the physically relevant case when , . Thus writing the heat flux in terms of gives

 K(λ,θ)∇θ=K(λ,Θ(w))∇Θ(w)=K(λ,w)∇w, % where  K(λ,w):=K(λ,Θ(w))cv(Θ(w)). (16)

Moreover, the terms and obviously do not play any role in the minimization (13a) and can be omitted. Thus we may rewrite (13) in terms of as follows:

 minimize∫Ω(ϕ\thinspace∙ν+12m⋅∇um−h(t)⋅m+ϰ2∣∣∇Δ−1(λ(t)−L\thinspace∙ν)∣∣2)dxsubject tom={\rm id}\thinspace∙ν,   on Ω,div(μ0∇um−χΩm)=0 on Rd,ν∈Yp(Ω;Rd), m∈Lp(Ω;Rd), um∈H1(Rd),⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭for t∈[0,T], (17a) Missing dimension or its units for \hskip (17b) Missing dimension or its units for \hskip (17c) (K(λ,w)∇w)⋅n+bΘ(w)=bθext on Σ. (17d)

Eventually, we complete this transformed system by the initial conditions

 λ(0,⋅)=λ0,w(0,⋅)=w0:=ˆcv(θ0)  on  Ω, (18)

where is the initial phase field value, and is the initial temperature.

Now we are ready to define a weak solution to our problem. We denote by the set of time-dependent Young measures, i.e., the set of maps . We again refer to Appendix for details on Young measures.

###### Definition 3.1 (Weak solution[5]).

The triple

 (ν,λ,w)∈(Yp(Ω;Rd))[0,T]×W1,q([0,T];Lq(Ω;Rd+1))×L1([0,T];W1,1(Ω))

such that and is called a weak solution to (17) if it satisfies:

1. The minimization principle: For all in and all

 G(t,ν,λ,Θ(w))≤G(t,~ν,λ,Θ(w)). (19)
2. The magnetostatic equation: For a.a. and all

 μ0∫Rd∇um⋅∇φdx=∫Ωm⋅∇φdx. (20)
3. The flow rule: For any

 ∫Q ((Θ(w)−θc)→a⋅(φ−%.λ\large.λ)+δ∗S(φ)+ϵq|φ|q+ϰ∇Δ−1(λ−L\thinspace∙ν)⋅∇Δ−1(φ−\LARGE.λ))dxdt Missing dimension or its units for \hskip (21)
4. The enthalpy equation: For any

 ∫Q(K(λ,w)∇w⋅∇φ−w\LARGE.φ)dxdt+∫ΣbΘ(w)φdSdt=∫Ωw0φ(0)dx Missing dimension or its units for \hskip (22)
5. The initial conditions in (18): and .

Data qualifications:
The following the data qualification are needed in [5] to prove the existence of weak solutions; cf. [5]:

 isothermal part of the anisotropy energy: ϕ∈C(Rd) and ∃cA1,cA2>0, p>4:  cA1(1+|⋅|p)≤ϕ(⋅)≤cA2(1+|⋅|p), (23b) dissipation function: δ∗S∈C(Rd+1) positively homogeneous, and ∃c1,D,c2,D>0:  c1,D(|⋅|)≤δ∗S(⋅)≤c2,D(|⋅|), (23c) external magnetic field: h∈C1([0,T];L2(Ω;Rd)), (23d) specific heat capacity: cv∈C(R) and, with q from (5), ∃c1,θ,c2,θ>0, ω1≥ω≥q′, c1,θ(1+θ)ω−1≤cv(θ)≤c2,θ(1+θ)ω1−1, (23e) heat conduction tensor: K∈C(Rd+1×R;Rd×d) and Extra open brace or missing close brace (23f) external temperature: θext∈L1(Σ),  θext≥0,  and  b∈L∞(Σ),  b≥0, (23g) initial conditions: ν0∈Yp(Ω;Rd) solving % (???) ,λ0∈Lq(Ω;Rd+1),w0=ˆcv(θ0)∈L1(Ω) with θ0≥0. (23h)

The following theorem is proved in [5].

###### Theorem 3.1.

Let (23) hold. Then at least one weak solution to the problem (17) in accord with Definition 3.1 does exist. Moreover, some of these solutions satisfies also

 w∈Lr([0,T];W1,r(Ω))∩W1,1(I;W1,∞(Ω)∗) with 1≤r

The proof of the Theorem 3.1 in [5] exploits the following time-discrete approximations which also create basis for our fully discrete solution. Given and we call the triple the discrete weak solution of (17) subject to boundary condition (17d) at time-level , , if it satisfies:

In (25d), we denoted by and respectively suitable approximation of the original initial conditions and such that

 λ0,τ→λ0  strongly in Lq(Ω;Rd+1), and  ∥λ0,τ∥L2q(Ω;Rd+1)≤Cτ−1/(2q+1), (26a) w0,τ→w0  strongly in L1(Ω), and  w0,τ∈L2(Ω). (26b)

Moreover and are defined in such a way that their piecewise constant interpolants

 [¯θext,τ,¯bτ](t):=(θkext,τ,bkτ,) for  (k−1)τ

satisfy

 ¯θext,τ→θext  strongly in L1(Σ) and  ¯bτ\lx@stackrel∗⇀b  weakly* in L∞(Σ). (27)

We introduce the notion of piecewise affine interpolants and defined by

 [λτ,wτ](t):=t−(k−1)ττ(λkτ,wkτ)+kτ−tτ(λk−1τ,wk−1τ) for  t∈[(k−1)τ,kτ]

with . In addition, we define the backward piecewise constant interpolants , , and by

 Missing or unrecognized delimiter for \big (28)

Finally, we also need the piecewise constant interpolants of delayed enthalpy and magnetization , defined by

 [w––τ(t),m––τ(t)]:=[wk−1τ,{\rm id}\thinspace∙νk−1τ] for  (k−1)τ

### 3.1 Energetics

In this section we summarize some basic energetic estimates available for our model. First we define the purely magnetic part of the Gibbs free energy as

 Extra open brace or missing close brace (30)

The purely magnetic part of the Gibbs energy satisfies (see [5, Formula (4.19)]) the following energy inequality

 Missing dimension or its units for \hskip (31)

with .

As is a minimizer of (25a), the partial sub-differential of the cost functional with respect to has to be zero at . This condition holds at each time level and, thus, summing up for gives

 ∫tℓ0∫Ω(δ∗S(\LARGE.λτ)+ϵq|\LARGE.λτ|q)dxdt≤∫tℓ0(ϰ⟨⟨¯λτ−L\thinspace∙¯ντ,vτ−\LARGE.λτ⟩⟩ Missing dimension or its units for \hskip (32)

where is an arbitrary test function such that is piecewise constant on the intervals and for every .

Hence, for we get the energy balance of the thermal part of the Gibbs energy, namely

 ∫tℓ0∫Ω(δ∗S(\LARGE.λτ)+ϵq|\LARGE.λτ|q)dxdt≤∫tℓ0(−ϰ⟨⟨¯λτ−L\thinspace∙¯ντ,\LARGE.λτ⟩⟩−∫Ω(Θ(w––τ)−θc)→a⋅\LARGE.λτ+2qτ|¯λτ|2q−2¯λτ\LARGE.λτ)dt . (33)

This inequality couples the dissipated energy and temperature evolution.

## 4 Numerical approximations and computational examples

Dealing with a numerical solution, we have to find suitable spatial approximations for , , , and in each time step. In our numerical method, we require that (4) is satisfied which means that knowing the Young measure we can easily calculate the momenta . We present a spatial discretization of involved quantities in each time step.

The domain of the ferromagnetic body is discretized by a regular triangulation in triangles (in 2D) or in tetrahedra (in 3D) for which will be called elements. The triangulations are nested, i.e., that , so that the discretizations are finer as increases. Let us now describe the approximation.

Young measure. Young measures are parametrized (by probability measures supported on . Hence, we need to handle their discretization in as well as in . Our aim is to approximate a general Young measure by a convex combination of a finite number of Dirac measures (atoms) supported on such that this convex combination is elementwise constant. Let us now describe a rigorous procedure how to achieve this goal. We first omit the time discretization parameter and discuss the discretization of the Young measure in . In order to approximate a Young measure , we follow [9, 22] and define for the following projection operator ( denotes the -dimensional Lebesgue measure)

 [Π1ℓz](x,s)=1Ld(△)∫△z(~x,s)d~x  if x∈△∈Tℓ .

Notice that is elementwise constant in the -variable. We now turn to a discretization of in terms of large cubes in , i.e., for we consider a cube (i.e. we call it “a cube” even if ) which is discretized into smaller cubes with the edge length for some . Corners of small cubes are called nodal points. We define elements on the cube which consist of tensorial products of affine functions in each spatial variable of . In this way, we find basis functions for such that and for all . Moreover, if is the -th nodal point then , where is the Kronecker symbol. Further, each can be continuously extended to and such an extended function can even vanish at infinity, i.e., it belongs to . This construction defines a projector as

 [Π2α,nz](x,s):=(n+1)d∑i=1z(x,si)fi(s) .

Finally, we define , so that

 [Πℓ,α,nz](x,s):=1Ln(△)(n+1)d∑i=1∫△z(~x,si)vi(s)d~x %ifx∈△∈Tℓ .

If we now take and denote we calculate

 ∫Ω∫Rd[Πlz](x,s)νx(ds)dx=∫Ω∫Rdz(x,s)[νl]x(ds)dx , (34)

where for

 [νl]x:=(n+1)d∑i=1ξi,l(x)δsi , (35)

with

 ξi,l(x):=1Ld(△)∫△∫Rdfi(s)νx(ds)dx , x∈△∈Tℓ .

Let us denote the subset of Young measures from which are in the form of (35) by . Notice that and that . Hence, the projector corresponds to approximation of by a spatially piecewise constant Young measure which can be written as a convex combination of Dirac measures (atoms). We refer to [29] for a thorough description of various kinds of Young measure approximations. In order to indicate that the measure is time-dependent we write in the -th time-step

 [νkl,τ]x:=(n+1)d∑i=1ξki,l,τ(x)δsi .

Magnetostatic potential. Following [8], we simplify the calculation of the reduced Maxwell system in magnetostatics by assuming that the magnetostatic potential vanishes outside a large bounded domain . Hence, given , we solve the Poisson problem on with homogeneous Dirichlet boundary condition on . The set is discretized by an outer triangulation that contains the triangulation of the ferromagnetic magnetic body. Then, the magnetostatic potential

 umkl,τ∈P10(^Tℓ) (36)

in the -th time-step is approximated in the space of scalar nodal and elementwise linear functions defined on the triangulation and satisfying zero Dirichlet boundary conditions on the triangulation boundary . For illustration, see Figure 1. The magnetization vector

 mkl,τ∈P0(Tℓ)d (37)

in the -th time-step is approximated in the space of vector and elementwise constant functions. Another numerical approaches to solutions of magnetostatics using e.g. BEM are also available [3].

Enthalpy. The enthalpy

 wkℓ,τ∈P1(Tℓ) (38)

in the -th time-step is approximated in the space of scalar nodal and elementwise linear functions.

Having time and spatial discretizations we can set up an algorithm to solve the problem which is just (25) with additional spatial discretization. Finally, we apply the spatial discretization just described and we arrive at the following problem.

Given spatially discretized boundary condition (17d) and we solve:

where is calculated via (34) and is a piecewise affine approximation of . There is no initial condition for as it is now fully determined by .

In computations, several simplifications were taken to account. First of all, we assume

 d=2,q=2. (40)

In view of (4), the macroscopic magnetization is elementwise constant and it is the first moment of . As the anisotropy energy density is minimized for a given temperature on a sphere in we put the support of the Young measure on this sphere and its vicinity to decrease the number of variables in our problem. In what follows, the number of Dirac atoms in is denoted by . It is then convenient to work in polar coordinates where is the radius and the corresponding angle of the -th atom. Hence, we have

 mkl,τ=λk1,l,τ=pkτN∑i=1ξki,l,τri(cos(φi),sin(φi)),λk2,l,τ=(pkτ)2N∑i=1ξ