Bifurcation analysis field potential in isotropic materials defined by semi-spherical layers using finite elements

In this document, the finite element method is developed in three dimensions to find the potential field in a region composed of hemispherical shells defining subregions with different materials, the numerical solution is made by FlexPDE software version 7.12 and the general process to solve is shown Problems of three-dimensional differential equations using the Galerkin criterion to approximate functions by the linear combination of functions of form, thus, the differential problem is reduced to solve a set of algebraic equations to find the coefficients of the linear combination.


Introduction
In the solution of differential equations by means of numerical analysis, the numerical finite element method allows a quick response in comparison with analytical solution methods dividing the domain of the function into subdomains with predefined geometries in each of which some known functions are defined called form functions or "Shape Functions" in such a way that border conditions are met, this reduces the problem to finding a solution in finite subdomains or elements, hence it gets its name [1,2].

Steps for the solution of a differential equation using finite elements
To solve a problem with finite elements, the most important part that depends on the engineer is the topology of the elements with which the domain is going to be subdivided as this affects the approximation of the solution. The domain will be divided into an equivalent system of elements with associated nodes, choosing the type of element appropriately to get as close as possible to the geometry of the domain. (Figure 1).

Select a function of form
Select a function of form in each element we define a function of defined form in the domain of each element and zero outside the domain, the most used functions are linear, quadratic and cubic polynomials because they are simpler to work with the formulation of finite elements, even so, trigonometric series can also be used [4]. The selected linear trigonal element is composed of 4 nodes whose function is (Equation (1)) [1].
For the four nodes located at the vertices of the tetrahedron, four equations are generated (Equation (3)).
The set of equations shown above can be represented as a matrix (Equation (4) and Equation (5)).
The scalar ϕ 9 is the value of the function at the point (x 9 , y 9 , z 9 ); by clearing α from the Equation (5) we have Equation (6).
The matrix [1 x y z]C :$ defines the set of functions in such a way Equation (8).
To facilitate the solution of the system, a system of volumetric coordinates can be obtained by subdividing the element into four tetrahedrons which result from the union of the four vertices of the tetrahedron with a point P contained within the element and dividing it by the total volume of the tetrahedron, this is Equation (10) [1]: Each V 9 is the volume of each of the tetrahedrons that result from dividing the element, such as , from there, the linear form functions are defined as Equation (11): Using the area coordinates, the volume integrals can be easily calculated using the following Equation (12) [1].

Calculate the stiffness matrix and mass matrix for the elements
The matrices of stiffness and mass are obtained from the set of equations resulting from the development of the weak formulation in each element once the boundary conditions have been replaced in these equations, obtaining these matrices represents a great human effort, but a computer can calculate them quickly and accurately. From the formulation obtained in literal b, we have for each element (Equation (13)) [2].
For the case of tetrahedral elements, there are four functions of form per element, and in total a set of 4 * N equations where N is the number of elements in which the domain is subdivided, such as the function f(x, y, z) is approximate as Equation (14).

Numerical experiments
The solution of the Laplace equation is proposed to find the potential field in a region defined by hemispherical shells in which each one has different materials (Equation (16)).
Being V the potential the region and the parameter ϵ depends on the material. As shown in Figure 2, each shell is composed of different materials represented by different colors and constants K. The boundary conditions are: On the outer surface (red surface in Figure 2) the potential is zero, and on the inner surface, (blue-green surface) the potential is one, the software chosen for the solution of this problem is FlexPDE version 7.12, in this software, the definition of the regions is done in parts, first the inner surface is defined (lower half of the red surface) and then a layer is defined, which must always be between surfaces, the layers are the space between two surfaces, and that is where each constant K is defined, so it is built from bottom to top surface and layer, until reaching the upper surface (upper half of the red surface), the construction process of the surfaces is shown in Figure 3 construction of the surfaces that define the domain of the problem [3,4]. The development of the Equation (17) is as follows.  In Figure 4, three-dimensional mesh showing the different surfaces that are each defined by a different color you can see in different color each surface declared in exPDE, this has nothing to do with the value of K for each region, you can also see the tetragonal elements that subdivide the domain whose surfaces are triangular [7,8].
In each individual element is defined a function that delivers the potential in each internal point, the union of the solutions in all the elements is the approximation of the solution in the whole domain when making a cross section in Z = 0.1, the equipotential lines in said surface can be observed ( Figure 5) [9,10].  It is appreciated as in the interior material (with a K = 3) there are more abrupt changes, and in the middle (with a K=50), the solution is more uniform, this could be exploited by increasing the size of the elements in the intermediate region and decreasing the size of the elements of the innermost region without altering the global quantity of finite elements to increase the resolution in the region where the solution varies more with the position, because it is not generally known a priori in which regions require more resolution due to large variabilities of the solution could be carried out in two stages [11,12].   Figure 6 shows the electric vector field in a portion of the surface at Z = 0.15, again the differences in the behavior of the solution according to the material are observed, which shows the utility of the method in circumstances with diverse geometries such like this case in which the domain is defined by semi cylindrical shells with different materials whose analytical solution represents a greater computational effort with a maximum error of 1.329% compared with the solution with finite elements [13].

Conclusion
As can be seen in Figure 5, the electric field depends inversely on the permittivity of the material, this can be exploited by decreasing the size of the elements in the regions with less permittivity and increasing them otherwise.