 |
Logbook
(07-04-2025)
Static problems
|
|
|
Logbook
(07-04-2025)
Static problems
|
|
Classes:
1) BatchSSOLIAXI
2) SolverSSOLIAXI
3) ExactSolutionSSOLIAXI_A
4) ExactSolutionSSOLIAXI_B
5) SettingsSSOLIAXI
Files:
1) ssol-i-axi/src/main.cpp
2) ssol-i-axi/include/solver.hpp
3) ssol-i-axi/src/solver.cpp
4) ssol-i-axi/include/exact_solution.hpp
5) ssol-i-axi/src/exact_solution.cpp
6) ssol-i-axi/src/static_scalar_input.cpp
7) ssol-i-axi/include/settings.hpp
In the most general case a three-dimensional problem formulated in terms of the vector magnetic potential, \(\vec{A}\), is described by the curl-curl partial differential equation,
\[ \vec{\nabla} \times \bigg( \frac{1}{\mu} \vec{\nabla} \times \vec{A} \bigg) = \vec{J}_f. \]
If, however, a problem exhibits a rotation symmetry and the magnetic vector potential is expected to have the following form in the cylindrical coordinate system:
\[ \vec{A}(r,z) = 0\hat{r} + A(r,z) \hat{\phi} + 0 \hat{z}, \]
the curl-curl equation can be replaced by the following div-grad equation:
\[ - \vec{\nabla} \cdot \bigg(\frac{1}{\mu'} \vec{\nabla} A' \bigg) = J_f, \]
where
\[ \mu' = r \mu \]
and
\[ A' = r A, \]
see Section 3.4.2 Scalar axisymmetric problem.
This numerical experiment illustrates the application of the StaticScalarSolver::Solver to a problem formulated in terms of the scaled magnetic vector potential, \(A'\). To this end, we will obtain numerical solutions to a problem on a set of progressively refined meshes and compare them to the closed-form analytical solution.
In this numerical experiment the problem of the ssol-i numerical experiment is reused. This time, however, we do not use the three-dimensional problem domain and construct a two-dimensional axisymmetric domain as shown in the figure below.
We would like to calculate the scaled magnetic vector potential, \(A'\), and scaled magnetic field, \(\vec{B}'\).
The program operates in two modes, "Dirichlet" and "Exact". In the "Dirichlet" mode the homogeneous Dirichlet boundary condition,
\[ A' = 0, \]
is applied to both segments of the boundary, \(\Gamma_{D1}\) and \(\Gamma_{D2}\). In the "Exact" mode the homogeneous Dirichlet boundary condition on \(\Gamma_{D2}\) is replaced with
\[ A' = \eta, \]
where \(\eta\) is computed by evaluating the closed-form analytical expression of the solution given below. The "Exact" mode stitches off the simulation errors associated with truncation of the problem domain. This allows observation of the convergence rates free of truncation errors. The program can be toggled between the two modes by modifying Setting::type_of_bc.
To derive a closed-form analytical expression for the scaled magnetic field we make the following substitution into the analytical solution to the problem of ssol-i numerical experiment:
\[ r \rightarrow s \]
and multiply the result by the distance to the axis of rotation symmetry, \(r\),
\[ \vec{B}' = r \vec{B}(s,\theta) = \left\{ \begin{aligned} &r\dfrac{2}{3}\mu_0 K_0 a\big(\cos{\theta}\hat{s}-\sin{\theta}\hat{\theta}\big)&&\text{if }&& s \le a\\ &r\dfrac{2}{3}\mu_0 K_0 \frac{a^4}{s^3}\big(\cos{\theta}\hat{s}+\frac{1}{2}\sin{\theta}\hat{\theta} \big)&&\text{if }&& s\ge a. \\ \end{aligned} \right. \]
In the last equation \(a\) is the radius of the coil, \(s\) is the distance to the origin, and \(\theta\) is the polar angle. That is, here \(s\) and \(\theta\) are spherical coordinates while \(r\) and \(z\) are a cylindrical coordinates, see the figure above.
This time around the magnetic vector potential is gauged implicitly by the rotation symmetry of the problem. For this reason, we can expect a unique solution in terms of \(A'\). Consequently, we can compute the \(L^2\) and \(H^1\) error norms. To do so, we need close-form expressions for \(A'\) and its gradient. Note that in the current context by the gradient we mean the following differential operator:
\[ \vec{\nabla} = \dfrac{\partial}{\partial r} \hat{r} + \dfrac{\partial}{\partial z} \hat{z}, \]
To derive an expression for the equation for the scaled magnetic vector potential we make the following substitution in the expression for \(\vec{A}\) presented in Section 2.4.6:
\[ r \rightarrow s \]
and multiply the result by \(r\),
\[ A' = r A(s,\theta)=\left\{ \begin{aligned} &r\dfrac{1}{3}\mu_0 K_0 a s \sin{\theta} &&\text{if } && s \le a\\ &r\dfrac{1}{3}\mu_0 K_0 a^4 \dfrac{1}{s^2}\sin{\theta}&&\text{if }&& s\ge a. \end{aligned} \right. \]
Next, we note that
\[ \sin(\theta) = \dfrac{r}{s} \]
and, thus,
\[ A' = \left\{ \begin{aligned} &\dfrac{1}{3}\mu_0 K_0 a r^2 &&\text{if } && s \le a\\ &\dfrac{1}{3}\mu_0 K_0 a^4 \dfrac{r^2}{s^3}&&\text{if }&& s\ge a. \end{aligned} \right. \]
Next, we note that
\[ \dfrac{\partial}{\partial r} \dfrac{r^2}{s^3} = \dfrac{\partial}{\partial r} r^2(r^2+z^2)^{-\frac{3}{2}} = 2r(r^2+z^2)^{-\frac{3}{2}} -\dfrac{3}{2} r^2(r^2+z^2)^{-\frac{5}{2}}2r= 2r(r^2+z^2)^{-\frac{3}{2}} -3 r^3(r^2+z^2)^{-\frac{5}{2}} = \dfrac{2r}{s^3} - \dfrac{3r^3}{s^5} \]
and
\[ \dfrac{\partial}{\partial z} \dfrac{r^2}{s^3} = \dfrac{\partial}{\partial z} r^2(r^2+z^2)^{-\frac{3}{2}} = -\dfrac{3}{2} r^2(r^2+z^2)^{-\frac{5}{2}}2z = - \dfrac{3r^2 z}{s^5}. \]
Therefore,
\[ \vec{\nabla}A' = \left\{ \begin{aligned} &\dfrac{2}{3}\mu_0 K_0 a r \hat{r} &&\text{if } && s \le a\\ &\dfrac{1}{3}\mu_0 K_0 a^4 \bigg[\bigg(\dfrac{2r}{s^3} - \dfrac{3r^3}{s^5}\bigg)\hat{r} - \dfrac{3r^2 z}{s^5}\hat{z} \bigg] &&\text{if }&& s\ge a. \end{aligned} \right. \]
The figure below illustrates the mesh.
The radius \(b\) of the curved segment of the boundary of the mesh is relatively small. This is adequate in the "Exact" mode of operation. One needs to increase the radius if "Dirichlet" mode is used.
The following boundary value problem is solved in this numerical experiment:
\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \bigg(\dfrac{1}{\mu'_{0}} \vec{\nabla} A' \bigg)= 0 & \text{in} & \Omega & \text{(i)},\\ \text{(e)} & A' = 0 & \text{on} & \Gamma_{D1} & \text{(ii)},\\ \text{(e)} & A' = \eta & \text{on} & \Gamma_{D2} & \text{(iii)},\\ \text{(e)} & A_{+} = A_{-} & \text{on} & \Gamma_{I1} & \text{(iv)},\\ \text{(n)}&\dfrac{1}{\mu'_0}\hat{n}\cdot\bigg(\vec{\nabla} A'_{+}\bigg) -\dfrac{1}{\mu'_0}\hat{n}\cdot\bigg(\vec{\nabla}A'_{-}\bigg)= -K_f&\text{on}&\Gamma_{I1}&\text{(v)}, \end{array} \end{equation}
where
\begin{equation} \mu'_0 = r \mu_0. \end{equation}
The parameter \(\eta\) equals zero in the "Dirichlet" mode. In the "Exact" mode it is computed by evaluating the exact expression for \(A'\) given above.
First, a class template (SolverSSOLIAXI) derived from StaticScalarSolver::Solver is used to obtain the numerical solution, \(A'\). After that, the numerical solution is fed to objects derived from StaticScalarSolver::ProjectAphiToBrz class template to calculate the scaled magnetic field \(\vec{B}'\) as
\[ \vec{B}' = - \vec{\nabla}\overset{V}{\times} A'. \]
We will not compute the scaled auxiliary vector field \(\vec{H}'\). In this particular problem \(\vec{H}'\) cannot be modeled by the Nedelec finite elements as the necessary condition \(K_f =0\) does not hold.
The SSOLIAXI experiment is implemented in accordance with the base code structure. The build process generates one executable file: ssol-i-axi-circle. To rebuild it change into ssol-i-axi/build/Release directory and execute the following:
Then the executable file must be executed again. This can be done by changing into the ssol-i-axi/bin/Release directory and executing ssol-i-axi-circle,
This will generate various files in the ssol-i-axi/bin/Release/Data directory. Among the generated files there are vtu files that can be viewed with a help of ParaView software package of the Kitware, Inc. The data files in tex and txt format contain the convergence tables. The ssol-i-axi/bin/Release directory also contains circle.gpi file. It can be used to visualize the calculated potential along the \(r\) axis.
Note that executable files require a set of meshes to be present in the ssol-i-axi/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into ssol-i-axi/gmsh directory and executing the following:
This will generate a set of globally refined meshes in ssol-i-axi/gmsh/data.
The SettingsSSOLIAXI 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 experiment was executed under the following conditions: \(a = 0.5[m]\), \(b=1.0[m]\), \(K_0=1.0[A \cdot m^{-1}]\).
All convergence tables were generated in the "Exact" mode, see above.
The figure below illustrates a plot of the scaled magnetic vector potential, \(A'\), the output of an object of the SolverSSOLIAXI type.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}/\mu_0\) | \(\alpha_{L^2}\) | \(\|e\|_{H^1}/\mu_0\) | \(\alpha_{H^1}\) |
|---|---|---|---|---|---|---|---|
| 1 | 15 | 1960 | 2031 | 2.96e-05 | - | 1.78e-03 | - |
| 1 | 16 | 2250 | 2326 | 2.63e-05 | 1.70 | 1.66e-03 | 1.01 |
| 1 | 17 | 2560 | 2641 | 2.35e-05 | 1.71 | 1.56e-03 | 1.01 |
| 1 | 18 | 2890 | 2976 | 2.12e-05 | 1.71 | 1.47e-03 | 1.01 |
| 2 | 15 | 1960 | 7981 | 1.33e-07 | - | 2.38e-05 | - |
| 2 | 16 | 2250 | 9151 | 1.08e-07 | 3.00 | 2.07e-05 | 2.00 |
| 2 | 17 | 2560 | 10401 | 8.90e-08 | 3.00 | 1.82e-05 | 2.00 |
| 2 | 18 | 2890 | 11731 | 7.42e-08 | 3.00 | 1.61e-05 | 2.00 |
| 3 | 15 | 1960 | 17851 | 1.73e-09 | - | 4.40e-07 | - |
| 3 | 16 | 2250 | 20476 | 1.31e-09 | 3.98 | 3.58e-07 | 2.99 |
| 3 | 17 | 2560 | 23281 | 1.03e-09 | 3.85 | 2.95e-07 | 2.99 |
| 3 | 18 | 2890 | 26266 | 8.11e-10 | 3.86 | 2.46e-07 | 2.99 |
The following notations were used in the table header:
The figure below illustrates a plot of the calculated scaled magnetic field, \(\vec{B}'\), the output of an object of the StaticScalarSolver::ProjectArToBrz type.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}/\mu_0\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 1 | 15 | 1960 | 3990 | 1.78e-03 | - |
| 1 | 16 | 2250 | 4575 | 1.66e-03 | 1.01 |
| 1 | 17 | 2560 | 5200 | 1.56e-03 | 1.01 |
| 1 | 18 | 2890 | 5865 | 1.47e-03 | 1.01 |
| 2 | 15 | 1960 | 15820 | 2.38e-05 | - |
| 2 | 16 | 2250 | 18150 | 2.07e-05 | 2.00 |
| 2 | 17 | 2560 | 20640 | 1.82e-05 | 2.00 |
| 2 | 18 | 2890 | 23290 | 1.61e-05 | 2.00 |
| 3 | 15 | 1960 | 35490 | 4.40e-07 | - |
| 3 | 16 | 2250 | 40725 | 3.58e-07 | 2.99 |
| 3 | 17 | 2560 | 46320 | 2.95e-07 | 2.99 |
| 3 | 18 | 2890 | 52275 | 2.46e-07 | 2.99 |
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 scaled magnetic vector potential, \(A'\), 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 scaled magnetic field, \(\vec{B}'\) equals
\[ \alpha_{L^2} \approx p. \]
This is to be expected as \(\vec{B}'\) is related to \(A'\) by a first-order derivative.
1.9.1