 |
Logbook
(07-04-2025)
Static problems
|
|
|
Logbook
(07-04-2025)
Static problems
|
|
Classes:
1) BatchMMSAXI
2) SolverMMSAXI
3) ExactSolutionMMSAXI_PHI
4) ExactSolutionMMSAXI_E
5) ExactSolutionMMSAXI_D
6) SettingsMMSAXI
Files:
1) mms-axi/src/main.cpp
2) mms-axi/include/solver.hpp
3) mms-axi/src/solver.cpp
4) mms-axi/include/exact_solution.hpp
5) mms-axi/src/exact_solution.cpp
6) mms-axi/src/static_scalar_input.cpp
7) mms-axi/include/settings.hpp
By default, all problems in electromagnetics are three-dimensional. As three-dimensional problems require more computer resources than two-dimensional problems, it is always beneficial to recast a three-dimensional problem to a two-dimensional problem. The last is possible if the initial three-dimensional problem exhibits a symmetry. We can distinguish two kinds of symmetries with this respect: the translation symmetry and the rotation symmetry.
The translation symmetry arises in problems that describe straight, infinitely long cable-like structures with a cross section that does not change along the structure. Any cross section made by a plane perpendicular to the axis of the structure presents the same picture. We can choose one cross section and use it as a two-dimensional problem domain. Such domains are called two-dimensional planar problem domains. The two-dimensional version of the Floating conductor (flc/) numerical experiment exemplifies the application of the two-dimensional planar domain.
The axisymmetric symmetry arises in problem domains that can be generated by rotating a two-dimensional structure around an axis. In such cases, we can reduce the dimensionality of the problem by choosing this two-dimensional structure as an equivalent two-dimensional problem domain. Such domains are called two-dimensional axisymmetric problem domains.
When we program a two- dimensional problem, deal.II automatically assumes that the problem domain is planar two-dimensional. Strictly speaking, deal.II makes no assumptions leaving this choice to the user. It is convenient, however, to think that it does. If we wish to implement an axisymmetric problem domain, we need to do an extra effort:
In this numerical experiment we will use the method of manufactured solutions to test if the code implemented by StaticScalarSolver::Solver, StaticScalarSolver::ProjectPHItoE, and StaticScalarSolver::ProjectPHItoD class templates can solve axisymmetric problems. In the three-dimensional version of this experiment we will solve a problem in which the manufactured solution and the problem domain exhibit the cylindrical symmetry. In the two-dimensional version of the experiment we will solve the same problem in an equivalent two-dimensional axisymmetric problem domain. In this numerical experiment we presume that the permittivity is a smooth function of spatial coordinates. We will also set the surface free-charge density, \(\kappa_f\), to zero.
We proceed in the cylindrical coordinate system \((r,\phi,z)\) and assume (somewhat arbitrary) that the manufactured solution reads
\[ \Phi_M(r, z) = \cos{(kr)} + \cos{(kz)}. \]
Next, we assume (somewhat arbitrary) that the coefficients on the left-hand side of the boundary value problem have the following form
\[ \epsilon = \epsilon_0 (r z^2 + 1), \]
\[ \gamma = \epsilon_0 \big( \sqrt{r^2 + z^2} + 2 \big). \]
Then the gradient of the manufactured solution and the right-hand side of the partial differential equation can be calculated by differentiating the manufactured solution:
\[ \vec{\nabla} \Phi_M(r, z) = - k \sin{(kr)} \hat{r} - k \sin{(kz)} \hat{z}, \]
and
\begin{equation} \rho_M = -\vec{\nabla} \cdot \bigg(\epsilon \vec{\nabla} \Phi_M \bigg)= \epsilon_0 k \bigg\{ \bigg( 2 z^2 + \frac{1}{r} \bigg) \sin{(kr)} + k (r z^2 + 1) \cos(kr) + 2 r z \sin(kz) + k (r z^2 + 1) \cos(kz) \bigg\}. \end{equation}
The manufactured electric field and manufactured displacement are evaluated as
\[ \vec{E}_M = - \vec{\nabla} \Phi_M = k \sin{(kr)} \hat{r} + k \sin{(kz)} \hat{z} \]
and
\[ \vec{D}_M = - \epsilon \vec{\nabla} \Phi_M = \epsilon_0 \big(r z^2 + 1\big) \big( k \sin{(kr)} \hat{r} + k \sin{(kz)} \hat{z} \big). \]
For a problem domain we choose a cylinder centered at the origin as shown in the figure below.
We apply the Dirichlet boundary condition on the two flat surfaces of the cylinder and the Robin boundary condition on the curved surface. We label the boundaries of the cylinder in accordance with the boundary ID convention. The manufactured solutions above are expressed in the cylindrical coordinate system. In three-dimensional version of the experiment they must be expressed in the Cartesian coordinate system. We recall that
\[ r = \sqrt{x^2 + y^2} \]
and
\[ \hat{r} = \cos{\phi} \hat{x} + \sin{\phi} \hat{y}. \]
Therefore, in the Cartesian coordinate system
\[ \Phi_M(x, y, z) = \cos{\bigg(k\sqrt{x^2+y^2}\bigg)} + \cos{(kz)}, \]
\[ \vec{\nabla} \Phi_M = -k \cos{\phi} \sin{\bigg(k\sqrt{x^2 + y^2}\bigg)} \hat{x} -k \sin{\phi} \sin{\bigg(k\sqrt{x^2 + y^2}\bigg)} \hat{y} -k \sin{(kz)} \hat{z}, \]
\[ \vec{E}_M = k \cos{\phi} \sin{\bigg(k\sqrt{x^2 + y^2}\bigg)} \hat{x} + k \sin{\phi} \sin{\bigg(k\sqrt{x^2 + y^2}\bigg)} \hat{y} + k \sin{(kz)} \hat{z}, \]
and
\[ \vec{D}_M = \epsilon_0 \big( z^2 \sqrt{x^2 + y^2} + 1\big) \bigg[ k \cos{\phi} \sin{\bigg(k\sqrt{x^2 + y^2}\bigg)} \hat{x} + k \sin{\phi} \sin{\bigg(k\sqrt{x^2 + y^2}\bigg)} \hat{y} + k \sin{(kz)} \hat{z} \Bigg]. \]
,
In the two-dimensional version of the experiment the manufactured solution remains the same. The three-dimensional problem domain and the manufactured solution described above exhibit rotation symmetry with the \(z\) axis being the axis of symmetry. Consequently, we can choose an \( rz \) plane for the problem domain as shown in the figure below.
Similarly to the three-dimensional case, we apply the Dirichlet boundary conditions to the two cross sections of the flat surfaces of the cylinder and the Robin boundary condition to the cross section of the curved boundary of the cylinder. The homogeneous Neumann boundary condition is the appropriate boundary condition on the axis of the rotation symmetry. Note that the manufactured electric field on the axis of symmetry has no radial component,
\[ \begin{array}{lcr} \vec{E}_M = k \sin{(kz)} \hat{z} & \text{if} & r = 0. \end{array} \]
Thus, the homogeneous Neumann boundary condition is a reasonable choice. Recall that the homogeneous Neumann boundary condition can be used on the axis of symmetry in all axisymmetric problems (the electric field is always aligned with the axis of rotation symmetry). We label the boundaries of the cylinder in accordance with the boundary ID convention.
The two-dimensional version of this experiment solves the following boundary value problem by finite element method:
\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big( \epsilon \vec{\nabla} \Phi \big)= \rho_M & \text{in} & \Omega & \text{(i)},\\ \text{(e)} &\Phi = \eta & \text{on} & \Gamma_{D1} & \text{(ii)},\\ \text{(n)}&\epsilon\hat{n}\cdot\vec{\nabla}\Phi+\gamma\Phi=\sigma &\text{on}&\Gamma_{R1}&\text{(iii)},\\ \text{(n)}&\epsilon\hat{n}\cdot\vec{\nabla}\Phi=0&\text{on}&\Gamma_{R2}&\text{(iv)}.\\ \end{array} \end{equation}
The three-dimensional version of this experiment solves the following boundary value problem by finite element method:
\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big( \epsilon \vec{\nabla} \Phi \big)= \rho_M & \text{in} & \Omega & \text{(i)},\\ \text{(e)} &\Phi = \eta & \text{on} & \Gamma_{D1} & \text{(ii)},\\ \text{(n)}&\epsilon\hat{n}\cdot\vec{\nabla}\Phi+\gamma\Phi=\sigma &\text{on}&\Gamma_{R1}&\text{(iii)}.\\ \end{array} \end{equation}
The right-hand sides of the Dirichlet and the Robin boundary conditions are implemented by evaluating the following expressions on the corresponding boundaries:
\[ \eta = \Phi_M, \]
\[ \sigma = \epsilon \hat{n} \cdot \vec{\nabla} \Phi_M + \gamma \Phi_M. \]
This boundary value problems above are special cases of the general static scalar boundary value problem.
This experiment is implemented in accordance with the base code structure. The mms-axi/ program is, essentially, a modified version of the mms/ program. The build process generates two executable files: mms-axi-cylinder_2d and mms-axi-cylinder_3d. That is, as their names suggest, one file for the two-dimensional version of the experiment and one file for the three-dimensional version of the experiment. To rebuild them change into mms-axi/build/Release directory and execute the following:
Then all executable files must be executed again. This can be done by changing into the mms-axi/bin/Release directory and executing run-all script there,
This will generate various files in the mms-axi/bin/Release/Data directory. Among the generated files there are vtk files that can be viewed with a help of VisIt software package of the Lawrence Livermore National Laboratory. The data files in tex and txt format contain the convergence tables. Alternatively, the two executable files can be run one-by one. In this case, however, there will be no files containing aggregated tables in the mmsi-axi/bin/Release/Data/ directory.
Note that executable files require a set of meshes to be present in the mms-axi/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into mms-axi/gmsh directory and executing the following:
This will generate a set of refined meshes in mms-axi/gmsh/data.
The SettingsMMSAXI 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 (manufactured) solutions into the vtk files next to the numerical solution. The last can be useful when debugging.
All the simulations described below have been done with a setting \(k=\pi\).
The two figures below illustrate the calculated electric potentials and list the corresponding convergence data.
mms-axi-cylinder-3d 
mms-axi-cylinder-2d 
The following notations are used in the headers of the tables that describe the convergence of the electric potential, \(\Phi\):
The three-dimensional plot above can be derived by rotating the two-dimensional plot about the vertical axis. Therefore, we can conclude that the two-dimensional plot looks exactly as we would expect from the solution of the axisymmetric problem. The convergence orders in the two-dimensional axisymmetric case fit the theoretical predictions:
\[ \alpha_{L^2} \approx p + 1 \]
and
\[ \alpha_{H^1} \approx p. \]
The two figures below illustrate the calculated electric field and list the corresponding convergence data.
mms-axi-cylinder-3d 
mms-axi-cylinder-2d 
The following notations are used in the headers of the tables that describe the convergence of the electric field, \( \vec{E} \):
The two figures below illustrate the calculated displacement and list the corresponding convergence data.
mms-axi-cylinder-3d 
mms-axi-cylinder-2d 
The following notations are used in the headers of the tables that describe the convergence of the displacement, \( \vec{D} \):
1.9.1