 |
Logbook
(07-04-2025)
Static problems
|
|
|
Logbook
(07-04-2025)
Static problems
|
|
Classes:
1) BatchABC
2) SolverABC
3) ExactSolutionABC_PHI
4) SettingsABC
Files:
1) abc/src/main.cpp
2) abc/include/solver.hpp
3) abc/src/solver.cpp
4) abc/include/exact_solution.hpp
5) abc/src/exact_solution.cpp
6) abc/src/static_scalar_input.cpp
7) abc/include/settings.hpp
Many electrostatic systems are not shielded by a grounded conducting enclosure. The electrostatic field induced by such systems extends to infinity. Consequently, the problem domains for modeling such systems must be unbounded. A straightforward application of the finite element method to an unbounded domain would require a mesh with an infinite amount of cells. The last, of course, is impossible. There exists a number of methods for overcoming this problem. Most of these methods require a modification to the standard finite element method. There are, however, two methods of dealing with unbounded domains that can be implemented in the framework of the standard finite element method. These methods are: the method of truncation and the method of the first-order asymptotic boundary condition (ABC).
The method of truncation is the most simple and straightforward method. One just needs to truncate the infinite domain up to a certain boundary and apply the homogeneous Dirichlet boundary condition,
\[ \Phi = 0, \]
or the homogeneous Neumann boundary condition,
\[ \epsilon \hat{n} \cdot \vec{\nabla} \Phi=0, \]
to the outer boundary of the truncated domain.
The ABC requires a truncation of the problem as well. In the case of the ABC the outer boundary of the truncated domain must have a shape of a circle in two dimensions and that of a sphere in three dimensions. Both, the circle and the sphere, must be centered at the origin. The following Robin boundary condition is applied to the outer boundary
\[ \epsilon \hat{n} \cdot \vec{\nabla} \Phi + \frac{\epsilon}{R_{\infty}} \Phi =0, \]
where \(R_{\infty}\) is the radius of the boundary.
In this numerical experiment we will test both methods, i.e., the method of truncation and ABC. We will do so by means of the method of exact solutions. That is, we will solve two textbook problems, two-dimensional and three-dimensional, with a help of the scalar solver on progressively refined meshes and calculate \(L^2\) and \(H^1\) error norms. We will apply the homogeneous Dirichlet, Neumann, and Robin (ABC) conditions to the outer boundary of the problem domain and compare the results.
The following textbook problem is solved in the two-dimensional version of the experiment. (Section 6.1.8. Two parallel cylindrical conductors)
The system consists of two identical infinitely long cylindrical conductors of radius \(R\). Both conductors are located on the \(xy\)-plane parallel to the \(z\)-axis. The distance between each conductor and the \(z\)-axis equals \(x_0\). The space around the conductors is empty. The conductor located along the line \(x=-x_0\) is at a potential of \(\Phi=-\Phi_0\) Volts. The conductor located along the line \(x=x_0\) is at a potential of \(\Phi=\Phi_0\) Volts. The figure below illustrates the system and the corresponding problem domain.
We would like to calculate the electrostatic potential.
As soon as the system exhibits a translation symmetry along the \(z\) axis, we reduce the problem domain to a two-dimensional planar domain. We remove both conductors from the problem domain and model them by setting up Dirichlet boundary conditions on the surfaces that delineate them. We truncate the space up to a circle of radius \(R_{\infty}\). Therefore, the problem domain consists of three circular boundaries and the empty space between them.
This problem can be solved by the analytical methods. The solution to this problem reads
\begin{equation} \Phi = \frac{\lambda'}{2} \ln\Bigg[\frac{(x+d)^2+y^2}{(x-d)^2+y^2}\Bigg], \end{equation}
\begin{equation} \vec{\nabla} \Phi = \lambda' \Bigg[\frac{(x+d)\hat{x}+y\hat{y}}{(x+d)^2+y^2}- \frac{(x-d)\hat{x}+y\hat{y}}{(x-d)^2+y^2} \Bigg], \end{equation}
where
\begin{equation} d=\sqrt{x_0^2-R^2}, \end{equation}
and
\begin{equation} \lambda' = \frac{ 1 }{\ln \bigg[\frac{x_0}{R}+\sqrt{\frac{x_0^2}{R^2}-1}\bigg]} \Phi_0. \end{equation}
The following textbook problem is solved in the three-dimensional version of the experiment. (Section 6.1.7. Isolated sphere in free space)
An isolated conducting sphere of radius \(a\) is located in an free unbounded space. The sphere is at potential of \(\Phi=\Phi_0[V]\).
We would like to calculate the electrostatic potential.
As the system exhibits a spherical symmetry, we can reduce the problem domain to a one-dimensional domain. Our intention, however, is to test the solver in three-dimensional space. Therefore, we will leave the problem domain as it is: a regular three-dimensional domain. We remove the conductor from the problem domain and model it by setting up Dirichlet boundary conditions on the surface that delineates it. We truncate the space up to a sphere of radius \(R_{\infty}\). Therefore, the problem domain consists of two concentric spheres and the empty space between them.
The solution to this problem reads
\begin{equation} \Phi = a \frac{1}{r} \Phi_0, \end{equation}
\begin{equation} \vec{\nabla} \Phi = -a\frac{1}{r^2}\Phi_0\hat{r}. \end{equation}
Which boundary value problem solved in this numerical experiment depends on the macro definition BC_TYPE__ set in CMakeLists.txt. If BC_TYPE__ is set to 1, the following boundary value problems are solved in the two- and three- dimensional versions of the experiment, respectively:
(2D)
\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big( \epsilon_0 \vec{\nabla} \Phi \big)= 0 & \text{in} & \Omega & \text{(i)},\\ \text{(e)} &\Phi = \Phi_0 & \text{on} & \Gamma_{D1} & \text{(ii)},\\ \text{(e)} &\Phi = 0 & \text{on} & \Gamma_{D2} & \text{(iii)},\\ \text{(e)} &\Phi =-\Phi_0 & \text{on} & \Gamma_{D3} & \text{(iv)}, \end{array} \end{equation}
and
(3D)
\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big( \epsilon_0 \vec{\nabla} \Phi \big)= 0 & \text{in} & \Omega & \text{(i)},\\ \text{(e)} &\Phi = \Phi_0 & \text{on} & \Gamma_{D1} & \text{(ii)},\\ \text{(e)} &\Phi = 0 & \text{on} & \Gamma_{D2} & \text{(iii)}. \end{array} \end{equation}
In the case BC_TYPE__=2, the following boundary value problems are solved in two- and three- dimensional versions of the experiment, respectively:
(2D)
\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big( \epsilon_0 \vec{\nabla} \Phi \big)= 0 & \text{in} & \Omega & \text{(i)},\\ \text{(e)} &\Phi = \Phi_0 & \text{on} & \Gamma_{D1} & \text{(ii)},\\ \text{(n)}&\epsilon_0\hat{n}\cdot\vec{\nabla}\Phi+ \big(\epsilon_0/R_{\infty}\big) \Phi = 0&\text{on}&\Gamma_{R1}&\text{(iii)},\\ \text{(e)} &\Phi = -\Phi_0 & \text{on} & \Gamma_{D3} & \text{(iv)}, \end{array} \end{equation}
and
(3D)
\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big( \epsilon_0 \vec{\nabla} \Phi \big)= 0 & \text{in} & \Omega & \text{(i)},\\ \text{(e)} &\Phi = \Phi_0 & \text{on} & \Gamma_{D1} & \text{(ii)},\\ \text{(n)}&\epsilon_0\hat{n}\cdot\vec{\nabla}\Phi+ \big(\epsilon_0/R_{\infty}\big) \Phi = 0&\text{on}&\Gamma_{R1}&\text{(iii)}, \end{array} \end{equation}
where \(R_{\infty}\) is the radius of the circle (in 2D) or a sphere (in 3D) that represent the outer boundary of the problem domain. In the case BC_TYPE__=0, no boundary condition is applied to the outer boundary of the problem domain. Recall, that applying no boundary is as good as applying the homogeneous Neumann boundary condition. In effect the following two boundary value problem are solved
(2D)
\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big( \epsilon_0 \vec{\nabla} \Phi \big)= 0 & \text{in} & \Omega & \text{(i)},\\ \text{(e)} &\Phi = \Phi_0 & \text{on} & \Gamma_{D1} & \text{(ii)},\\ \text{(n)}&\epsilon_0\hat{n}\cdot\vec{\nabla} \Phi = 0&\text{on}&\Gamma_{R1}&\text{(iii)},\\ \text{(e)} &\Phi = -\Phi_0 & \text{on} & \Gamma_{D3} & \text{(iv)}, \end{array} \end{equation}
and
(3D)
\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big( \epsilon_0 \vec{\nabla} \Phi \big)= 0 & \text{in} & \Omega & \text{(i)},\\ \text{(e)} &\Phi = \Phi_0 & \text{on} & \Gamma_{D1} & \text{(ii)},\\ \text{(n)}&\epsilon_0\hat{n}\cdot\vec{\nabla}\Phi = 0&\text{on}&\Gamma_{R1}&\text{(iii)}, \end{array} \end{equation}
if BC_TYPE__ is set to zero in CMakeLists.txt.
The three conditions are summarized in the table below.
| BC_TYPE__ | BC on conductors | Boundary IDs on conductors | BC at infinity | Boundary ID at infinity |
|---|---|---|---|---|
| 0 | Dirichlet | 1, 3 | Neumann | 0 |
| 1 | Dirichlet | 1, 3 | Dirichlet | 5 |
| 2 | Dirichlet | 1, 3 | ABC | 2 |
All the 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 build process generates three executable files for the two- dimensional version of the experiment: abc-ppc-dirichlet, abc-ppc-neumann, and abc-ppc-abc. The names of the files consist of three words each. The first word, abc, is the name of the experiment. The second word, ppc, is an acronym for "pair of parallel conductors". It codifies the configuration of the problem domain. The third word codifies the boundary condition applied to the outer boundary of the problem domain. That is, each executable file applies one of the following boundary condition on the outer boundary: the Dirichlet, Neumann, or ABC. Similarly, the build process generates three files for the three-dimensional version of the experiment: abc-shell-dirichlet, abc-shell-neumann, and abc-shell-abc. The word shell in the names of the files codifies the configuration of the problem domain. Other words in the names of the files and their meaning are exactly the same as in the two-dimensional version of the experiment.
To rebuild the executable files change into abc/build/Release directory and execute the following:
Then all executable files must be executed again. This can be done by changing into the abc/bin/Release directory and executing run-all script there,
This will generate various files in the abc/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 abc/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into abc/gmsh directory and executing the following:
This will generate a set of globally refined meshes in abc/gmsh/data.
The SettingsABC 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 vtu files next to the numerical solution. The last can be useful when debugging.
Inside each executable file the radius of the outer boundary is varied as the following
\[ R_{\infty} = 2 m R_{\text{mid}}, \text{ } \text{ } m = 1,2,3,4,5, \]
where
\[ R_{\text{mid}} = 1.5 (R+x_0). \]
The parameters \(R\) and \(x_0\) in the last equation are shown in the figure that depicts the two-dimensional configuration, see above. The error norms are calculated only within limited area (or volume in 3D) of the problem domains. The error norms outside this limited area (or volume in 3D) are set to zero. The last is done by the class StaticScalarSolver::Weight. The figure below gives an example of the \(||e||_{L^2}\) norms calculated for the two- and three- dimensional versions of the experiment under condition \(m=1\).
Although the radius of the outer boundary increases proportionally to \(m\), the size of the limited area (or volume in 3D) in which the error norms are calculated stays fixed. In the two- dimensional version of the experiment it is a square with a side of \( 4x_0\). In the three-dimensional version of the experiment it is a spherical shell with the inner radius of \(a\) and the outer radius of \(R_{\text{mid}} = 1.5 (R+x_0)\). The outer radius of the three- dimensional shell depends on the parameters of the two- dimensional domain, i.e., \(R\) and \(x_0\), so the dimensions of the domains are roughly of the same order of magnitude.
All simulations are repeated for three values of the degree of the interpolating polynomials, \(p=1\), \(2\), and \(3\). The degree of mesh refinement is varied as discussed on this page. The degree of refinement is codified by the parameter \(r=1\), \(2\), and \(3\). The radius of the outer boundary of the problem domain is varied as discussed above. Consequently, resultant error norms are functions of: the radius of the outer boundary (represented by m), the degree of the interpolating polynomials ( \(p\)), the degree of mesh refinement (represented by \(r\)), the type of boundary condition applied on the outer boundary (Dirichlet, Neumann, ABC). All the data is saved in the abc/bin/Release/Data directory.
The experiment was executed under the following conditions: \(x_0=0.1[m]\), \(R=0.05[m]\), \(a=0.1[m]\), \(\Phi_0=1[V]\).
The figure below illustrates the electrostatic potential calculated for \(m=1\).
The two plots below show the \(L^2\) error norms for the least computationally demanding case ( \(p=1\) and \(r=1\)), and for the most computationally demanding case ( \(p=3\) and \(r=3\)). The cases in between these two extremes can be found in here. These plots suggest that the first-order ABC outperforms the Dirichlet and Neumann boundary conditions. The superiority of the first-order ABC is more pronounced for higher degrees of the interpolating polynomials and on finer meshes.
Similar results has been obtained for the \(H^1\) error norms as the plots below illustrate. The rest of the plots can be found in here.
The figure below illustrates the electrostatic potential calculated for \(m=1\).
The two plots below illustrate \(L^2\) norms. Here again, the most and the least demanding cases with respect to the computational load are shown. The rest of the plots can be found in here.
Similar results has been obtained for the \(H^1\) error norms.
The rest of the plots can be found in here.
In all data sets that has been produced by this numerical experiment, the first-order ABC outperforms the method of truncation, i.e., the Dirichlet and Neumann boundary conditions. This is especially true if higher-order polynomials are used on finer meshes.
1.9.1 
