 |
Logbook
(07-04-2025)
Static problems
|
|
|
Logbook
(07-04-2025)
Static problems
|
|
Classes:
1) BatchRHOAXI
2) SolverRHOAXI
3) ExactSolutionRHOAXI_PHI
4) SettingsRHOAXI
Files:
1) rho-axi/src/main.cpp
2) rho-axi/include/solver.hpp
3) rho-axi/src/solver.cpp
4) rho-axi/include/exact_solution.hpp
5) rho-axi/src/exact_solution.cpp
6) rho-axi/src/static_scalar_input.cpp
7) rho-axi/include/settings.hpp
In the preceding numerical experiments that deal with axisymmetric problems the volume free-charge density, \(\rho_f\), was a continuous function of spatial coordinates. In this numerical experiment we will test if the StaticScalarSolver::Solver class template can solve axisymmetric problems in which the volume free-current density \(\rho_f\) is discontinuous. To do so, we will solve two problems of the Volume charge (rho/) numerical experiment on axisymmetric domains with a help of StaticScalarSolver::Solver class template and compare the results to closed-form analytical solutions.
We convert the first problem of the Volume charge (rho/) numerical experiment into an axisymmetric problem domain as shown in the figure below. 
We expect that the electrostatic field will have only radial component and, thus, will be aligned with the vertical segments of the boundary. Furthermore, we can show that the electrostatic field on the axis of rotation symmetry equals zero. To do so, we can draw a cylindrical Gaussian box coaxial with the \(z\) axis in the figure above. The net flux of the electric field through the surface of the box equals \(Q_f/\epsilon_0\), where \(Q_f\) is the total free charge enclosed by the box. If we drive the radius of the Gaussian box to zero, \(Q_f\) will be driven to zero implying that the net flux of the electrostatic field through the Gaussian box is driven to zero as well. If we take into account the cylindrical symmetry of the problem, we have no choice but to conclude that the electrostatic field in this configuration equals zero on the axis of rotation symmetry. Therefore, its tangential (as well as normal) component equals zero. Consequently, we can conclude that the tangential component of the electrostatic field equals zero on the three segments of the boundary labeled as \(\Gamma_{R1}\) in the figure above. This allows us to apply the homogeneous Neumann boundary condition,
\begin{equation} \epsilon\hat{n}\cdot\vec{\nabla}\Phi=0, \end{equation}
to these three segments of the boundary. To apply this boundary condition, we mark the three segments of the boundary, \(\Gamma_{R1}\), with ID=0 in the rho-axi/gmsh/cylinder.geo file. The program will read in the zero boundary ID's, and will do nothing with respect to application of the boundary conditions to these three segments of the boundary. The homogeneous Neumann boundary condition will be applied implicitly as it is encoded in the first term of the functional. We apply the Dirichlet boundary condition \(\Phi=0\) on the fourth segment of the boundary, \(\Gamma_{D1}\) in the figure above, as dictated by the formulation of the problem. To do so, we mark the boundary segment \(\Gamma_{D1}\) with ID=1 in rho-axi/gmsh/cylinder.geo file. The program will read in this boundary ID and will apply the Dirichlet boundary condition.
The second problem of the Volume charge rho/ numerical experiment exhibits spherical symmetry. Despite the spherical symmetry we will use the cylindrical coordinate system to derive an axisymmetric domain (see Section 2.3. Two-dimensional problems). We derive an axisymmetric problem domain as illustrated in the figure below. 
We expect that the electrostatic field will be aligned with the vertical segment of the boundary, \(\Gamma_{R1}\) in the figure above. For this reason we will apply the homogeneous Neumann boundary condition,
\begin{equation} \epsilon\hat{n}\cdot\vec{\nabla}\Phi=0, \end{equation}
to this segment. We mark this boundary segment with ID=0 in the rho-axi/gmsh/sphere.geo file. The program will read in the zero boundary ID, and will do nothing with respect to application of the boundary conditions to this segment of the boundary. The homogeneous Neumann boundary condition will be applied implicitly as it is encoded in the first term of the functional. We will apply the Dirichlet boundary condition on the curved segment of the boundary, i.e., \(\Phi = 0[V]\) on \(\Gamma_{D1}\) as dictated by the formulation of the problem. The program will read in the boundary ID and will apply the homogeneous Dirichlet boundary condition.
This numerical experiment solves the following boundary value problem.
\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla}\cdot\big(\epsilon_0\vec{\nabla}\Phi\big)=\rho_f & \text{in} & \Omega & \text{(i)},\\ \text{(n)} & \epsilon \hat{n} \cdot \vec{\nabla} \Phi = 0 & \text{on} & \Gamma_{R1} & \text{(ii)},\\ \text{(e)} &\Phi = 0 & \text{on} & \Gamma_{D1} & \text{(iii)}. \end{array} \end{equation}
This boundary value problem is a special case of the general static scalar boundary value problem.
This experiment is implemented in accordance with the base code structure. The build process generates two executable files: rho-axi-cylinder and rho-axi-sphere. That is, one executable file for the first problem above and one for the second. To rebuild them change into rho-axi/build/Release directory and execute the following:
Then all executable files must be executed again. This can be done by changing into the rho-axi/bin/Release directory and executing run-all script there,
This will generate various files in the rho-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.
Note that executable files require a set of meshes to be present in the rho-axi/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into rho-axi/gmsh directory and executing the following:
This will generate a set of globally refined meshes in rho-axi/gmsh/data.
The SettingsRHOAXI 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 program invokes first- and second-order mapping when solving the first and the second problems, respectively.
Both problems above were solved numerically with the following settings (see SettingsRHOAXI): \(\rho_f=\epsilon_0\), \(a=0.5[m]\) and \(b=1.0[m]\).
The figures below illustrate the numerically calculated electrostatic potential. The plot on the right allows a comparison between the results obtained on the two-dimensional planar domain (rho/ numerical experiment) and on the two-dimensional axisymmetric domain (rho-axi/ numerical experiment).
The corresponding convergence table is shown below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) | \(\|e\|_{H^1}\) | \(\alpha_{H^1}\) |
|---|---|---|---|---|---|---|---|
| 1 | 7 | 72 | 91 | 1.72e-04 | - | 7.00e-03 | - |
| 1 | 8 | 98 | 120 | 1.27e-04 | 1.96 | 5.98e-03 | 1.02 |
| 1 | 9 | 128 | 153 | 9.79e-05 | 1.96 | 5.22e-03 | 1.02 |
| 1 | 10 | 162 | 190 | 7.77e-05 | 1.97 | 4.63e-03 | 1.02 |
| 2 | 7 | 72 | 325 | 1.45e-06 | - | 1.13e-04 | - |
| 2 | 8 | 98 | 435 | 9.13e-07 | 2.98 | 8.30e-05 | 1.98 |
| 2 | 9 | 128 | 561 | 6.13e-07 | 2.98 | 6.37e-05 | 1.98 |
| 2 | 10 | 162 | 703 | 4.31e-07 | 2.99 | 5.04e-05 | 1.99 |
| 3 | 7 | 72 | 703 | 3.52e-08 | - | 4.02e-06 | - |
| 3 | 8 | 98 | 946 | 1.91e-08 | 3.96 | 2.55e-06 | 2.96 |
| 3 | 9 | 128 | 1225 | 1.13e-08 | 3.97 | 1.71e-06 | 2.97 |
| 3 | 10 | 162 | 1540 | 7.06e-09 | 3.97 | 1.21e-06 | 2.98 |
The following notations were used in the table header:
The figures below illustrate the numerically calculated electrostatic potential. The plot on the right allows a comparison between the results obtained in the three-dimensional domain (rho/ numerical experiment) and in the two-dimensional axisymmetric domain (rho-axi/ numerical experiment).
The corresponding convergence table is shown below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) | \(\|e\|_{H^1}\) | \(\alpha_{H^1}\) |
|---|---|---|---|---|---|---|---|
| 1 | 7 | 360 | 391 | 1.90e-04 | - | 7.81e-03 | - |
| 1 | 8 | 490 | 526 | 1.40e-04 | 1.97 | 6.71e-03 | 0.98 |
| 1 | 9 | 640 | 681 | 1.08e-04 | 1.97 | 5.89e-03 | 0.98 |
| 1 | 10 | 810 | 856 | 8.53e-05 | 1.98 | 5.24e-03 | 0.99 |
| 2 | 7 | 360 | 1501 | 4.69e-06 | - | 3.65e-04 | - |
| 2 | 8 | 490 | 2031 | 2.97e-06 | 2.97 | 2.70e-04 | 1.97 |
| 2 | 9 | 640 | 2641 | 1.99e-06 | 2.98 | 2.07e-04 | 1.98 |
| 2 | 10 | 810 | 3331 | 1.40e-06 | 2.98 | 1.64e-04 | 1.98 |
| 3 | 7 | 360 | 3331 | 1.62e-07 | - | 1.78e-05 | - |
| 3 | 8 | 490 | 4516 | 8.81e-08 | 3.95 | 1.13e-05 | 2.95 |
| 3 | 9 | 640 | 5881 | 5.19e-08 | 3.96 | 7.61e-06 | 2.96 |
| 3 | 10 | 810 | 7426 | 3.25e-08 | 3.97 | 5.36e-06 | 2.97 |
The following notations were used in the table header:
1.9.1