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

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

$$\Delta \varphi = 0 \quad \mathrm{for } \boldsymbol{x} \not\in \Omega$$ (3.35)

with the boundary conditions at the boundary $\Gamma$ of $\Omega$

$$\mathrm{Div }\varphi = 0$$ (3.36)

and

$$\mathrm{Div }\frac{\partial \varphi }{\partial \boldsymbol{n}} = -\boldsymbol{n} \cdot \boldsymbol{J} \quad.$$ (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

$$\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,$$ (3.38)

which can again be written in matrix-vector format as

$$\boldsymbol{d} = D \cdot \boldsymbol{u}$$ (3.39)

with

$$D_{j,3j+k}=\int_\Omega u_{j,k} \nabla_k \eta_j \quad,$$ (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

$$\Delta \varphi _1 = \nabla \cdot \boldsymbol{J}$$ (3.41)

with the boundary condition

$$\frac{\partial \varphi _1}{\partial \boldsymbol{n}}= \boldsymbol{n} \cdot \boldsymbol{J} \quad.$$ (3.42)

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

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

$$\Delta \varphi _2=0$$ (3.43)

with

$$\mathrm{Div }\varphi _2 = \varphi _1$$ (3.44)

and

$$\mathrm{Div }\frac{\partial \varphi _2}{\partial \boldsymbol{n}} = 0 \quad.$$ (3.45)

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

Potential theory tells us that

$$\varphi _2(\boldsymbol{x})=\frac{1}{4\pi}\int_\Gamma \varph... ...{y})}{\partial \boldsymbol{n}(\boldsymbol{y})} d{a} \quad,$$ (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

$$\varphi _2(\boldsymbol{x})=\frac{1}{4\pi} \int_\Gamma \var... ...mbol{x})}{4 \pi}-1 \right) \varphi _1(\boldsymbol{x}) \quad,$$ (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

$$\mbox{\boldmath$\varphi $} _2 = B \mbox{\boldmath$\varphi $} _1$$ (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]:

$$\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,$$ (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].

$$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)$$ (3.50)

$$\mbox{\boldmath$\rho$}_i = \boldsymbol{r}_j -\boldsymbol{r}$$ (3.51)

$$s_{i} = \vert\mbox{\boldmath$\rho$}_{j+1}-\mbox{\boldmath$\rho$}_j\vert$$ (3.52)

$$\eta_{i} = \mbox{\boldmath$\hat{\eta}$}_j \cdot \mbox{\boldmath$\rho$}_j$$ (3.53)

$$\zeta = \mbox{\boldmath$\hat\zeta$}\cdot\mbox{\boldmath$\rho$}_j$$ (3.54)

$$\gamma_{ij} = \mbox{\boldmath$\hat\xi_{i+1}$} \cdot \mbox{\boldmath$\hat\xi$}_j$$ (3.55)

$$P_j = \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

$$S_t = 2 \cdot \mathrm{sgn}(\zeta) \cdot \arccos \left( \... ...oldmath$\rho$}_1\cdot\mbox{\boldmath$\rho$}_2) } } \right).$$ (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)$.