 |
Logbook
(07-04-2025)
Static problems
|
|
|
Logbook
(07-04-2025)
Static problems
|
|
Classes:
1) BatchSSOLI
2) SolverSSOLI
3) ExactSolutionSSOLI_B
4) SettingsSSOLI
Files:
1) ssol-i/src/main.cpp
2) ssol-i/include/solver.hpp
3) ssol-i/src/solver.cpp
4) ssol-i/include/exact_solution.hpp
5) ssol-i/src/exact_solution.cpp
6) ssol-i/src/static_vector_input.cpp
7) ssol-i/include/settings.hpp
In this numerical experiment we will calculate the magnetic field, \(\vec{B}\), induced by a thin spherical coil. The coil consists of an infinitesimally thin layer of free current flowing on a spherical surface. The current has only azimuthal component. The figure below illustrates the coil.
The calculation is done in two steps:
The Bossavit's diagrams in the figure below illustrate these two steps.
This numerical experiment is carried out in the three-dimensional space. As implied by the Bossavit's diagrams above, see also this discussion, the magnetic vector potential, \(\vec{A}\), belongs to the \(H(\text{curl})\) function space and, thus, must be modeled by the FE_Nedelec finite elements. Similarly, the magnetic field, \(\vec{B}\), belongs 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 result 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 norm. As soon as we do not gauge the vector potential explicitly, its conservative part will remain unknown. For this reason, it does not make much sense to visualize it or evaluate it by computing error norms. We will evaluate the computed vector potential indirectly by computing \(L^2\) error norm of the magnetic field derived from it.
The following textbook problem is considered (Section 6.2.4. Thin spherical coil).
A free current flows on a surface of a sphere of a radius \(a\). The sphere is shown in blue color in the figure above. The surface free-current density on the sphere has the following form in the spherical coordinate system
\[ \vec{K}_f = K_0 a \sin(\theta) \hat{\phi}, \]
where \( \theta \) is the polar angle, \(\phi\) is the azimuthal angle, and \( K_0 \) is a scalar.
The permeability equals \(\mu_0\) at all points in space.
We would like to calculate the magnetic field, \(\vec{B}\) inside and outside the sphere.
The figure below illustrates the geometry of the coil and the corresponding three-dimensional problem domain, \(\Omega\).
The Robin boundary condition (the first-order ABC),
\[ \frac{1}{\mu_0} \hat{n} \times \bigg(\vec{\nabla} \times \vec{A}\bigg) + \frac{1}{\mu_0 r} \hat{n} \times \bigg(\hat{n} \times \vec{A}\bigg) = 0, \]
is applied to the boundary of the problem domain.
The problem has the following closed-form analytical solution.
\[ \vec{B} = \left\{ \begin{aligned} &\dfrac{2}{3}\mu_0 K_0 a\big(\cos(\theta)\hat{r}-\sin(\theta)\hat{\theta}\big)&&\text{if }&& r < a\\ &\dfrac{2}{3}\mu_0 K_0 \frac{a^4}{r^3}\big(\cos(\theta)\hat{r}+\frac{1}{2}\sin(\theta)\hat{\theta} \big)&&\text{if }&& r>a. \\ \end{aligned} \right. \]
The figure below illustrates the mesh.
The main elements of the mesh are listed below.
The following boundary value problem is solved in this numerical experiment with a help of the StaticVectorSolver::Solver1 class template:
\begin{equation} \begin{array}{lrcll} \text{ } & \vec{\nabla}\times\bigg(\dfrac{1}{\mu_0} \vec{\nabla}\times\vec{A}\bigg) + \eta^2 \vec{A} = 0 & \text{in} & \Omega & \text{(i)}, \\ \text{(n)}& \dfrac{1}{\mu_0} \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}&\text{(iii)}, \\ \text{(n)}&\dfrac{1}{\mu_{0}}\hat{n}\times \bigg( \vec{\nabla} \times \vec{A}_{+} \bigg) - \dfrac{1}{\mu_{0}}\hat{n}\times \bigg( \vec{\nabla} \times \vec{A}_{-} \bigg) = \vec{K}_f & \text{on} & \Gamma_{I1} & \text{(iv)}. \end{array} \end{equation}
where
\[ \gamma = \dfrac{1}{\mu_0 r}. \]
Then the numerically calculated magnetic vector potential, \(\vec{A}\), is converted into a magnetic field as
\[ \vec{B} = \vec{\nabla} \times \vec{A}. \]
with a help of the algorithm of the StaticVectorSolver::ProjectAtoB class template.
This numerical experiment is implemented in accordance with the base structure described in here.
The simulations are done by algorithms described by two class templates, StaticVectorSolver::Solver1 and StaticVectorSolver::ProjectAtoB. These two class templates do not share the input parameters declared in static_vector_input.hpp. For this reason, there is no need to invoke the stage mechanism. In this numerical experiment we keep the stage parameter at its default value, i.e., we simply ignore it.
The build process generates one executable file: ssol-i-sphere. To rebuild it change into ssol-i/build/Release directory and execute the following:
Then the executable file must be executed again. This can be done as the following:
This will generate vtu files in the ssol-i/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-i/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into ssol-i/gmsh directory and executing the following:
This will generate a set of globally refined meshes in ssol-i/gmsh/data directory.
The SettingsSSOLI 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.
To work properly, the StaticVectorSolver::Solver1 needs some way of discriminating the cell faces that belong to the interface \(\Gamma_{I1}\), so it can distribute the surface free-current density, \(\vec{K}_f\), on it. This discrimination is accomplished by invoking the mechanism of cell and face user IDs which is available in deal.II. The StaticVectorSolver::Solver1 class template calculates integrals related to \(\vec{K}_f\) only if the user cell ID and the user face ID, do not equal zero. Furthermore, StaticVectorSolver::Solver1 passes the values of the user cell and face IDs to StaticVectorSolver::FreeSurfaceCurrent::value_list, so the algorithm implemented by this function knows which cell face is currently being processed. Therefore, from the user's standpoint \(\vec{K}_f\), is handled in two places in the code: (i) in the private function SolverSSOLI::mark_materials where cells and faces on the interface are marked; and (ii) StaticVectorSolver::FreeSurfaceCurrent, where proper values of \(\vec{K}_f\) are calculated.
Note that each face on an interface always belongs to two cells. Therefore, a care must be taken that each face on the interface is handled only ones. This is done by marking a face on the interface only if it belongs to a cell center of which is located within the sphere of radius \(a\). The code snippet below explains this in more detail.
Then the StaticVectorSolver::FreeSurfaceCurrent<3>::value_list calculates the values of \(\vec{K}_f\):
The simulation was done with the following settings: \(d_1 = 0.1[m]\), \(a = 1.2[m]\), \(d_2 = 2.0[m]\), \(d_3 = 3.0[m]\), and \(K_0=1.0[A \cdot m^{-1}]\).
The first figure on this page illustrates a vector plot of the computed the magnetic field, \(\vec{B}\).
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 | 2375 | 7200 | 3.00e-07 | - |
| 0 | 7 | 4104 | 12420 | 2.48e-07 | 1.05 |
| 0 | 8 | 6517 | 19698 | 2.12e-07 | 1.02 |
| 0 | 9 | 9728 | 29376 | 1.84e-07 | 1.07 |
| 1 | 6 | 2375 | 57300 | 2.45e-08 | - |
| 1 | 7 | 4104 | 98928 | 1.94e-08 | 1.28 |
| 1 | 8 | 6517 | 156996 | 1.52e-08 | 1.58 |
| 1 | 9 | 9728 | 234240 | 1.30e-08 | 1.17 |
| 2 | 6 | 2375 | 193050 | 5.16e-10 | - |
| 2 | 7 | 4104 | 333396 | 2.74e-10 | 3.48 |
| 2 | 8 | 6517 | 529200 | 1.40e-10 | 4.34 |
| 2 | 9 | 9728 | 789696 | 8.73e-11 | 3.55 |
The following notations were used in the table header:
1.9.1