 |
Logbook
(07-04-2025)
Static problems
|
|
|
Logbook
(07-04-2025)
Static problems
|
|
Classes:
1) BatchINT
2) SolverINT
3) ExactSolutionINT_PHI
4) ExactSolutionINT_E
5) ExactSolutionINT_D
6) SettingsINT
Files:
1) int/src/main.cpp
2) int/include/solver.hpp
3) int/src/solver.cpp
4) int/include/exact_solution.hpp
5) int/src/exact_solution.cpp
6) int/src/static_scalar_input.cpp
7) int/include/settings.hpp
In the Method of manufactured solutions (mms/) numerical experiment we have tested the class templates
assuming that the permittivity \(\epsilon\) is a smooth function of spatial coordinates. In this numerical experiment we will test if these class templates can handle problems in which the permittivity is discontinuous on the interfaces between dissimilar dielectric materials. We will do so by means of the method of exact solutions. That is, we will solve two textbook problems that contain interfaces between dissimilar dielectric materials on progressively refined meshes and calculate error norms. As a consequence of the mesh refinement the error norms should converge to zero.
The following textbook problem is solved in the two-dimensional version of the experiment. (Section 6.1.3. Cylindrical interface between dielectrics)
The system consists of two conductors and two layers of dissimilar dielectric materials. The conductors and the two layers of dielectric materials are shaped as infinitely long cylindrical tubes. All the tubes are coaxial. The two dielectric materials fill all the space between the conductors forming one interface. The interface is shaped as an infinitely long cylinder coaxial with the tubes. The inner conductor is at a potential of \(\Phi=\Phi_0\) Volts. The outer conductor is grounded. The figure below illustrates the system and the corresponding problem domain.
We would like to calculate the electric potential, electric field, and displacement.
As soon as the system exhibits translation symmetry, we can reduce the problem domain to a two-dimensional planar domain. We remove the conductors from the problem domain and model them by setting up Dirichlet boundary conditions on the surfaces that delineate them. Consequently, the problem domain consists of two ring-shaped dissimilar dielectric materials with Dirichlet conditions applied to the inner and outer boundaries.
This problem can be solved separately for each subdomain by the analytical methods. The solution reads
\begin{equation} \begin{aligned} &\Phi_1 = -\alpha\epsilon_2 \bigg[ \ln(r/d) - \frac{\epsilon_1}{\epsilon_2}\big[\ln(b/d)\big] \bigg]\Phi_0, \\ &\Phi_2 = -\alpha\epsilon_1 \bigg[ \ln(r/b) \bigg]\Phi_0, \end{aligned} \end{equation}
\begin{equation} \begin{aligned} &\vec{\nabla} \Phi_1 = -\alpha\epsilon_2 \frac{1}{r}\Phi_0\hat{r}, \\ &\vec{\nabla} \Phi_2 = -\alpha\epsilon_1 \frac{1}{r}\Phi_0\hat{r}, \end{aligned} \end{equation}
\begin{equation} \begin{aligned} & \vec{E}_1 = - \vec{\nabla} \Phi_1, \\ & \vec{E}_2 = - \vec{\nabla} \Phi_2, \end{aligned} \end{equation}
\begin{equation} \begin{aligned} & \vec{D}_1 = - \epsilon_1 \vec{\nabla} \Phi_1, \\ & \vec{D}_2 = - \epsilon_2 \vec{\nabla} \Phi_2, \end{aligned} \end{equation}
where
\begin{equation} \alpha=\frac{1}{\epsilon_1\ln(b/d) + \epsilon_2\ln(d/a)} \end{equation}
is a constant.
The following textbook problem is solved in the three-dimensional version of this numerical experiment. (Section 6.1.4. Spherical interface between dielectrics)
The system consists of two conductors and two layers of dissimilar dielectric materials. The conductors and the two layers of dielectric materials are shaped as spherical shells. All the shells are concentric. The dielectric materials fill all the space between the conductors forming one interface. The interface is shaped as a sphere concentric with the shells. The inner conductor is at a potential of \(\Phi=\Phi_0\) Volts. The outer conductor is grounded. The figure below illustrates the system and the problem domain.
We would like to calculate the electric potential, electric field, and displacement.
Due to the spherical symmetry of the problem we can reduce it to one dimension. Our intension, 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 conductors from the problem domain and model them by setting up Dirichlet boundary conditions on the surfaces that delineate them. Consequently, the problem domain consists of two dissimilar dielectric materials shaped as spherical shells with Dirichlet conditions applied to the inner and outer boundaries.
This problem can be solved by analytical methods. The solution reads
\begin{equation} \begin{aligned} &\Phi_1=-\alpha\epsilon_2\Bigg[\frac{1}{r}-\frac{1}{d}-\frac{\epsilon_1} {\epsilon_2}\Bigg(\frac{1}{b} - \frac{1}{d}\Bigg) \Bigg]\Phi_0, \\ &\Phi_2=-\alpha\epsilon_1\Bigg[\frac{1}{r}-\frac{1}{b}\Bigg]\Phi_0, \end{aligned} \end{equation}
\begin{equation} \begin{aligned} &\vec{\nabla} \Phi_1 = \alpha\epsilon_2 \frac{1}{r^2}\Phi_0\hat{r}, \\ &\vec{\nabla} \Phi_2 = \alpha\epsilon_1 \frac{1}{r^2}\Phi_0\hat{r}, \end{aligned} \end{equation}
\begin{equation} \begin{aligned} & \vec{E}_1 = - \vec{\nabla} \Phi_1, \\ & \vec{E}_2 = - \vec{\nabla} \Phi_2, \end{aligned} \end{equation}
\begin{equation} \begin{aligned} & \vec{D}_1 = - \epsilon_1 \vec{\nabla} \Phi_1, \\ & \vec{D}_2 = - \epsilon_2 \vec{\nabla} \Phi_2, \end{aligned} \end{equation}
where
\begin{equation} \alpha=\frac{1}{\epsilon_1\bigg[\dfrac{1}{b}-\dfrac{1}{d}\bigg]+ \epsilon_2\bigg[\dfrac{1}{d}-\dfrac{1}{a}\bigg]} \end{equation}
is a constant.
Both versions of this numerical experiment, the two- and three- dimensional, solve the following boundary value problem by the finite element method:
\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big( \epsilon \vec{\nabla} \Phi \big)= 0 & \text{in} & \Omega & \text{(i)},\\ \text{(e)} &\Phi = \eta_1 & \text{on} & \Gamma_{D1} & \text{(ii)},\\ \text{(e)} &\Phi = \eta_2 & \text{on} & \Gamma_{D2} & \text{(iii)},\\ \text{(e)} &\Phi_{+} = \Phi_{-} & \text{on} & \Gamma_{I1} & \text{(iv)},\\ \text{(n)}&\epsilon_{+}\hat{n}\cdot\vec{\nabla}\Phi_{+}-\epsilon_{-}\hat{n}\cdot\vec{\nabla}\Phi_{-}=0&\text{on}&\Gamma_{I1}&\text{(v)}. \end{array} \end{equation}
This boundary value problem is a special case of the general static scalar boundary value problem. The values of the parameters \(\eta_1\) and \(\eta_2\) are calculated as
\[ \eta_1 = \Phi_0 \]
and
\[ \eta_2 = 0, \]
where \(\Phi_0\) is a constant.
First, the StaticScalarSolver::Solver class template is used to obtain the numerical solution, \(\Phi\). After that, the numerical solution is fed to the class templates StaticScalarSolver::ProjectPHItoE and StaticScalarSolver::ProjectPHItoD to calculate the electric field and displacement as
\[ \vec{E} = - \vec{\nabla} \Phi \]
and
\[ \vec{D} = - \epsilon \vec{\nabla} \Phi. \]
This experiment is implemented in accordance with the base code structure. This time we use second-order mapping and attach a spherical manifold to all faces of the mesh. That is, we declare the spherical manifold in the SolverINT template,
and attach it to the mesh as the following:
Then we pass 2 as the second parameter to the constructor of SolverINT class template to enable second-order mapping,
This value of the mapping degree is then distributed to the projectors,
The build process generates two executable files: int-ring and int-shell. That is, one executable file for the two-dimensional experiment, int-ring, and one executable file for the three-dimensional experiment, int-shell. To rebuild them change into int/build/Release directory and execute the following:
Then all executable files must be executed again. This can be done by changing into the int/bin/Release directory and executing run-all script there,
This will generate various files in the int/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 int/bin/Release/Data directory also contains gpi files. They can be used to visualize the calculated potential along the \(x\) axis.
Note that executable files require a set of meshes to be present in the int/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into int/gmsh directory and executing the following:
This will generate a set of globally refined meshes in int/gmsh/data.
The SettingsINT 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.
Both, two- and three- dimensional versions of the experiment were executed under the following conditions: \(a = 0.3[m]\), \(b=1.0[m]\), \(d=0.65[m]\), \(\epsilon_1 = 32\epsilon_0\big[\frac{C^2}{Nm^2}\big]\), and \(\epsilon_2 = 4\epsilon_0\big[\frac{C^2}{Nm^2}\big]\).
The figure below illustrates the numerically calculated electric potential. As expected, the potential is continuous on the interface, i.e., at \(x = 0.65[m]\). Its derivative, i.e., the electric field, is discontinuous.
The corresponding convergence table is shown below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) | \(\|e\|_{H^1}\) | \(\alpha_{H^1}\) |
|---|---|---|---|---|---|---|---|
| 1 | 9 | 1024 | 1088 | 8.12e-04 | - | 5.20e-02 | - |
| 1 | 10 | 1296 | 1368 | 6.42e-04 | 2.00 | 4.62e-02 | 1.00 |
| 1 | 11 | 1600 | 1680 | 5.20e-04 | 2.00 | 4.16e-02 | 1.00 |
| 1 | 12 | 1936 | 2024 | 4.30e-04 | 2.00 | 3.78e-02 | 1.00 |
| 2 | 9 | 1024 | 4224 | 5.80e-06 | - | 8.58e-04 | - |
| 2 | 10 | 1296 | 5328 | 4.08e-06 | 3.00 | 6.79e-04 | 2.00 |
| 2 | 11 | 1600 | 6560 | 2.97e-06 | 3.00 | 5.50e-04 | 2.00 |
| 2 | 12 | 1936 | 7920 | 2.23e-06 | 3.00 | 4.55e-04 | 2.00 |
| 3 | 9 | 1024 | 9408 | 3.49e-07 | - | 3.34e-05 | - |
| 3 | 10 | 1296 | 11880 | 2.18e-07 | 4.00 | 2.34e-05 | 3.00 |
| 3 | 11 | 1600 | 14640 | 1.43e-07 | 3.99 | 1.71e-05 | 3.00 |
| 3 | 12 | 1936 | 17688 | 9.77e-08 | 4.01 | 1.29e-05 | 3.00 |
The following notations were used in the table header:
The figure below illustrates the calculated potential. The shape of the plot is in agreement with the equations listed above.
The corresponding convergence table is shown below.
| p | r | cells | dofs | \(\|e\|_{L^2}\) | \(\alpha_{L^2}\) | \(\|e\|_{H^1}\) | \(\alpha_{H^1}\) |
|---|---|---|---|---|---|---|---|
| 1 | 5 | 3072 | 3474 | 9.15e-03 | - | 2.41e-01 | - |
| 1 | 6 | 6000 | 6622 | 5.89e-03 | 1.97 | 1.93e-01 | 1.00 |
| 1 | 7 | 10368 | 11258 | 4.10e-03 | 1.98 | 1.61e-01 | 1.00 |
| 1 | 8 | 16464 | 17670 | 3.02e-03 | 1.99 | 1.38e-01 | 1.00 |
| 2 | 5 | 3072 | 26146 | 1.80e-04 | - | 1.30e-02 | - |
| 2 | 6 | 6000 | 50442 | 9.21e-05 | 3.01 | 8.37e-03 | 1.96 |
| 2 | 7 | 10368 | 86450 | 5.33e-05 | 3.01 | 5.84e-03 | 1.98 |
| 2 | 8 | 16464 | 136474 | 3.35e-05 | 3.00 | 4.30e-03 | 1.98 |
| 3 | 5 | 3072 | 86450 | 1.08e-05 | - | 8.49e-04 | - |
| 3 | 6 | 6000 | 167462 | 4.47e-06 | 3.94 | 4.43e-04 | 2.92 |
| 3 | 7 | 10368 | 287786 | 2.18e-06 | 3.96 | 2.59e-04 | 2.94 |
| 3 | 8 | 16464 | 455198 | 1.18e-06 | 3.97 | 1.64e-04 | 2.96 |
The following notations were used in the table header:
The figure below presents the simulated electric field, \(\vec{E}\), and its magnitude, \( \big| \vec{E} \big| \).
The interface conditions for the electric field suggest that the tangential component of the electric field must be continuous on the interface. The normal component of the electric field must be discontinuous. The magnitude of the jump of the normal component is proportional to the magnitude of the total surface charge density (bound charge density plus free charge density) on the interface. In the two problems above the free surface charge density, \(\kappa_f\), is zero, implying that the jump in the normal component of the electric field is proportional to the magnitude of the bound surface charge density. This field behavior on interfaces between dissimilar materials can be modeled by the Nedelec finite elements.
The table below presents the corresponding convergence data.
| p | r | cells | dofs | \(\|e\|_{L^2}/\epsilon_0\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 1 | 9 | 1024 | 2112 | 5.20e-02 | - |
| 1 | 10 | 1296 | 2664 | 4.62e-02 | 1.00 |
| 1 | 11 | 1600 | 3280 | 4.16e-02 | 1.00 |
| 1 | 12 | 1936 | 3960 | 3.78e-02 | 1.00 |
| 2 | 9 | 1024 | 8320 | 8.58e-04 | - |
| 2 | 10 | 1296 | 10512 | 6.79e-04 | 2.00 |
| 2 | 11 | 1600 | 12960 | 5.50e-04 | 2.00 |
| 2 | 12 | 1936 | 15664 | 4.55e-04 | 2.00 |
| 3 | 9 | 1024 | 18624 | 3.34e-05 | - |
| 3 | 10 | 1296 | 23544 | 2.34e-05 | 3.00 |
| 3 | 11 | 1600 | 29040 | 1.71e-05 | 3.00 |
| 3 | 12 | 1936 | 35112 | 1.29e-05 | 3.00 |
The following notations were used in the table header:
The figure below illustrates the calculated electric field. Here again, the radial component of the field, i.e., the component normal to the interface, is discontinuous.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}/\epsilon_0\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 1 | 5 | 3072 | 10000 | 2.41e-01 | - |
| 1 | 6 | 6000 | 19220 | 1.93e-01 | 1.00 |
| 1 | 7 | 10368 | 32856 | 1.61e-01 | 1.00 |
| 1 | 8 | 16464 | 51772 | 1.38e-01 | 1.00 |
| 2 | 5 | 3072 | 76832 | 1.30e-02 | - |
| 2 | 6 | 6000 | 148840 | 8.37e-03 | 1.96 |
| 2 | 7 | 10368 | 255792 | 5.84e-03 | 1.98 |
| 2 | 8 | 16464 | 404600 | 4.30e-03 | 1.98 |
| 3 | 5 | 3072 | 255792 | 8.49e-04 | - |
| 3 | 6 | 6000 | 496860 | 4.43e-04 | 2.92 |
| 3 | 7 | 10368 | 855432 | 2.59e-04 | 2.94 |
| 3 | 8 | 16464 | 1354836 | 1.64e-04 | 2.96 |
The following notations were used in the table header:
The figure below presents the simulated displacement, \(\vec{D}\), and its magnitude, \( \big| \vec{D} \big| \).
The interface conditions for the displacement suggest that in absence of surface free-charge density, \(\kappa_f\), on the interface the normal component of the displacement must be continuous. This field behavior on interfaces between dissimilar materials can be modeled by the Raviart-Thomas finite elements.
The table below presents the corresponding convergence data.
| p | r | cells | dofs | \(\|e\|_{L^2}/\epsilon_0\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 1 | 9 | 1024 | 2112 | 1.62e-01 | - |
| 1 | 10 | 1296 | 2664 | 1.37e-01 | 1.47 |
| 1 | 11 | 1600 | 3280 | 1.17e-01 | 1.48 |
| 1 | 12 | 1936 | 3960 | 1.02e-01 | 1.48 |
| 2 | 9 | 1024 | 8320 | 1.66e-02 | - |
| 2 | 10 | 1296 | 10512 | 1.31e-02 | 1.98 |
| 2 | 11 | 1600 | 12960 | 1.07e-02 | 1.99 |
| 2 | 12 | 1936 | 15664 | 8.82e-03 | 1.99 |
| 3 | 9 | 1024 | 18624 | 3.28e-04 | - |
| 3 | 10 | 1296 | 23544 | 2.27e-04 | 3.11 |
| 3 | 11 | 1600 | 29040 | 1.64e-04 | 3.10 |
| 3 | 12 | 1936 | 35112 | 1.22e-04 | 3.09 |
The following notations were used in the table header:
The figure below illustrates the calculated electric field. Here again, the radial component of the field, i.e., the component normal to the interface, is continuous.
The corresponding convergence table is presented below.
| p | r | cells | dofs | \(\|e\|_{L^2}/\epsilon_0\) | \(\alpha_{L^2}\) |
|---|---|---|---|---|---|
| 1 | 5 | 3072 | 9600 | 1.82e+00 | - |
| 1 | 6 | 6000 | 18600 | 1.36e+00 | 1.29 |
| 1 | 7 | 10368 | 31968 | 1.07e+00 | 1.29 |
| 1 | 8 | 16464 | 50568 | 8.81e-01 | 1.29 |
| 2 | 5 | 3072 | 75264 | 2.74e-01 | - |
| 2 | 6 | 6000 | 146400 | 1.78e-01 | 1.93 |
| 2 | 7 | 10368 | 252288 | 1.25e-01 | 1.95 |
| 2 | 8 | 16464 | 399840 | 9.23e-02 | 1.96 |
| 3 | 5 | 3072 | 252288 | 1.80e-02 | - |
| 3 | 6 | 6000 | 491400 | 9.11e-03 | 3.04 |
| 3 | 7 | 10368 | 847584 | 5.20e-03 | 3.08 |
| 3 | 8 | 16464 | 1344168 | 3.22e-03 | 3.10 |
The following notations were used in the table header:
An inspection of the convergence tables above suggests that the convergence order of the \(L^2\) error norm can be described as
\[ \alpha_{L^2} \approx p + 1 \]
in the case of the electric potential and as
\[ \alpha_{L^2} \approx p \]
in the case of the electric field and displacement with the notable exception of the displacement modeled by the finite elements of the lowermost order. In the case of the displacement and \(p=1\) the convergence order is somewhat higher than \(p\).
From the results above we can conclude that the class templates StaticScalarSolver::Solver, StaticScalarSolver::ProjectPHItoE, and StaticScalarSolver::ProjectPHItoD are able to solve problems with discontinuity of permittivity on interfaces between dissimilar dielectric materials.
1.9.1