next up previous contents
Next: 3.5 Effective Field Up: 3. Finite Element Micromagnetics Previous: 3.3.3 Zeeman Energy   Contents


3.4 Demagnetizing Field and Magnetostatic Energy

The demagnetizing field is a little more complicated to handle, because it is an ``open boundary problem'' with one of its boundary conditions at infinity. In order to overcome this problem Fredkin and Koehler [41,14,42] proposed a hybrid finite element/boundary element method, which requires no finite elements outside the magnetic domain $\Omega$.

Since we assume no free currents in our system, we can calculate the demagnetizing field using a magnetic scalar potential $\varphi (\boldsymbol{x})$. It has to satisfy

$\displaystyle \Delta \varphi$ $\textstyle =$ $\displaystyle \nabla \cdot \boldsymbol{J}(\boldsymbol{x}) \quad \mathrm{for } \boldsymbol{x} \in \Omega$ (3.34)
$\displaystyle \Delta \varphi$ $\textstyle =$ $\displaystyle 0 \quad \mathrm{for } \boldsymbol{x} \not\in \Omega$ (3.35)

with the boundary conditions at the boundary $\Gamma$ of $\Omega$
\begin{displaymath}
\mathrm{Div }\varphi = 0
\end{displaymath} (3.36)

and
\begin{displaymath}
\mathrm{Div }\frac{\partial \varphi }{\partial \boldsymbol{n}} = -\boldsymbol{n} \cdot \boldsymbol{J} \quad.
\end{displaymath} (3.37)

In addition it is required that $\varphi \to 0$ for $\vert\boldsymbol{x}\vert \to \infty$. The weak formulation of $\nabla \cdot \boldsymbol{J}$ is simply given by
\begin{displaymath}
\int_\Omega \nabla \cdot \boldsymbol{J}  d{v}  =
\int_\O...
...ega \sum_j \sum_{k}^{\{x,y,z\}} u_{j,k} \nabla_k \eta_j \quad,
\end{displaymath} (3.38)

which can again be written in matrix-vector format as
\begin{displaymath}
\boldsymbol{d} = D \cdot \boldsymbol{u}
\end{displaymath} (3.39)

with
\begin{displaymath}
D_{j,3j+k}=\int_\Omega u_{j,k} \nabla_k \eta_j \quad,
\end{displaymath} (3.40)

where $k$ stands for the three Cartesian components $\{x,y,z\}$.

The main idea now is to split the magnetic scalar potential $\varphi $ into $\varphi _1$ and $\varphi _2$. Then the problem can be reformulated for these potentials as

\begin{displaymath}
\Delta \varphi _1 = \nabla \cdot \boldsymbol{J}
\end{displaymath} (3.41)

with the boundary condition
\begin{displaymath}
\frac{\partial \varphi _1}{\partial \boldsymbol{n}}= \boldsymbol{n} \cdot \boldsymbol{J} \quad.
\end{displaymath} (3.42)

In addition $\varphi _1(\boldsymbol{x})=0$ for $\boldsymbol{x} \not\in \Omega$.

As a result, we find for $\varphi _2$

\begin{displaymath}
\Delta \varphi _2=0
\end{displaymath} (3.43)

with
\begin{displaymath}
\mathrm{Div }\varphi _2 = \varphi _1
\end{displaymath} (3.44)

and
\begin{displaymath}
\mathrm{Div }\frac{\partial \varphi _2}{\partial \boldsymbol{n}} = 0 \quad.
\end{displaymath} (3.45)

It is required that $\varphi _2 \to 0$ for $\vert\boldsymbol{x}\vert \to \infty$.

Potential theory tells us that

\begin{displaymath}
\varphi _2(\boldsymbol{x})=\frac{1}{4\pi}\int_\Gamma \varph...
...{y})}{\partial \boldsymbol{n}(\boldsymbol{y})}  d{a}  \quad,
\end{displaymath} (3.46)

where $G(\boldsymbol{x}, \boldsymbol{y})=1/\vert\boldsymbol{x} - \boldsymbol{y}\vert$ is the Green function.

$\varphi _1$ can be easily calculated using the standard finite element method as explained in Sec. 2.

The (numerically expensive) evaluation of Eq. (3.46) in all $\Omega$ can be avoided by just calculating the boundary values of $\varphi _2$ on $\Gamma$ and then solving the Dirichlet problem Eq. (3.43) with the given boundary values. For $\boldsymbol{x} \to \Gamma$ Eq. (3.46) is given by

\begin{displaymath}
\varphi _2(\boldsymbol{x})=\frac{1}{4\pi}
\int_\Gamma \var...
...mbol{x})}{4 \pi}-1
\right)
\varphi _1(\boldsymbol{x}) \quad,
\end{displaymath} (3.47)

where $S(\boldsymbol{x})$ denotes the solid angle subtended by $\Gamma$ at $\boldsymbol{x}$. Upon triangulation of the surface $\Gamma$ of the domain $\Omega$ with triangular elements (which we naturally get from a triangulation of $\Omega$ with tetrahedral elements) and discretization of $\varphi _1$ and $\varphi _2$ we can rewrite Eq. (3.47) as
\begin{displaymath}
\mbox{\boldmath$\varphi $} _2 = B \mbox{\boldmath$\varphi $} _1
\end{displaymath} (3.48)

with the boundary matrix $B$, which is a dense matrix with a size of $n_b \times n_b$ elements, where $n_b$ is the number of nodes on the surface $\Gamma$.

The discretization of the scalar double layer operator in Eq. (3.47) has been derived by Lindholm [43]:

\begin{displaymath}
\int_\Gamma \varphi _1(\boldsymbol{y})
\frac{\partial G(\b...
..._{t \in \Gamma} \sum_{i=1}^{3} L_{t,i} \varphi _{1,t,i} \quad,
\end{displaymath} (3.49)

where $t$ runs over all triangles on the surface $\Gamma$ of the domain $\Omega$ and $i$ runs over the three nodes of each triangle.

In order to calculate the matrix entries of $B$ element by element (rather triangle by triangle) we use the local coordinates defined in Fig. 3.1.

Figure 3.1: Local coordinate system and various vectors required for the discretization of the boundary integral Eq. (3.47) [43].
\begin{figure}\center
\scalebox{0.7}{%
\input{fig/fem/lindholm.pstex_t}}\end{figure}

\begin{displaymath}
L_{t,i}=\frac{s_{i+1}}{8\pi\vert t\vert}
\left(
\eta_{i+1} S_t - \zeta \sum_{j=1}^3 \gamma_{ij} P_j
\right)
\end{displaymath} (3.50)


$\displaystyle \mbox{\boldmath$\rho$}_i$ $\textstyle =$ $\displaystyle \boldsymbol{r}_j -\boldsymbol{r}$ (3.51)
$\displaystyle s_{i}$ $\textstyle =$ $\displaystyle \vert\mbox{\boldmath$\rho$}_{j+1}-\mbox{\boldmath$\rho$}_j\vert$ (3.52)
$\displaystyle \eta_{i}$ $\textstyle =$ $\displaystyle \mbox{\boldmath$\hat{\eta}$}_j \cdot \mbox{\boldmath$\rho$}_j$ (3.53)
$\displaystyle \zeta$ $\textstyle =$ $\displaystyle \mbox{\boldmath$\hat\zeta$}\cdot\mbox{\boldmath$\rho$}_j$ (3.54)
$\displaystyle \gamma_{ij}$ $\textstyle =$ $\displaystyle \mbox{\boldmath$\hat\xi_{i+1}$} \cdot \mbox{\boldmath$\hat\xi$}_j$ (3.55)
$\displaystyle P_j$ $\textstyle =$ $\displaystyle \ln\frac{\rho_j+\rho_{j+1}+s_j}{\rho_j+\rho_{j+1}-s_j}$ (3.56)

$\vert t\vert$ denotes the area of triangle $t$ and $S_t$ the solid angle subtended by triangle $t$ at the ``observation point'' $\boldsymbol{r}$, which is given by
\begin{displaymath}
S_t =
2 \cdot \mathrm{sgn}(\zeta) \cdot \arccos
\left(
\...
...oldmath$\rho$}_1\cdot\mbox{\boldmath$\rho$}_2)
}
}
\right).
\end{displaymath} (3.57)

In order to calculate the demagnetizing field, we have to perform the following steps:


Initialization

  1. Discretize Eq. (3.41).
  2. Calculate the boundary matrix in Eq. (3.48).

Solution

  1. Solve Eq. (3.41) for a given magnetization distribution $\boldsymbol{J}$ using the standard FE method.
  2. Calculate $\varphi _2$ on the boundary $\Gamma$ using Eq. (3.48) to get the values for the Dirichlet boundary conditions.
  3. Calculate $\varphi _2$ in the whole domain $\Omega$ using Eq. (3.43) with Dirichlet boundary values.
  4. Calculate $\boldsymbol{H}_\mathrm{demag}=-\nabla(\varphi _1+\varphi _2)$.


next up previous contents
Next: 3.5 Effective Field Up: 3. Finite Element Micromagnetics Previous: 3.3.3 Zeeman Energy   Contents
Werner Scholz 2003-06-08