Introduction

In this numerical experiment we will continue testing

class templates. This time we will test them on axisymmetric problem domains. Axisymmetric problem domains arise in problems that exhibit rotation symmetry. See introduction to the Axisymmetric - method of manufactured solutions (mms-axi/) numerical experiment for more details. This time we will test if the three class templates mentioned above can solve problems with an interface between dissimilar materials on axisymmetric problem domains. Both problems considered in the Interface between dielectrics (int/) numerical experiment exhibit rotation symmetry. We will reuse them. That is, in this numerical experiment we will solve the two problems of the Interface between dielectrics (int/) numerical experiment on two-dimensional axisymmetric problem domains.

Implementation

The first problem

The first problem of the Interface between dielectrics (int/) numerical experiment exhibits cylindrical symmetry. We convert the problem into an axisymmetric problem domain as illustrated 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, $\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 these two segments of the boundary. We mark these boundaries with ID=0 in the int-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 two 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 will apply the Dirichlet boundary conditions on the horizontal segments of boundary: $\Phi_0[V]$ on $\Gamma_{D1}$ and $0[V]$ on $\Gamma_{D2}$ as dictated by the formulation of the problem. The program will read in the boundary IDs and will apply the boundary conditions $\Phi = \Phi_0[V]$ and $\Phi = 0[V]$ to the boundaries with ID=1 and ID=3, respectively.

The second problem

The second problem of the Interface between dielectrics (int/) 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 two vertical segments 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 these two segments. We mark these boundary segments with ID=0 in the int-axi/gmsh/sphere.geo The program will read in the zero boundary ID's, and will do nothing with respect to application of the boundary conditions to these two 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 will apply the Dirichlet boundary conditions on the curved segments of the boundary: $\Phi_0[V]$ on $\Gamma_{D1}$ and $0[V]$ on $\Gamma_{D2}$ as dictated by the formulation of the problem. The program will read in the boundary IDs and will apply the boundary conditions $\Phi = \Phi_0[V]$ and $\Phi = 0[V]$ to the boundaries with ID=1 and ID=3, respectively.

The boundary value problem

The two problems above are described by the following boundary value problem.

\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 electrostatic field and displacement as

\[ \vec{E} = - \vec{\nabla} \Phi \]

and

\[ \vec{D} = - \epsilon \vec{\nabla} \Phi. \]

The program

This experiment is implemented in accordance with the base code structure. The build process generates two executable files: int-axi-cylinder and int-axi-sphere. That is, one executable file for a problem: int-axi-cylinder for the first problem, and int-axi-sphere for the second problem. To rebuild them change into int-axi/build/Release directory and execute the following:

user@computer .../int-axi/build/Release$ ./clean

user@computer .../int-axi/build/Release$ ./build

Then all executable files must be executed again. This can be done by changing into the int-axi/bin/Release directory and executing run-all script there,

user@computer .../int-axi/build/Release$ cd ../../bin/Release

user@computer .../int-axi/bin/Release$ ./run-all

This will generate various files in the int-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 int-axi/gmsh/data directory. If they are missing, they can be generated anew. This can be done by changing into int-axi/gmsh directory and executing the following:

user@computer .../int-axi/gmsh$ ./clean

user@computer .../int-axi/gmsh$ ./build

This will generate a set of globally refined meshes in int-axi/gmsh/data.

The SettingsINTAXI 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.

Simulation results

Both 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]$.

Electrostatic scalar potential

The first problem

The figure below illustrates the numerically calculated electric potential.

The plot on the right compares the result of this experiment with the result of the two-dimensional version of the Interface between dielectrics (int/) numerical experiment. As expected, both simulations yield the same result.

The convergence table generated by this numerical experiment is shown below.

p	r	cells	dofs	$\\|e\\|_{L^2}$	$\alpha_{L^2}$	$\\|e\\|_{H^1}$	$\alpha_{H^1}$
1	9	128	153	2.28e-04	-	1.45e-02	-
1	10	162	190	1.81e-04	2.00	1.29e-02	1.00
1	11	200	231	1.46e-04	2.00	1.16e-02	1.00
1	12	242	276	1.21e-04	2.00	1.06e-02	1.00
2	9	128	561	1.79e-06	-	2.66e-04	-
2	10	162	703	1.26e-06	2.99	2.10e-04	1.99
2	11	200	861	9.20e-07	2.99	1.70e-04	1.99
2	12	242	1035	6.92e-07	3.00	1.41e-04	1.99
3	9	128	1225	2.90e-08	-	6.29e-06	-
3	10	162	1540	1.82e-08	3.96	4.43e-06	2.97
3	11	200	1891	1.19e-08	3.98	3.24e-06	2.98
3	12	242	2278	8.16e-09	3.98	2.44e-06	2.98

The following notations were used in the table header:

p - the degree of the interpolating Lagrange polynomials that constitute the shape functions.
r - the number of nodes on transfinite lines, see discussion in here.
cells - the total amount of active cells.
dofs - the amount of degrees of freedom in the active cells.
$\|e\|_{L^2}$ - the $L^2$ error norm.
$\|e\|_{H^1}$ - the $H^1$ error norm.
$\alpha_{L^2}$ - the order of convergence of the $L^2$ error norm.
$\alpha_{H^1}$ - the order of convergence of the $H^1$ error norm.

The second problem

The figure below illustrates the numerically calculated electric potential.

The plot on the right compares the result of this experiment with the result of the three-dimensional version of the Interface between dielectrics (int/) numerical experiment. As expected, both simulations yield the same result.

The convergence table generated by this numerical experiment is shown below.

p	r	cells	dofs	$\\|e\\|_{L^2}$	$\alpha_{L^2}$	$\\|e\\|_{H^1}$	$\alpha_{H^1}$
1	5	128	153	5.65e-03	-	1.48e-01	-
1	6	200	231	3.64e-03	1.97	1.19e-01	0.99
1	7	288	325	2.54e-03	1.98	9.90e-02	0.99
1	8	392	435	1.87e-03	1.98	8.49e-02	0.99
2	5	128	561	1.26e-04	-	9.07e-03	-
2	6	200	861	6.48e-05	2.99	5.88e-03	1.95
2	7	288	1225	3.75e-05	3.00	4.11e-03	1.96
2	8	392	1653	2.36e-05	3.00	3.03e-03	1.97
3	5	128	1225	6.56e-06	-	6.47e-04	-
3	6	200	1891	2.76e-06	3.88	3.40e-04	2.89
3	7	288	2701	1.35e-06	3.92	1.99e-04	2.92
3	8	392	3655	7.36e-07	3.93	1.27e-04	2.94

The following notations were used in the table header:

p - the degree of the interpolating Lagrange polynomials that constitute the shape functions.
r - the number of nodes on transfinite lines, see discussion in here.
cells - the total amount of active cells.
dofs - the amount of degrees of freedom in the active cells.
$\|e\|_{L^2}$ - the $L^2$ error norm.
$\|e\|_{H^1}$ - the $H^1$ error norm.
$\alpha_{L^2}$ - the order of convergence of the $L^2$ error norm.
$\alpha_{H^1}$ - the order of convergence of the $H^1$ error norm.

Electric field

The first problem

The figure below illustrates the numerically calculated electric field. As expected, the normal component of the electric field is discontinuous on the interface.

The corresponding convergence table is shown below.

p	r	cells	dofs	$\\|e\\|_{L^2}/\epsilon_0$	$\alpha_{L^2}$
1	9	128	280	1.45e-02	-
1	10	162	351	1.29e-02	1.00
1	11	200	430	1.16e-02	1.00
1	12	242	517	1.06e-02	1.00
2	9	128	1072	2.66e-04	-
2	10	162	1350	2.10e-04	1.99
2	11	200	1660	1.70e-04	1.99
2	12	242	2002	1.41e-04	1.99
3	9	128	2376	6.29e-06	-
3	10	162	2997	4.43e-06	2.97
3	11	200	3690	3.24e-06	2.98
3	12	242	4455	2.44e-06	2.98

The following notations were used in the table header:

p - the degree of the FE_Q finite elements that were used to model the electrostatic potential, $\Phi$. The degree of the FE_Nedelec finite elements that were used to model the electrostatic field, $\vec{E}$, equals $p-1$.
r - the number of nodes on transfinite lines, see discussion in here.
cells - the total amount of active cells.
dofs - the amount of degrees of freedom in the active cells.
$\|e\|_{L^2}$ - the $L^2$ error norm.
$\alpha_{L^2}$ - the order of convergence of the $L^2$ error norm.

The second problem

The figure below illustrates the numerically calculated electric field. As expected, the normal component of the electric field is discontinuous on the interface.

The corresponding convergence table is shown below.

p	r	cells	dofs	$\\|e\\|_{L^2}/\epsilon_0$	$\alpha_{L^2}$
1	5	128	280	1.48e-01	-
1	6	200	430	1.19e-01	0.99
1	7	288	612	9.90e-02	0.99
1	8	392	826	8.49e-02	0.99
2	5	128	1072	9.07e-03	-
2	6	200	1660	5.88e-03	1.95
2	7	288	2376	4.11e-03	1.96
2	8	392	3220	3.03e-03	1.97
3	5	128	2376	6.47e-04	-
3	6	200	3690	3.40e-04	2.89
3	7	288	5292	1.99e-04	2.92
3	8	392	7182	1.27e-04	2.94

The following notations were used in the table header:

p - the degree of the FE_Q finite elements that were used to model the electrostatic potential, $\Phi$. The degree of the FE_Nedelec finite elements that were used to model the electrostatic field, $\vec{E}$, equals $p-1$.
r - the number of nodes on transfinite lines, see discussion in here.
cells - the total amount of active cells.
dofs - the amount of degrees of freedom in the active cells.
$\|e\|_{L^2}$ - the $L^2$ error norm.
$\alpha_{L^2}$ - the order of convergence of the $L^2$ error norm.

Displacement

The first problem

The figure below illustrates the numerically calculated displacement. As expected, the normal component of the displacement is continuous on the interface.

The corresponding convergence table is shown below.

p	r	cells	dofs	$\\|e\\|_{L^2}/\epsilon_0$	$\alpha_{L^2}$
1	9	128	280	6.54e-02	-
1	10	162	351	5.53e-02	1.42
1	11	200	430	4.76e-02	1.43
1	12	242	517	4.15e-02	1.44
2	9	128	1072	5.66e-03	-
2	10	162	1350	4.49e-03	1.96
2	11	200	1660	3.65e-03	1.97
2	12	242	2002	3.03e-03	1.97
3	9	128	2376	1.26e-04	-
3	10	162	2997	8.71e-05	3.10
3	11	200	3690	6.28e-05	3.11
3	12	242	4455	4.67e-05	3.11

The following notations were used in the table header:

p - the degree of the FE_Q finite elements that were used to model the electrostatic potential, $\Phi$. The degree of the FE_RaviartThomas finite elements that were used to model the displacement, $\vec{D}$, equals $p-1$.
r - the number of nodes on transfinite lines, see discussion in here.
cells - the total amount of active cells.
dofs - the amount of degrees of freedom in the active cells.
$\|e\|_{L^2} / \epsilon_0$ - the normalized $L^2$ error norm. The error norm is normalized by the permittivity of free space, so its values are closer to unity.
$\alpha_{L^2}$ - the order of convergence of the normalized $L^2$ error norm.

The second problem

The figure below illustrates the numerically calculated displacement. As expected, the normal component of the displacement is continuous on the interface.

The corresponding convergence table is shown below.

p	r	cells	dofs	$\\|e\\|_{L^2}/\epsilon_0$	$\alpha_{L^2}$
1	5	128	280	1.28e+00	-
1	6	200	430	9.47e-01	1.33
1	7	288	612	7.38e-01	1.37
1	8	392	826	5.96e-01	1.39
2	5	128	1072	2.23e-01	-
2	6	200	1660	1.46e-01	1.89
2	7	288	2376	1.03e-01	1.92
2	8	392	3220	7.64e-02	1.94
3	5	128	2376	1.50e-02	-
3	6	200	3690	7.70e-03	3.00
3	7	288	5292	4.42e-03	3.05
3	8	392	7182	2.75e-03	3.08

The following notations were used in the table header:

p - the degree of the FE_Q finite elements that were used to model the electrostatic potential, $\Phi$. The degree of the FE_RaviartThomas finite elements that were used to model the displacement, $\vec{D}$, equals $p-1$.
r - the number of nodes on transfinite lines, see discussion in here.
cells - the total amount of active cells.
dofs - the amount of degrees of freedom in the active cells.
$\|e\|_{L^2} / \epsilon_0$ - the normalized $L^2$ error norm. The error norm is normalized by the permittivity of free space, so its values are closer to unity.
$\alpha_{L^2}$ - the order of convergence of the normalized $L^2$ error norm.

Table of Contents

Introduction

Implementation

The first problem

The second problem

The boundary value problem

The program

Simulation results

Electrostatic scalar potential

The first problem

The second problem

Electric field

The first problem

The second problem

Displacement

The first problem

The second problem