 |
Logbook
(07-04-2025)
Static problems
|
|
|
Logbook
(07-04-2025)
Static problems
|
|
Classes:
1) BatchSLDII
2) SolverSLDII
3) ExactSolutionSLDII_THETA
4) ExactSolutionSLDII_H
5) SettingsSLDII
Files:
1) sld-ii/src/main.cpp
2) sld-ii/include/solver.hpp
3) sld-ii/src/solver.cpp
4) sld-ii/include/exact_solution.hpp
5) sld-ii/src/exact_solution.cpp
6) sld-ii/src/static_scalar_input.cpp
7) sld-ii/include/settings.hpp
In some problems a portion of the auxiliary field \(\vec{H}\) is solenoidal. This happens when the free currents, \(\vec{J}_f\), cannot be excluded from the problem domain. In such cases, the field lines of \(\vec{H}\) tend to curl around the free currents \(\vec{J}_f\). The last makes application of the total magnetic scalar potential utterly impossible as one cannot describe a solenoidal field by a gradient of a scalar potential. There is, however, a way around (see Section 3.2. Reduced magnetic scalar potential): we split the vector filed \(\vec{H}\) into two components,
\[ \vec{H} = \vec{H}_S + \vec{H}_C, \]
where \(\vec{H}_S\) and \(\vec{H}_C\) is the solenoidal and conservative parts of the vector field \(\vec{H}\). The solenoidal part, \(\vec{H}_S\), is the field induced by the free currents, \(\vec{J}_f\), in free space. There is a number of ways of calculating the solenoidal part of \(\vec{H}\). The Biot-Sawart law is one of them. In any case, we assume that \(\vec{H}_S\) is known in advance and is not a part of the problem. We derive the conservative part, \(\vec{H}_C\), as a negative gradient of a reduced magnetic scalar potential, \(\Theta\):
\[ \vec{H}_C = - \vec{\nabla} \Theta. \]
The reduced magnetic scalar potential, in turn, is described by the static scalar boundary value problem which can be solved with a help of the StaticScalarSolver::Solver.
This numerical experiment is made as an illustration of application of the reduced magnetic scalar potential, \(\Theta\). The illustration is done by solving the two problems that were considered in the Magnetostatic shield -1 (sld-i) numerical experiment. The free current density, \(\vec{J}_f\), in these two problems is outside the problem domain. So, strictly speaking, we do not have a good reason for resorting to the reduced magnetic scalar potential. However, as soon as we will concentrate on calculating the H-field induced by magnetization, it does not really matter if the applied H-field is solenoidal or conservative. In this numerical experiment we will split the total \(\vec{H}\) field as the following
\[ \vec{H} = \vec{H}_0 + \vec{H}_C, \]
where \(\vec{H}_0\) is the known applied field (it is not solenoidal), and concentrate on calculating \(\vec{H}_C\), which is the auxiliary H-field induced by magnetization.
In the two-dimensional version of the experiment we will reuse the cylindrical magnetic shield problem as well as the corresponding problem domain and the mesh. We express the exact closed-form analytical solution and its gradient as
\[ \Theta = \Psi + H_0 x \]
and
\[ \vec{\nabla} \Theta = \vec{\nabla} \Psi + H_0, \]
where \(\Psi\) and \(\vec{\nabla} \Psi \) are the exact solutions presented here.
In the three-dimensional version of the experiment we will reuse the spherical magnetic shield problem as well as the corresponding problem domain and the mesh. We express the exact closed-form analytical solution and its gradient as
\[ \Theta = \Psi + H_0 z \]
and
\[ \vec{\nabla} \Theta = \vec{\nabla} \Psi + H_0, \]
where \(\Psi\) and \(\vec{\nabla} \Psi \) are the exact solutions presented here.
In the approach based on the total magnetic scalar potential, (it is described in the Magnetostatic shield - 1 (sld-i) numerical experiment) the information on the sources of the magnetic field is introduced into the problem domain by means of the Dirichlet boundary condition. This time we introduce the information on the sources of the magnetic field by means of the reduced magnetic surface charge density, \(\kappa_r\). The reduced magnetic surface charge density is computed as
\[ \kappa_r = - (\mu_{+}-\mu_{-})\big(\hat{n}\cdot\vec{H}_0\big). \]
As soon as the permeability is constant inside and outside the shield, the reduced magnetic volume charge density, \( \rho_r \), equals zero everywhere. The H-field induced by the magnetization vanish at infinity. Therefore, the reduced magnetic scalar potential at infinity is constant. We are free to choose any constant, so we choose \(\Theta=0\) at infinity. Then the boundary value problem solved in the two- and three- dimensional versions of this experiment reads
\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big( \mu \vec{\nabla} \Theta \big)=0 & \text{in} & \Omega & \text{(i)}\\ \text{(e)} &\Theta = \eta & \text{on} & \Gamma_{D1} & \text{(ii)}\\ \text{(e)} &\Theta_{+} = \Theta_{-} & \text{on} & \Gamma_{I1} & \text{(iii)}\\ \text{(n)}&\mu_{+}\hat{n}\cdot\vec{\nabla}\Theta_{+}-\mu_{-}\hat{n}\cdot\vec{\nabla}\Theta_{-}= - \kappa_r & \text{on} & \Gamma_{I1}&\text{(iv)}\\ \text{(e)} &\Theta_{+} = \Theta_{-} & \text{on} & \Gamma_{I2} & \text{(v)}\\ \text{(n)}&\mu_{+}\hat{n}\cdot\vec{\nabla}\Theta_{+}-\mu_{-}\hat{n}\cdot\vec{\nabla}\Theta_{-}= - \kappa_r & \text{on} & \Gamma_{I2}&\text{(vi)}. \end{array} \end{equation}
The program can operate in two modes: "Dirichlet" and "Exact". The reduced magnetic scalar potential vanishes at infinity. For this reason it is appropriate to apply the homogeneous Dirichlet boundary condition,
\[ \Theta = 0, \]
to entire boundary in both problems. The program in the "Dirichlet" mode does exactly this. Both two- and three- dimensional problem domains in this numerical experiment are truncated. That is, the boundary of the problem domain is at a finite distance from the origin. The domains of the initial problems are unbounded. For this reason, application of the Dirichlet boundary condition will lead to truncation errors. We would like to observe the convergence rate without the truncation errors. In general, the truncation errors can be driven to be insignificantly small by increasing the size of the problem domain, i.e., moving the boundary farther away from the origin. This, however, will increase the simulation time. Instead, we will keep the sizes of the problem domains relatively small and use the exact boundary conditions. That is, we will apply the Dirichlet boundary condition and use the known exact solutions to compute the right-hand sides \(\eta\). This will allow us to switch off the truncation errors and observe the convergence rates free of the truncation error. This is what program does if switched into the "Exact" mode. The user can switch the program between the two modes by setting Setting::type_of_bc.
First, a class template (SolverSLDII) derived from StaticScalarSolver::Solver is used to obtain the numerical solution, \(\Theta\). After that, the numerical solution is fed to an object derived from the class template StaticScalarSolver::ProjectTHETAtoH (this is just another name for StaticScalarSolver::ProjectPSItoH) to calculate the auxiliary vector field \(\vec{H}_C\) as
\[ \vec{H}_C = - \vec{\nabla} \Theta. \]
This time we do not compute the magnetic field as the reduced magnetic field does not make any sense. We cannot compute a reduced magnetic field simply as
\[ \vec{B}_C = \mu \vec{H}_C = -\mu\vec{\nabla}\Theta. \]
Indeed, equations (iv) and (vi) in the boundary value problem above imply that the normal component of \(-\mu\vec{\nabla}\Theta\) is discontinuous on interfaces. A normal component of a magnetic field must be continuous on interfaces. Strictly speaking, \(-\mu\vec{\nabla}\Theta\) is not a magnetic field: it exhibits a wrong behavior on interfaces. So, the "reduced" magnetic field does not make much sense. One must compute the total magnetic field by projecting \(\mu(\vec{H}_S - \vec{\nabla}\Theta)\) onto \(H(\text{div})\) function space. The result should be identical to the magnetic fields computed in Magnetostatic shield - 1 (sld-i) numerical experiment.
The computer program of this numerical experiment differs from the computer program of the Magnetostatic shield - 1 numerical experiment at three points:
To work properly, the StaticScalarSolver::Solver needs a method for discriminating the faces that belong to the interfaces \(\Gamma_{I1}\) and \(\Gamma_{I2}\) so it can distribute the reduced magnetic surface charge density, \(\kappa_r\), on them. This discrimination is accomplished by invoking the mechanism of cell and face user ID which is available in deal.II. The StaticScalarSolver::Solver class template calculates the integral related to \(\kappa_r\) only if both user IDs, the user cell ID and the user face ID, do not equal zero. Furthermore, StaticScalarSolver::Solver passes the values of the user cell and face IDs to StaticScalarSolver::FreeSurfaceCharge::value_list so the algorithm implemented by this function knows which interface is currently being processed. Therefore, from the user's standpoint \(\kappa_r\) is handled in two places in the code: (i) the function SolverSLDII::mark_materials() marks the relevant cells and faces with user IDs when mesh is loaded, (ii) StaticScalarSolver::FreeSurfaceCharge::value_list, where proper values of \(\kappa_r\) are calculated based on cell and face user IDs. That is, we first mark the cells and faces,
and then choose the proper combinations of the cell user IDs and face user IDs,
Both, two- and three- dimensional versions of the experiment were executed under the following conditions: \(a = 0.2[m]\), \(b=0.4[m]\), \(d_2=0.8[m]\), \(d_3=2.0[m]\), \(\mu = 4.0 \mu_0\). All convergence tables were generated in the "Exact" mode, see above.
The figure below illustrates a contour plot of the reduced magnetic scalar potential, \(\Theta\), the output of an object of the SolverSLDII type.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) | \(\|e\|_{H^1}\) | \(\alpha_{H^1}\) |
|---|---|---|---|---|---|---|---|
| 1 | 10 | 2916 | 2953 | 2.36e-04 | - | 1.97e-02 | - |
| 1 | 11 | 3600 | 3641 | 1.92e-04 | 1.98 | 1.77e-02 | 0.98 |
| 1 | 12 | 4356 | 4401 | 1.59e-04 | 1.99 | 1.62e-02 | 0.98 |
| 1 | 13 | 5184 | 5233 | 1.33e-04 | 1.99 | 1.48e-02 | 0.98 |
| 2 | 10 | 2916 | 11737 | 8.04e-06 | - | 8.41e-04 | - |
| 2 | 11 | 3600 | 14481 | 5.89e-06 | 2.95 | 6.84e-04 | 1.96 |
| 2 | 12 | 4356 | 17513 | 4.44e-06 | 2.96 | 5.67e-04 | 1.96 |
| 2 | 13 | 5184 | 20833 | 3.43e-06 | 2.96 | 4.78e-04 | 1.97 |
| 3 | 10 | 2916 | 26353 | 2.63e-07 | - | 3.74e-05 | - |
| 3 | 11 | 3600 | 32521 | 1.74e-07 | 3.92 | 2.75e-05 | 2.93 |
| 3 | 12 | 4356 | 39337 | 1.20e-07 | 3.93 | 2.08e-05 | 2.94 |
| 3 | 13 | 5184 | 46801 | 8.49e-08 | 3.94 | 1.61e-05 | 2.95 |
The following notations were used in the table header:
The auxiliary H-field was calculated by feeding the numerically calculated reduced magnetic scalar potential, \(\Theta\), to an object of the type StaticScalarSolver::ProjectPSItoH. The figure below illustrates the calculated auxiliary H-field.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 1 | 10 | 2916 | 5868 | 1.97e-02 | - |
| 1 | 11 | 3600 | 7240 | 1.77e-02 | 0.98 |
| 1 | 12 | 4356 | 8756 | 1.62e-02 | 0.98 |
| 1 | 13 | 5184 | 10416 | 1.48e-02 | 0.98 |
| 2 | 10 | 2916 | 23400 | 8.41e-04 | - |
| 2 | 11 | 3600 | 28880 | 6.84e-04 | 1.96 |
| 2 | 12 | 4356 | 34936 | 5.67e-04 | 1.96 |
| 2 | 13 | 5184 | 41568 | 4.78e-04 | 1.97 |
| 3 | 10 | 2916 | 52596 | 3.74e-05 | - |
| 3 | 11 | 3600 | 64920 | 2.75e-05 | 2.93 |
| 3 | 12 | 4356 | 78540 | 2.08e-05 | 2.94 |
| 3 | 13 | 5184 | 93456 | 1.61e-05 | 2.95 |
The following notations were used in the table header:
The figure below illustrates a contour plot of a two-dimensional slice of the three-dimensional reduced magnetic scalar potential, \(\Theta\), the output of an object of the SolverSLDII type.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) | \(\|e\|_{H^1}\) | \(\alpha_{H^1}\) |
|---|---|---|---|---|---|---|---|
| 1 | 5 | 8505 | 8808 | 2.73e-03 | - | 7.37e-02 | - |
| 1 | 6 | 15851 | 16288 | 1.80e-03 | 2.01 | 5.98e-02 | 1.00 |
| 1 | 7 | 26533 | 27128 | 1.27e-03 | 2.02 | 5.03e-02 | 1.01 |
| 1 | 8 | 41175 | 41952 | 9.45e-04 | 2.03 | 4.34e-02 | 1.01 |
| 2 | 5 | 8505 | 69131 | 2.85e-04 | - | 1.14e-02 | - |
| 2 | 6 | 15851 | 128407 | 1.54e-04 | 2.98 | 7.57e-03 | 1.96 |
| 2 | 7 | 26533 | 214467 | 9.18e-05 | 3.00 | 5.39e-03 | 1.98 |
| 2 | 8 | 41175 | 332303 | 5.91e-05 | 3.01 | 4.02e-03 | 2.00 |
| 3 | 5 | 8505 | 232000 | 3.01e-05 | - | 1.66e-03 | - |
| 3 | 6 | 15851 | 431464 | 1.34e-05 | 3.89 | 9.12e-04 | 2.89 |
| 3 | 7 | 26533 | 721216 | 6.84e-06 | 3.94 | 5.51e-04 | 2.94 |
| 3 | 8 | 41175 | 1118104 | 3.82e-06 | 3.96 | 3.57e-04 | 2.96 |
The following notations were used in the table header:
The auxiliary H-field was calculated by feeding the numerically calculated three-dimensional reduced magnetic scalar potential, \(\Theta\), to an object of the type StaticScalarSolver::ProjectPSItoH. The figure below illustrates a two-dimensional slice of the calculated three-dimensional auxiliary H-field.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 1 | 5 | 8505 | 26060 | 7.37e-02 | - |
| 1 | 6 | 15851 | 48352 | 5.98e-02 | 1.00 |
| 1 | 7 | 26533 | 80700 | 5.03e-02 | 1.01 |
| 1 | 8 | 41175 | 124976 | 4.34e-02 | 1.01 |
| 2 | 5 | 8505 | 206182 | 1.14e-02 | - |
| 2 | 6 | 15851 | 383474 | 7.57e-03 | 1.96 |
| 2 | 7 | 26533 | 641022 | 5.39e-03 | 1.98 |
| 2 | 8 | 41175 | 993802 | 4.02e-03 | 2.00 |
| 3 | 5 | 8505 | 693456 | 1.66e-03 | - |
| 3 | 6 | 15851 | 1290684 | 9.12e-04 | 2.89 |
| 3 | 7 | 26533 | 2158560 | 5.51e-04 | 2.94 |
| 3 | 8 | 41175 | 3347628 | 3.57e-04 | 2.96 |
The following notations were used in the table header:
The convergence tables above suggest that the order of convergence of the \(L^2\) and \(H^1\) error norms of the reduced magnetic scalar potential, \(\Theta\), can be described by the following expressions:
\[ \alpha_{L^2} \approx p + 1 \]
and
\[ \alpha_{H^1} \approx p, \]
where p is the order of the interpolating Lagrange polynomials. This corresponds to the upper boundary of the expected convergence order. The convergence order computed for the \(\vec{H}\) field equals
\[ \alpha_{L^2} \approx p. \]
This is to be expected as \(\vec{H}\) is related to \(\Theta\) by a first-order derivative.
1.9.1