Micromagnetics simulation of ferromagnetic materials using large cells

Abstract

In this paper, numerical techniques for simulation of magnetic domain structures of ferromagnetic materials using micromagnetics with large cells are investigated to save the computation cost. The exchange energy is modified to be applicable to large cells considering the angle between two neighboring cells. The demagnetizing energy is also modified to take account of erroneous magnetic surface charges which appear around the domain walls as a result of discretization. The method is then applied to a simple model and the result shows that an ideal vortex magnetic domain structure and domain wall motion can be represented.

Keywords

Magnetic domain micromagnetics exchange energy demagnetization energy

1. Introduction

Several methods for mathematically modeling the hysteresis phenomenon in ferromagnetic materials such as electrical steel sheets have been proposed [1, 2, 3, 4, 5]. However, such methods often require a large number of parameters obtained from measurements. Moreover, these mathematical modeling methods do not consider the magnetic domain structures and the domain wall motion under an applied magnetic field. Therefore, a physical modeling method which considers the domain structures, such as micromagnetics [6], seems attractive.

However, it is difficult to apply micromagnetics to electrical steel sheets, because in a standard micromagnetic simulation, the model is required to be discretized into small elementary cells for which the cell size must be smaller than the so-called exchange length $l_{\textit{exch}}$ of nanometer scale for iron [4]. The reason behind it is that the magnetization $\bm{M}$ is assumed to be uniform within each cell. The assumption is reasonable for small cells because the exchange interaction between the dipole moments of atoms in a ferromagnetic material, which does not allow sharp variation of $\bm{M}$ in small dimensions, can be considered. However, for cells much larger than $l_{\textit{exch}}$ , the smooth variation of $\bm{M}$ is not guaranteed. Simulation of a model with dimensions 0.001 m $\times$ 0.001 m $\times$ 0.001 m discretized into small cells less than the exchange length requires a huge amount of memory and CPU time, which is far beyond the capacity of current personal computers and perhaps even super computers.

In this paper, the possibility of the simulation with cubic cells of relative large size (5 $\times$ 10 ${}^{-6}$ m) is investigated. The anisotropy and the Zeeman energies are treated the same as in small cell simulation. However, the representation of the exchange and the demagnetizing energies are modified to be applicable to large cells. The proposed method is then applied to a simple model with dimensions (155 $\times$ 10 ${}^{-6}$ m) $\times$ (155 $\times$ 10 ${}^{-6}$ m) $\times$ (35 $\times$ 10 ${}^{-6}$ m) under an external magnetic field, and the obtained domain structure and wall motion are shown.

2. Methods

2.1 Micromagnetics

According to micromagnetics, the distribution of $\bm{M}$ in a crystal of a ferromagnetic material is determined so that the total magnetic free energy $E_{t}$ given by the following equation is minimized [4].

$\displaystyle E_{t}=E_{\textit{exch}}+E_{d}+E_{ani}+E_{ext}$ (1)

where $E_{\textit{exch}}$ , $E_{d}$ , $E_{ani}$ and $E_{ext}$ are the exchange, demagnetizing, anisotropy and Zeeman energy, respectively. The representations of the energies are expressed as follows:

$\displaystyle E_{\textit{exch}}=\int\limits_{V}{A(\left|{\nabla m_{x}}\right|^% {2}+\left|{\nabla m_{y}}\right|^{2}+\left|{\nabla m_{z}}\right|^{2})}dv$ (2) $\displaystyle E_{d}=\int\limits_{V}{-\frac{1}{2}\mu_{0}}M_{s}\bm{m}\cdot\bm{H}% _{d}dv$ (3) $\displaystyle E_{ani}=\int\limits_{V}{\left[{K_{1}(m_{x}^{2}m_{y}^{2}+m_{x}^{2% }m_{z}^{2}+m_{y}^{2}m_{z}^{2})+K_{2}m_{x}^{2}m_{y}^{2}m_{z}^{2}}\right]}dv$ (4) $\displaystyle E_{ext}=\int\limits_{V}{-\mu_{0}M_{s}\bm{m}\cdot\bm{H}_{ext}}dv$ (5)

where $A$ is the exchange stiffness constant, $\bm{m}=\bm{M}/M_{s}$ is the direction vector of magnetization, $M_{s}$ is the saturation magnetization, $m_{x}$ , $m_{y}$ and $m_{z}$ are the Cartesian components of $\bm{m}$ , $\bm{H}_{d}$ is the demagnetizing field arising from the internal dipole-dipole interaction, $K_{1}$ and $K_{2}$ are the first and second anisotropy constants, $\bm{H}_{ext}$ is the external applied field and $\mu$ ${}_{0}$ is the permeability of free space and $V$ is the volume of the body [7].

When $\bm{M}$ lies parallel everywhere in the body with the effective field $\bm{H}_{\textit{eff}}$ as defined below, $E_{t}$ becomes the minimum.

$\bm{H}_{\textit{eff}}=-\frac{\partial E_{t}}{\mu_{0}\partial\bm{M}}=\bm{H}_{% \textit{exch}}+\bm{H}_{d}+\bm{H}_{\textit{ani}}+\bm{H}_{\textit{ext}}$ (6)

By neglecting the terms of orders greater than fourth (i.e. $K_{2}=$ 0), $\bm{H}_{\textit{ani}}$ is:

$\bm{H}_{ani}=-\frac{2K_{1}}{\mu_{0}M_{s}}\left[{m_{x}(1-m_{x}^{2}){\rm{\bf i}}% +m_{y}(1-m_{y}^{2}){\rm{\bf j}}+m_{z}(1-m_{z}^{2}){\rm{\bf k}}}\right]$ (7)

where i, j and k are unit vectors along $x, y$ and $z$ axes, respectively.

When the model is uniformly discretized into a grid of tightly-packed identical small cubic cells, an appropriate representation of the exchange energy which covers large angle deviation is given as [8]:

$\displaystyle E_{\textit{exch}}=\frac{A}{2}\sum\limits_{i=1}^{N}{\sum\limits_{% j=1}^{nn}{\left({\frac{\theta_{ij}}{d}}\right)^{2}V_{i}}}=\frac{Ad}{2}\sum% \limits_{i=1}^{N}{\sum\limits_{j=1}^{nn}{\theta_{ij}^{2}}}$ (8) $\displaystyle\bm{H}_{\textit{exch}}=-\frac{1}{\mu_{0}d^{3}M_{s}}\frac{\partial E% _{\textit{exch}}}{\partial\bm{m}_{i}}=\frac{2A}{\mu_{0}d^{2}M_{s}}\sum\limits_% {j=1}^{nn}{\frac{\theta_{ij}}{\sin\theta_{ij}}\bm{m}_{j}}$ (9)

where $\theta_{ij}$ is the angle between $\bm{m}_{i}$ and $\bm{m}_{j}$ ( $\bm{m}_{i}$ . $\bm{m}_{j}=\cos\theta_{ij}$ ), $d$ is distance between two neighboring cells, $N$ is the total number of the cells, $V_{i}$ is the volume of cell $i$ , and nn means that the sum extends over the nearest neighboring cells.

As for $\bm{H}_{d}$ , from the formulation of $E_{d}$ in [9], the $x$ component of the volume-averaged demagnetizing field becomes:

$\displaystyle H_{d,x}(i,j,k)=\sum\limits_{i^{\prime}=1}^{I+1}{\sum\limits_{j^{% \prime}=1}^{J}{\sum\limits_{k^{\prime}=1}^{K}{DW_{xx}(i-i^{\prime},j-j^{\prime% },k-k^{\prime})}}}\sigma_{x}(i^{\prime},j^{\prime},k^{\prime})+\sum\limits_{i^% {\prime}=1}^{I}{\sum\limits_{j^{\prime}=1}^{J+1}{\sum\limits_{k^{\prime}=1}^{K% }{DW_{xy}(i-i^{\prime},j-j^{\prime},k-k^{\prime})}}}\sigma_{y}(i^{\prime},j^{% \prime},k^{\prime})+\sum\limits_{i^{\prime}=1}^{I}{\sum\limits_{j^{\prime}=1}^% {J}{\sum\limits_{k^{\prime}=1}^{K+1}{DW_{xz}(i-i^{\prime},j-j^{\prime},k-k^{% \prime})}}}\sigma_{z}(i^{\prime},j^{\prime},k^{\prime})$ (10)

where $i$ , $j$ and $k$ are the coordinates of the cell centers in the grid of $I\times J\times K$ cells in units of $d$ , and $\sigma_{x}$ , $\sigma_{y}$ and $\sigma_{z}$ are the charge densities of the cell surfaces perpendicular to $x$ , $y$ and $z$ axes respectively given by the following equations:

$\displaystyle\sigma_{x}(i,j,k)=m_{x}(i-1,j,k)-m_{x}(i,j,k)\ \ i=1\ldots I+1,j=% 1\ldots J,k=1\ldots K$ $\displaystyle\sigma_{y}(i,j,k)=m_{y}(i,j-1,k)-m_{y}(i,j,k)\ \ i=1\ldots I,j=1% \ldots J+1,k=1\ldots K$ (11) $\displaystyle\sigma_{z}(i,j,k)=m_{z}(i,j,k-1)-m_{z}(i,j,k)\ \ i=1\ldots I,j=1% \ldots J,k=1\ldots K+1$

And ${DW}_{xx}$ , ${DW}_{xy}$ and ${DW}_{xz}$ are matrices of interaction coefficient differences which can be derived from $W$ matrices in ref. [9]. The equations for $H_{d,y}$ and $H_{d,z}$ are obtained by permutation of coordinates in the Eq. (10).

Figure 1.

The method used in this paper for rotating $\bm{m}$ toward $\bm{h}_{\textit{eff}}$ in order to minimize the total free magnetic energy ( $i$ denotes the cell index).

For energy minimization, an initial distribution is given to $\bm{m}$ and $\bm{H}_{\textit{eff}}$ is calculated using Eqs (7), (9) and (10). For each cell, if the angle between $\bm{m}$ and $\bm{H}_{\textit{eff}}$ is not zero, it is slowly rotated toward $\bm{H}_{\textit{eff}}$ by the method shown in Fig. 1. The relaxation factor, $\beta$ ( $0<\beta<1$ ), determines how fast $\bm{m}$ moves toward $\bm{H}_{\textit{eff}}$ . A small $\beta$ may guarantee the convergence but it often increases the number of iterations required. On the other hand, with a large $\beta$ , convergence may not be achieved. Having tried different values of $\beta$ , we found $\beta=$ 0.5 a good choice. Once $\bm{m}$ rotated for all the cells, $\bm{H}_{\textit{eff}}$ is recalculated and the angle between $\bm{H}_{\textit{eff}}$ and $\bm{m}$ is evaluated for every cell. The procedure is repeated till in all cells the angles fall below a defined tolerance or $E_{t}$ starts increasing due to being close to a shallow minimum.

Figure 2.

Variation of the magnetization angle between two neighboring cells where the rotation takes place within the region shown in gray.

2.2 Modifying the exchange energy representation

For large cells, the exchange energy is modified. The exchange energy calculation of two neighboring cells with small angle difference, which are considered to be within one domain, is treated differently from that with large angle difference, which are considered to be in different domains. For the two neighboring cells of small angle difference within one domain, the angle variation takes place within a depth of approximately 2 $l_{exch}$ in the interface of the two cells, as shown in Fig. 2. The exchange energy is then concentrated in the rotation region and should be modified as follows:

$E_{exch}^{\ast}=\frac{A}{2}\sum\limits_{i=1}^{N}{\sum\limits_{j=1}^{nn}{\left(% {\frac{\theta_{ij}}{2l_{\textit{exch}}}}\right)^{2}2}}l_{\textit{exch}}d^{2}=% \frac{Ad^{2}}{4l_{\textit{exch}}}\sum\limits_{i=1}^{N}{\sum\limits_{j=1}^{nn}{% \theta_{ij}^{2}}}$ (12)

For the two neighboring cells with large angle difference in different domains, the cells boundaries are considered as domain walls. The domain wall energy includes the exchange energy and anisotropic energy and is the result of a balance between them. The actual domain wall energy is two times of the exchange energy or the anisotropic energy as it reaches its minimum when the exchange and the anisotropic energies are equal. Thus the exchange energy for cells with large angle difference can be calculated according to the domain wall energy as follows, which is derived by using domain wall energies with random degree transition in the walls [10].

$\displaystyle E_{\textit{exch}}^{\ast}=2Ad^{2}\theta_{ij}^{2}/\delta$ (13) $\displaystyle\delta=\sqrt{\frac{32A\theta_{ij}^{2}}{K_{1}(4\theta_{ij}^{2}-(% \sin 4(\alpha+\theta_{ij}^{2})-\sin 4\alpha))}}$ (14)

where $\alpha$ is the angel between $\bm{m}_{i}$ and the easy axis.

In summary, if the angle difference $\theta_{ij}$ is smaller than a threshold value, $\theta_{w}$ , the cell and the neighboring one are considered belonging to the same domain and the modified exchange energy representation in Eq. (12) is used. When $\theta_{ij}>\theta_{w}$ , the cells are considered to be in different domains and the true exchange energy representation inside domain wall given by Eq. (13) is used.

Figure 3.

(a) A continuum model; the domain structure has zero demagnetizing energy. (b) The same model when discretized into cubic cells; the erroneous surface charges on the stair-like walls create erroneous demagnetizing field with large magnitudes in the neighboring cells.

2.3 Modifying the demagnetization energy representation

As a consequence of discretization, erroneous surface charges may appear around the domain walls as shown in Fig. 3. The erroneous charges create quite large erroneous demagnetizing field which severely affects the simulated domain structure. The demagnetizing field needs to be modified to take account of such errors. In order to suppress the error, any surface charge $\sigma$ on the walls are ignored. To implement the modification, for any surface except for the bounding surface of the model, the angle difference between the magnetization of both sides of the surface is evaluated. If the angle difference is larger than $\theta_{w}$ (the same criterion as in Section 2.2), the surface charge is neglected otherwise it is taken into account. This is equivalent of the following equations:

$\displaystyle\sigma_{x}^{\ast}(i,j,k)=0\qquad\quad\quad\ \ \textit{if}\ \ \ \ % \ {\Delta\theta_{x}}\geqslant\theta_{w}$

(15) $\displaystyle\sigma_{x}^{\ast}(i,j,k)=\sigma_{x}(i,j,k)\ \ \ \ \textit{if}\ \ % \ \ \ {\Delta\theta_{x}}<\theta_{w}$

where $\Delta\theta_{x}$ denotes the angle difference between $\bm{m}$ ( $i-1$ , $j$ , $k$ ) and $\bm{m}$ ( $i$ , $j$ , $k$ ). Similar treatments are required for $\sigma_{y}^{\ast}(i,j,k)$ and $\sigma_{z}^{\ast}(i,j,k)$ .

Figure 4.
Distributions of magnetic moments $\bm{M}$ s. (a) Initial state, (b) obtained by using the ordinary $E_{\textit{exch}}$ , (c) the first stage, obtained by using $E_{\textit{exch}}$ , $\theta_{w}=$ 180 deg., (d) the second stage, obtained by using $E_{\textit{exch}}$ , $\theta_{w}=$ 10 deg.

Figure 5.
Distributions of magnetic moments $\bm{M}$ s. (a) under applied fields with $H=$ 1000 A/m, (b) under applied fields with $H=$ $-$ 1000 A/m.

3. Results and discussion

The proposed method is applied to a simple model with the dimension of (155 $\times$ 10 ${}^{-6}$ m) $\times$ (155 $\times$ 10 ${}^{-6}$ m) $\times$ (35 $\times$ 10 ${}^{-6}$ m) under a zero applied field first. The ideal domain arrangement for zero resultant magnetic moments is a vortex domain structure [10]. The material parameters are the typical values for iron, i.e., the exchange stiffness constant $A=$ 2.07 $\times$ 10 ${}^{-11}$ J/m, the anisotropy constant $K_{1}=$ 4.8 $\times$ 10 ${}^{4}$ J/m ${}^{3}$ , and the saturation magnetization $M_{s}=$ 2.19 T. In the first stage at the initial state, the magnetic moments are oriented randomly as shown in Fig. 4a and no domain walls exist. Thus the distribution of $\bm{M}$ s is calculated by using modified $E_{\textit{exch}}$ * in Eq. (13) for all neighboring cells, namely $\theta_{w}=$ 180 deg. A rough vortex domain structure is formed as shown in Fig. 4c, in which the domain walls are not so obvious, whereas it is improved compared with the result in Fig. 4b obtained by using the ordinary $E_{\textit{exch}}$ . Then in the second stage, domain structure exists thus $\bm{M}$ s distribution is calculated by using modified $E_{\textit{exch}}$ * in Eq. (12) in neighboring cells with $\theta_{ij}$ less than 10 deg., namely only the exchange energies inside one domain are counted and modified. An ideal vortex domain structure with clear domain walls is formed as shown in Fig. 4d. The proposed method is also applied to the same model under applied magnetic fields with $H=$ 1000 A/m and $H=$ $-$ 1000 A/m. The material parameters are the typical values for iron the same as shown above. The obtained domain structure is shown in Fig. 5. Domain wall motion is seen as that the domain with $\bm{M}$ s in the same direction of the applied field is enlarged and the domain with $\bm{M}$ s in the opposite direction of the applied field becomes smaller.

4. Conclusions

In the micromagnetics simulation of ferromagnetic materials using large cells, the ideal vortex domain structure can be obtained by improving the calculation methods of exchange and demagnetizing energies. And domain wall motion is also realized by applying external fields.

References

Mayergoyz

I.D.

, Dynamic Preisach models of hysteresis, IEEE Transactions on Magnetics24(6) (1988), 2925–2927.

Boukhtache

Azoui

and Féliachi

, Optimized model for magnetic hysteresis in silicon-iron sheets by using the simulated annealing algorithm, International Journal of Applied Electromagnetics and Mechanics30(1–2) (2009), 1–7, DOI: 10.3233/JAE-2009-1012.

Francesco

R.F.

and Alessandro

, Hysteresis model identification by the Flock-of-Starlings optimization, International Journal of Applied Electromagnetics and Mechanics30(3–4) (2009), 321–331, DOI: 10.3233/JAE-2009-1032.

Guérin

Leconte

Dular

Geuzaine

and Mazauric

, Using a vector Jiles-Atherton static model to simulate the effect of hysteresis on the operation of fast-acting linear actuators, International Journal of Applied Electromagnetics and Mechanics52(1–2) (2016), 547–554, DOI: 10.3233/JAE-162078.

Matsuo

, Anisotropic vector hysteresis model using an isotropic vector play model, IEEE Transactions on Magnetics46(8) (2010), 3041–3044.

Brown

W.F.

, Micromagnetics, John Wiley & Sons Inc., 1963.

Kronmueller

and Faehnle

, Micromagnetism and the microstructure of ferromagnetic solids, Cambridge University Press, 2003.

Donahue

M.J.

and McMichael

R.D.

, Exchange energy representation in computational micromagnetics, Physica B233(4) (1997), 272–278.

Fabian

Kirchner

Williams

Heider

Leibl

and Hubert

, Three-dimensional micromagnetic calculations for magnetite using fft, Geophys. J. Int.124 (1996), 89–104.

10.

Kittel

, Physical theory of ferromagnetic domains, Reviews of Modern Physics21(4) (1949), 541–581.