Introduction

In this numerical experiment we will calculate the magnetic field, $\vec{B}$, induced by a thick spherical coil with a spherical magnetic core. The figure below illustrates a cross section of the coil.

The calculation is done in four stages:

Calculate the current vector potential, $\vec{T}$, for a given volume free-current density $\vec{J}_f$ (StaticVectorSolver::Solver1).
Check the calculated $\vec{T}$ by converting it back into $\vec{J}_f$ (StaticVectorSolver::ProjectTtoJ).
Use the calculated current vector potential to solve for $\vec{A}$ (StaticVectorSolver::Solver2).
Convert computed $\vec{A}$ into the magnetic field $\vec{B}$ (StaticVectorSolver::ProjectAtoB).

The Bossavit's diagrams in the figure below illustrate these four stages.

This numerical experiment is carried out in the three-dimensional space. As implied by the Bossavit's diagrams above, see also this discussion, the vector potentials, $\vec{T}$ and $\vec{A}$, belong to the $H(\text{curl})$ function space, and, thus, must be modeled by the FE_Nedelec finite elements. Similarly, $\vec{J}_f$ and $\vec{B}$ belong to the $H(\text{div})$ function space and, thus, must be modeled by the FE_RaviartThomas finite elements.

As per usual, we will execute the experiment on a family of progressively refined meshes and compare the numerical results to closed-form analytical expression of the exact solution. We will vary the degree of the finite elements and observe that the higher degree of the finite elements imply faster convergence of the $L^2$ error norms. As soon as we do not gauge the vector potentials ( $\vec{T}$ and $\vec{A}$) explicitly, their conservative parts will remain unknown. For this reason, it does not make much sense to visualize them or evaluate them by computing error norms. We will evaluate the computed vector potentials indirectly by computing $L^2$ error norms of the fields derived from them ( $\vec{J}_f$ and $\vec{B}$).

Implementation

The problem and the problem domain

The following textbook problem is considered (Section 6.2.6. Spherical coil with magnetic core).

Volume free-current density, $\vec{J}_f$, is confined to an interior of a spherical shell. This shell is shown in blue color in the figure above. The inner and outer radii of the shell equal $ a_2 $ and $ b_2 $, respectively. The free-current density can be expressed in a spherical coordinate system as
\[ \vec{J}_f = \left\{ \begin{aligned} & 0 &&\text{if }&& r < a_2\\ & K_0 r \sin(\theta) \hat{\phi} &&\text{if }&& a_2 \le r \le b_2\\ & 0 &&\text{if }&& r > b_2.\\ \end{aligned} \right. \]

There is a magnetic core of permeability $\mu_1$ inside the coil. The core is shaped as a spherical shell. This shell is shown in green color in the figure above.The inner and outer radii of the shell equal $a_1$ and $b_1$, respectively. The permeability of the entire space can be expressed as
\begin{equation} \mu = \left\{ \begin{aligned} & \mu_0 &\text{if } \text{ } & r < a_1\\ & \mu_1 &\text{if } \text{ } & a_1 \le r \le b_1\\ & \mu_0 &\text{if } \text{ } & r > b_1. \end{aligned} \right. \end{equation}

The shell filled with free current is concentric with the shell filled with magnetic material.

We would like to calculate the magnetic field, $\vec{B}$, induced by the coil observing the compatibility condition for the curl-curl equation.

The figure below illustrates the geometry of the coil and the corresponding three-dimensional problem domain.

The Dirichlet boundary condition,

\[ \hat{n} \times \vec{T} = 0, \]

is applied when computing the current vector potential. The boundary ID in this case equals $1$. This boundary ID is set in the geo files. The Robin boundary condition (first-order ABC),

\[ \frac{1}{\mu} \hat{n} \times \bigg(\vec{\nabla} \times \vec{A}\bigg) + \frac{1}{\mu r} \hat{n} \times \bigg(\hat{n} \times \vec{A}\bigg) = 0, \]

is applied when computing the magnetic vector potential. The boundary ID in this case equals $2$. This boundary ID is set by the code of the program. That is, during the simulation the boundary ID has two values: $1$ and $2$.

The exact solution

The problem has the following closed-form analytical solution.

The total magnetic field induced by the coil with the core can be expressed as a sum of the magnetic field induced by the free-current density in free space, $\vec{B}_J$, and the magnetic field induced by the magnetization of the magnetic core, $\vec{B}_{\mu}$,

\[ \vec{B} = \vec{B}_J +\vec{B}_{\mu}. \]

The term $ \vec{B}_J $ can be expressed in spherical coordinates as

\[ \vec{B}_J = \left\{ \begin{aligned} & \frac{1}{2} (b_2^2-a_2^2) \vec{F}_1(\theta) &&\text{if }&& r \le a_2\\ & \frac{1}{2} (b_2^2-r^2) \vec{F}_1(\theta) + \frac{1}{5} (r^5-a_2^5) \vec{F}_2(r, \theta) && \text{if } && a_2 \le r \le b_2\\ & \frac{1}{5} (b_2^5-a_2^5) \vec{F}_2(r, \theta) &&\text{if }&& r \ge b_2,\\ \end{aligned} \right. \]

where

\[ \begin{aligned} &\vec{F}_1(\theta) = \dfrac{2}{3}\mu_0 K_0\big(\cos(\theta)\hat{r}-\sin(\theta)\hat{\theta}\big), \\ &\vec{F}_2(r, \theta) = \dfrac{2}{3}\mu_0 K_0 \frac{1}{r^3}\big(\cos(\theta)\hat{r}+\frac{1}{2}\sin(\theta)\hat{\theta}\big). \end{aligned} \]

We can convert these expressions into Cartesian coordinate system by noting that (See Appendix A)

\[ r = \sqrt{x^2+y^2+z^2}, \]

\[ \cos(\theta) = \frac{z}{r}, \]

\[ \sin(\theta) = \frac{\sqrt{x^2+y^2}}{r}, \]

\[ \cos(\phi) = \frac{x}{\sqrt{x^2+y^2}}, \]

\[ \sin(\phi) = \frac{y}{\sqrt{x^2+y^2}}, \]

\[ \hat{r} = \frac{1}{r} (x\hat{i}+y\hat{j}+z\hat{k}), \]

and

\[ \hat{\theta}= \cos(\theta)\cos{\phi}\hat{i} + \cos(\theta)\sin{\phi}\hat{j} - \sin(\theta) \hat{k}. \]

The magnetic field induced by the magnetic core can be expressed in Cartesian coordinates as

\[ \vec{B}_{\mu} = - \mu\vec{\nabla} \Psi - \mu_0 H_0 \hat{k}, \]

where

\[ H_0 = \frac{1}{3} K_0 (b_2^2-a_2^2), \]

\begin{equation} \mu = \left\{ \begin{aligned} & \mu_0 &\text{if } \text{ } & r < a_1\\ & \mu_1 &\text{if } \text{ } & a_1 \le r \le b_1\\ & \mu_0 &\text{if } \text{ } & r > b_1, \end{aligned} \right. \end{equation}

and

\begin{equation} \vec{\nabla} \Psi = \left\{ \begin{aligned} &\delta_1 \hat{k} &\text{if } \text{ } & r \le a_1\\ &-3\gamma_1\frac{xz}{r^5}\hat{i} -3\gamma_1\frac{yz}{r^5}\hat{j} + \bigg(\beta_1 + \gamma_1\frac{1}{r^3} - 3\gamma_1\frac{z^2}{r^5} \bigg)\hat{k} &\text{if } \text{ } & a_1 \le r \le b_1\\ &-3\alpha_1\frac{xz}{r^5}\hat{i} -3\alpha_1\frac{yz}{r^5}\hat{j} + \bigg(-H_0 + \alpha_1\frac{1}{r^3} - 3\alpha_1\frac{z^2}{r^5} \bigg)\hat{k} &\text{if } \text{ } & r \ge b_1. \end{aligned} \right. \end{equation}

The following notations were used in the last expression:

\begin{equation} \begin{aligned} &\Omega = \frac{(\mu_r - 1)}{(\mu_r + 2)} \frac{a_1^3}{b_1^3}, \\ &\gamma_1 = \frac{-3b_1^3 H_0 \Omega}{(2\mu_r+1)-2 (\mu_r-1) \Omega}, \\ &\beta_1 = \frac{(2\mu_r+1) \gamma_1}{(\mu-1)a_1^3}, \\ &\alpha_1 = \frac{-b_1^3 H_0+2\mu_r\gamma_1-\mu_r b_1^3\beta_1}{2}, \\ &\delta_1 = \frac{\mu_r a_1^3\beta_1-2\mu_r\gamma_1}{a_1^3}. \end{aligned} \end{equation}

The relative permeability in the last five equations is computed as

\[ \mu_r = \frac{\mu_1}{\mu_0}. \]

The mesh

The figure below illustrates the mesh.

The main elements of the mesh are listed below.

The cube with a side of $2d_1$ in the middle of the mesh. The purpose of this cube is instrumental. This is the only way to construct a mesh that has spherical interfaces and contains no tetrahedral cells as the cells around the origin must be hexahedral.
The spherical shell nr.1. This shell represents the magnetic core. The permeability of this region is $\mu_1$. The permeability outside this region is $\mu_0$.
The spherical shell nr.2. This region contains the free current, $ \vec{J}_f = K_0 r \sin(\theta) \hat{\phi} $. There is no free current outside this region.
The sphere nr.1. The $L^2$ error norms are computed inside the region delineated by this surface.
The sphere nr.2. The boundary of the problem domain. Represents infinity. The boundary conditions (Dirichlet, Robin) is applied to this surface.

The boundary value problem

At the 0th stage of the experiment the current vector potential, $\vec{T}$, is calculated by solving the following boundary value problem:

\begin{equation} \begin{array}{lrcll} \text{ } & \vec{\nabla}\times\bigg(\vec{\nabla}\times\vec{T}\bigg) +\eta^2\vec{T} = \vec{\nabla}\times\vec{J}_f & \text{in} & \Omega & \text{(i)},\\ \text{(e)}& \hat{n} \times \vec{T} = 0 & \text{on} & \Gamma_{R1} & \text{(ii)}. \\ \end{array} \end{equation}
At the 1st stage of the experiment the current vector potential calculated at the 0th stage, $\vec{T}$, is converted back into free-current density as
\[ \vec{J}_{f} = \vec{\nabla} \times \vec{T} \]
for the purpose of evaluation.
At the 2nd stage of the experiment the following boundary value problem is solved:
\begin{equation} \begin{array}{lrcll} \text{ } & \vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A}\bigg) + \eta^2 \vec{A} = \vec{\nabla}\times\vec{T} & \text{in} & \Omega & \text{(i)},\\ \text{(n)}& \dfrac{1}{\mu}\hat{n}\times\bigg(\vec{\nabla}\times\vec{A}\bigg) + \gamma \hat{n}\times\bigg(\hat{n}\times\vec{A}\bigg) =0 &\text{on} & \Gamma_{R1} & \text{(ii)}, \\ \text{(e)}&\hat{n}\times\vec{A}_{+} = \hat{n}\times\vec{A}_{-}&\text{on}&\Gamma_{I1}\cup\Gamma_{I2}&\text{(iii)},\\ \text{(n)}&\dfrac{1}{\mu}_{+}\hat{n}\times \bigg( \vec{\nabla} \times \vec{A}_{+} \bigg) - \dfrac{1}{\mu}_{-}\hat{n}\times \bigg( \vec{\nabla} \times \vec{A}_{-} \bigg) = 0 & \text{on} & \Gamma_{I1}\cup\Gamma_{I2} & \text{(iv)}, \end{array} \end{equation}
where
\[ \gamma = \frac{1}{\mu r}. \]
The current vector potential, $\vec{T}$, on the right-hand side of (i) is computed numerically at the 0th stage.
At the 3rd stage of the experiment the magnetic vector potential calculated at the 2nd stage, $\vec{A}$, is converted into magnetic field as
\[ \vec{B} = \vec{\nabla} \times \vec{A}. \]

The program

This numerical experiment is implemented in accordance with the base structure described in here. All computations are done in four stages. The stages are labeled from 0 to 3. The following algorithms where used:

At the 0th stage the current vector potential, $\vec{T}$, is computed with a help of the StaticVectorSolver::Solver1 class template.
At the 1st stage the current vector potential computed at the 0th stage is converted back into volume free-current density, $\vec{J}_{f}$, with a help of the StaticVectorSolver::ProjectTtoJ class template.
At the 2nd stage the magnetic vector potential, $\vec{A}$, is computed with a help of the StaticVectorSolver::Solver2 class template. The current vector potential, $\vec{T}$, computed at the 0th stage is fed to the algorithm.
At the 3rd stage the magnetic vector potential computed at the 2nd stage is converted to the magnetic field, $\vec{B}$, with a help of the StaticVectorSolver::ProjectAtoB class template.

The Dirichlet and Robin (ABC) boundary conditions are applied at the 0th and the 2nd stages, respectively. The boundary ID must satisfy the boundary ID convension. For this reason, two boundary ID's are used. At the 0th stage the boundary ID equals $1$. This boundary ID is set in the geo file as the following:

Physical Surface(1) = {115, 479, 665, 297, 847, 1029};

It is loaded together with the mesh. The boundary ID equals $2$ at the second stage. The boundary ID is changed from $1$ to $2$ in the constructor of the SolverSSOLIII_A class:

...
 
  SolverSSOLIII_A(unsigned int p,
                unsigned int mapping_degree,
                unsigned int r,
                const Triangulation<3>& triangulation_T,
                const DoFHandler<3>& dof_handler_T,
                const Vector<double>& solution_T,
                std::string fname)
  : Solver2<3, 2>(p,
                  mapping_degree,
                  triangulation_T,
                  dof_handler_T,
                  solution_T,
                  2,
                  0.0,
                  fname,
                  nullptr,
                  SettingsSSOLIII::print_time_tables,
                  false,
                  true)
  , r(r)
  , fname(fname)
{
  for (auto cell : Solver2<3, 2>::triangulation_T.active_cell_iterators()) {
    for (unsigned int f = 0; f < GeometryInfo<3>::faces_per_cell; f++) {
      if (cell->face(f)->at_boundary() && (cell->face(f)->boundary_id() == 1))
        cell->face(f)->set_boundary_id(2);
    }
  }
 
  Solver2<3, 2>::run();
}
 
...

Note that the free-current density, $\vec{J}_f$, equals zero on the boundary by definition of the problem. The tangential component of the current vector potential, $\vec{T}$, equals zero on the boundary as we apply the homogeneous Dirichlet boundary condition when computing it, see the first boundary value problem above. For these two reasons we can neglect integrals $I_{b3-2}$ in the Recipe for static vector solver in 3D (current vect. potential) and in Recipe for static vector solver in 3D. To do so, we pass type_of_pde_rhs = 2 to the constructors of both, StaticVectorSolver::Solver1 and StaticVectorSolver::Solver2.

The build process generates one executable file: ssol-iii-sphere. To rebuild it change into ssol-iii/build/Release directory and execute the following:

user@computer .../ssol-iii/build/Release$ ./clean

user@computer .../ssol-iii/build/Release$ ./build

Then the executable file must be executed again. This can be done as the following:

user@computer .../ssol-iii/build/Release$ cd ../../bin/Release

user@computer .../ssol-iii/bin/Release$ ./ssol-iii-sphere

This will generate vtu files in the ssol-iii/bin/Release/Data directory. They can be viewed with a help of a fresh version of ParaView. The VisIt visualization software will not do as this numerical experiment uses second-degree mapping.

Note that executable files require mesh files to be present in the ssol-iii/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into ssol-iii/gmsh directory and executing the following:

user@computer .../ssol-iii/gmsh$ ./clean

user@computer .../ssol-iii/gmsh$ ./build

This will generate a set of globally refined meshes in ssol-iii/gmsh/data directory.

The SettingsSSOLIII class allows switching on three useful features: printing time tables on the computer screen, logging convergence data of the conjugate gradient solver, and saving the exact solution into the vtu files next to the numerical solution. The last can be useful when debugging.

Simulation results

The simulation was done with the following settings: $d_1 = 0.1[m]$, $a_1 = 0.3[m]$, $b_1=0.6[m]$, $a_2 = 0.9[m]$, $b_2 = 1.2[m]$, $d_2 = 2.0[m]$, $d_3 = 3.0[m]$, $K_0=1[A\cdot m^{-2}]$, and $\mu_1 = 4\mu_0$.

The first figure on this page illustrates a vector plot of the computed free-current density, $\vec{J}_f$, and the magnetic field, $\vec{B}$.

The free-current density

The figure below illustrates two slices of the magnitude of free-current density, $|\vec{J}_{f}|$, calculated with finite elements of degree $p = 2$ on the mesh sphere_r9.msh.

The table below illustrates the convergence of the $L^2$ norm of the error made in computing the free-current density, $\vec{J}_{f}$.

p	r	cells	dofs	$\\|e\\|_{L^2}$	$\alpha_{L^2}$
0	6	4625	13950	1.66e-01	-
0	7	7992	24084	1.38e-01	0.99
0	8	12691	38220	1.19e-01	0.99
0	9	18944	57024	1.04e-01	0.99
1	6	4625	111300	8.12e-04	-
1	7	7992	192240	4.97e-04	2.69
1	8	12691	305172	3.32e-04	2.61
1	9	18944	455424	2.37e-04	2.54
2	6	4625	375300	6.78e-04	-
2	7	7992	648324	3.94e-04	2.97
2	8	12691	1029294	2.49e-04	2.98
2	9	18944	1536192	1.67e-04	2.99

The following notations were used in the table header:

p - the degree of the FE_Nedelec finite elements that were used to model the current vector potential, $\vec{T}$, and the degree of the FE_RaviartThomas finite elements that were used to model the free-current density, $\vec{J}_f$.
r - the number of nodes on transfinite lines, see discussion in here.
cells - the total amount of active cells.
dofs - the amount of degrees of freedom in the active cells.
$\|e\|_{L^2}$ - the $L^2$ error norm.
$\alpha_{L^2}$ - the order of convergence of the $L^2$ error norm.

The magnetic field

The figure below illustrates two slices of the magnitude of magnetic field, $|\vec{B}|$, calculated with finite elements of degree $p = 2$ on the mesh sphere_r9.msh.

The table below illustrates the convergence of the $L^2$ norm of the error made in computing the magnetic field, $\vec{B}$.

p	r	cells	dofs	$\\|e\\|_{L^2}$	$\alpha_{L^2}$
0	6	4625	13950	8.84e-08	-
0	7	7992	24084	7.36e-08	1.00
0	8	12691	38220	6.30e-08	1.01
0	9	18944	57024	5.51e-08	1.00
1	6	4625	111300	4.41e-09	-
1	7	7992	192240	3.11e-09	1.91
1	8	12691	305172	2.23e-09	2.18
1	9	18944	455424	1.71e-09	1.96
2	6	4625	375300	1.84e-10	-
2	7	7992	648324	1.03e-10	3.21
2	8	12691	1029294	6.08e-11	3.40
2	9	18944	1536192	4.04e-11	3.07