A Finite Volume Method using a Quadtree Non-Uniform Structured Mesh for Modeling in Electrical Capacitance Tomography

An electric field solver based on a finite volume method using refined structural mesh is proposed to implement a quadtree structure and estimate the electric flux in the mesh cell. Numerical experiments were carried out using uniform and non-uniform meshes to assess quality of numerical modeling. The proposed method of verification of the quality of numerical calculations based on circular symmetry of the electrical capacitance tomography (ECT) probe allows to assess the effectiveness of mesh refinement and to reduce the number of mesh elements. Experiments showed that even a moderate level of mesh refinement is sufficient to significantly reduce the simulation error that occurs in modeling of cylindrical probes. The reduced number of mesh elements and applied implementation of the quadtree ensures high speed of forward problem calculations.


Introduction
Electrical capacitance tomography (ECT) has been described by a nonlinear relationship of capacitance measurements of a pair of electrodes surrounding an examined volume with the electric permittivity distribution inside this volume [1]. Using linear approximation, a forward problem has been formulated using the following linear equation where c is the capacitance, e is the elements describing spatial permittivity distribution in the nodes of the discrete mesh and S is the Jacobian matrix describing sensitivity of a given capacitance measurement on a small change of permittivity in a small volume of space, i.e. in an element of this mesh [2]. Forward problem solving enables possibility of measurement simulation and verification of the tomographic sensor design [3]. By solving the inverse problem, the spatial distribution of electric permittivity is reconstructed from capacitance measurements.
Electric field modeling in a tomographic capacitance sensor is necessary to determine the sensitivity matrix of a tomographic system [4]. The basic and simple approach to calculate numerically the spatial distribution of electric potential is to use the finite difference method (FDM) in a regular structured discretization mesh [5]. In ECT, the geometry of most sensors is cylindrical, so the implementation of the FDM using the Cartesian mesh introduces discretization errors into simulation.
The most advanced technique to treat a partial differential equation (PDE) for electric field is the finite element method (FEM) in which space is represented by an irregular unstructured mesh. The main advantage of the FEM over the FDM is the quality of approximation of complicated structures including round objects. The FEM is flexible but computationally expensive due to complex mesh requirement [6]. The FEM in reconstruction of ECT images already has been applied [7,8]. To decrease computational complexity of the FEM some modifications were proposed, like a hybridization of the FEM solution with a perturbative approach for small pixel discretization [9] or the application of the element-free Galerkin method [10].
Another less popular method to treat the PDEs numerically is the finite volume method (FVM). In this method, the space is divided into small volumes or cells that constitute the mesh of elements. The FVM is well suited for vector fields where the flux is conserved [11]. Volume integral of a divergence term in a PDE is converted using the divergence theorem to surface integral of a vector field. These terms are then treated as fluxes at the surfaces of each finite volume. The flux should be conserved in a discretization cell. If faces are aligned with the Cartesian axes, numerical schemes are simple to derive [12]. The FVM is attractive in engineering because of its robustness and cheap implementation [11]. The FVM has advantage in relation to the FDM because the latter has problems with approximation of the derivatives at the points where the coefficients are discontinuous. The FVM in Cartesian mesh for 2D ECT has been proposed [13]. The FVM in a cubic mesh was implemented for nonlinear reconstruction in 3D ECT [14]. The FVM in unstructured mesh of triangles was proposed for 2D ECT [15]. Numerical approximation for the FVM in polar coordinates in 3D ECT has already been done [16]. One of the methods of implementing the FVM is the use of a structural mesh described by a quadtree or octree structure. The use of a refined structural mesh for electric field modeling has been widely used in other areas of scientific computing [17].
In ECT, the computational efficiency of a forward problem solver (Jacobian matrix calculation) is of special importance. For strongly non-uniform permittivity distribution, image reconstruction based on linear approximation gives poor results. In a case of complex object imaging nonlinear problem should be solved [8]. Because an analytic solution of such a problem is unknown, the iterative method like Gauss-Newton ought to be applied in which the Jacobian matrix is calculated in each step of iterative procedure [18]. Thus, in ECT a field modeling method which has an acceptable computational cost can handle circular geometry of objects is required.
In this paper the implementation of the FVM using nonuniform structured mesh is presented, which is new in ECT modeling. The suited method for effective meshing using a quadtree (octree for 3D) in Cartesian coordinates and numerical approximation of electric flux in the mesh cell has been elaborated.

Methods
To implement the FVM, the MATLAB environment was used. The toolbox for numerical modeling and image reconstruction in electric capacitance tomography which was previously used [5] by our group was modified by adding non-uniform discretization.

Space Description
To describe the objects in the modelled space a library of geometrical primitives with the algebra of objects was developed. Complex objects are defined using other objects through basic operations like logical sum, difference and product (Fig. 1). Along with the geometric parameters, electrical parameters of the objects are defined, such as electric permittivity and potential. Using the description of the objects geometry, the space is divided into small elements forming fine mesh, henceforth referred to as reference mesh.

Meshing
To minimize the numerical approximation error, the discretization mesh should be finer in areas where electric field gradient is high: at object boundaries, between areas with different permittivity, at high curvature surfaces, sharp corners and at small features.

Preparation of the Input Map
The 3D or 2D map of object regions located in the area of interest is generated. The size of the matrix of this map corresponds to the size of the reference mesh with the assumed smallest element size. The unique object identifier is entered into the map points covered by this object (Fig. 2).
Additionally, this map is modified by adding the pattern (lattice) in the selected regions covered by the special objects added to the numerical model. In regions where this pattern is added, the bisection algorithm divides space into smaller segments. In such a way the space between the adjacent electrodes in the ECT sensor, where a large

The bisection algorithm
The aim of the bisection algorithm is to divide given square (cube in case of a 3D model) into smaller parts if the area (volume) is non-uniform or covered by more than one figure (solid). The mesh obtained with the help of the bisection algorithm is presented in Fig. 3.
The bisection algorithm starts from the element covering the entire area of interest. The minimum and maximum values of the input map are found. If the difference between the minimum and maximum value is greater than the assumed threshold, the square (cubic) element is divided into four (eight) and the new nodes are added to the quadtree (octree). The same procedure is repeated for the children nodes until the difference between the minimum and maximum values of each mesh element is smaller than assumed threshold or when each element reaches the smallest size assumed in the reference mesh.
A quadtree (octree) data structure is used to describe the elements of the mesh (Fig. 4). The root element represents the whole area (volume) of interest. A tree node has children if it describes the element divided into smaller elements. Tree leaves represent mesh elements and store information about their coordinates and size. The algorithm stores only the list of leaf nodes instead of all tree nodes.
To completely describe the discretization mesh, the list of neighbors is attached to each mesh element to form the adjacency list. Knowledge of the neighborhood of each mesh element is required to calculate both the distribution of potentials and electric field gradients.

FVM Implementation
In a discrete numerical model, the space is described in finite number of mesh points. In the FVM each node of the mesh is surrounded by a small volume. The Gauss flux theorem relates the electric charge contained in such a volume with the electric flux through the surface of this volume. This law expressed in a term of electric displacement field D is given by the equation where Q is a charge enclosed in the surface S, D is an electric displacement field through this surface and dA is an oriented, infinitesimally small element of the surface. If in the regular structured mesh, the mesh element is a cube, the electric displacement flux is integrated over the faces of this cube using the formula where e f is the electric permittivity at the face f , u f is the electric potential at this face, n f is the normal vector to the face, A f is the face surface and dim is the space dimension (2D or 3D). In irregular meshes one face of cubic element borders one large element or many small ones. The face is divided into surfaces in contact with the faces of adjacent cells (Fig. 5). The flux is summed over all surfaces in contact with neighboring cells. The contact area corresponds to the face of the smaller element and is determined by the formula where w i , w j f ð Þ are, respectively, the width of the given i-th element of the mesh and the f -th neighbor which is j-th element of the mesh, dim is space dimension. The value of electric potential gradient at the surface between the cells depends on the potential difference in these mesh cells and is inversely proportional to the distance between the centers of these elements. The component of the electric potential gradient which is perpendicular to the cube face is given by: where u i , u j f ð Þ are, respectively, the potential values in the i-th cell and the f -th neighbor which is j-th element of the mesh. The value of electric permittivity at the junction surface is approximated using the formula: where e i is the electric permittivity of the given element, w i is the width of the given element; e j f ð Þ is the electric permittivity in the neighbor element and w j f ð Þ is the width of the neighbor element.For i-th mesh element where the boundary conditions are not set and there is no free charge, the equation for electric displacement flux is given by After using the approximations (4), (5) and (6) to substitute the parameters in Eq. (7), it becomes For i-th mesh element where the Dirichlet boundary conditions are set, the equation for potential has the form The equations for all mesh cells form a system of linear equations with a sparse matrix A and the column vector of constant terms B. The structure of the system matrix for the example mesh (Fig. 3) is shown in Fig. 6.

Solver
The nonzero coefficients of the matrix of the linear system for electric potential in mesh cells are calculated using the above approximations. The vector of constant terms is determined using the value of boundary conditions which are set on the measuring electrodes and on the shielding elements (guard Fig. 5 The discrete finite volume stencil for a given cell (big square) for flux calculation through the surface in contact with neighboring cell electrodes and external shield). One electrode from the set of L electrodes becomes a driving electrode and the others are sensing electrodes. In one measurement cycle, each electrode becomes a driving electrode (e.g. 10 V). Voltage equal to zero is driven on the measuring electrodes. With subsequent measurements the geometry of boundary conditions do not change, thus the state matrix also does not change. Only the potential value in the boundary points changes. A new vector of constant terms B i is generated for each i-th measurement. A linear system (11) with augmented matrix of constant terms with the number of columns corresponding to the L number of driving electrodes is solved to speed up calculations.
The backslash operator in MATLAB is used to solve the equation for potential distribution (u ¼ AnB).

Sensitivity Matrix Calculation
To calculate the sensitivity of m-th capacitance measurement to a permittivity change in the n-th small element of space the reciprocity theorem is used which gives the formula where C m is the capacitance between electrodes A and B, U AB is the potential difference between electrodes, e n is the permittivity value in the mesh element, u n and w n are potential values in this element adequately when A is a driving electrode and B is a sensing electrode and vice versa, and X n is the volume of the element [4]. Electric potential gradients ru n and rw n in the mesh cells are calculated using the Green-Gauss method [19]. The partial derivatives (12) calculated for all mesh elements n ¼ 1; ::; N and all measurements m ¼ 1; ::M form the Jacobian (sensitivity) matrix with the elements S mn ¼ dC m =de n . Using the sensitivity matrix S and linearization given by Eq. (1), it is possible to determine the capacitance measurements for the assumed permittivity distribution e.

Assessment of Mesh Refinement Quality
The trade-off between the computational cost and the numerical error is a key problem when choosing the mesh density. Modeling the ECT sensor with too low mesh density provides the largest error in determining the distribution of the electric field in the area between adjacent electrodes. This is due to the large curvature of the electrode edges, sharp corners, and a small space with a high difference of electric potential. For this reason, the relative error in determining the capacitance is much greater for the adjacent electrodes than for the opposite electrodes, although the capacitance of the latter is definitely smaller. In the case of a cylindrical ECT sensor with a uniform permittivity distribution inside it, the simulated capacitance values for all adjacent electrode pairs should be the same due to circular symmetry. In a coarse square mesh, these capacitance values differ due to the different discretization of the electrodes depending on the angular position of the adjacent electrode pair relative to the center of the sensor. The finer the mesh, the smaller discrepancy in the capacitance value of adjacent electrode pairs. This property is used to determine the acceptable level of mesh refinement which should be increased until the difference in capacitance values for adjacent electrodes become below the given threshold. We assume here that the level of mesh refinement determined for the most problematic area of the numerical model will also be sufficient for discretization of objects being placed inside the tomographic probe.

Results
To verify our elaborated method, the example tomographic sensor with 8 electrodes in one ring was modelled in a uniform and non-uniform discretization mesh. The quality of calculations was evaluated for similar number of elements in both meshes, i.e. similar computational cost (Table 1).
Obtained calculation times presented in Table 2 were achieved using computer equipped with 64 GB of RAM and six-core Intel Core i5-8500 CPU with 3 GHz base speed. The non-uniform mesh obtained with the help of the bisection algorithm is presented in Fig. 7. The mesh is densest in the space surrounding the electrodes, where large electric field gradient is expected. Elements size in that area is 1 512 of whole model width. In the center of the probe and beyond its screen the mesh is rarer and its elements there are two or four times larger.
The distribution of electric potential was calculated for both meshes: uniform and non-uniform. The potential in the sensor when 10 V is driven on the selected electrode are presented in Fig. 8. The outer screen and other electrodes are connected to ground. Figure 9 shows the electric field distribution (modulus of field vector) expressed in volt per meter. The image is zoomed into the area of driving electrode. The non-uniform mesh is denser exactly in the place where the field gradient is larger (near the edge of electrode).
The maps of spatial distribution of sensitivity for the measurement of mutual capacitance of opposite electrodes (1 and 5) calculated for both meshes are presented in Fig. 10. The integral of this distribution is proportional to the capacitance value for this pair of electrodes and corresponds to one of the minimum values in the capacitance plot shown in Fig. 11. The maps of the spatial distribution of sensitivity for the adjacent electrodes (7 and 8) are shown in Fig. 12. The integral of this distribution corresponds to one of the maximum value in the capacitance plot shown in Fig. 11. In the case of adjacent electrodes, a  small area with a small number of elements but very high sensitivity determines capacitance value for adjacent electrodes. The values in this area are two orders higher (Fig. 12) than the values for the opposite electrodes (Fig. 10). The mesh density in this area determines the accuracy of numerical calculation of the capacitance value. The plot of simulated capacitances using both meshes is shown in Fig. 11. The quality of both meshes was verified by assessment of capacitance value for adjacent electrodes. The capacitance value for all these pairs should be equal due to circular symmetry of the sensor. The relative standard deviation (RSD) of the capacitance value for adjacent electrodes was analyzed. The simulations were repeated for different orientations (angles) and different positions of the sensor relative to the center of the mesh to obtain independence from the position of the sensor elements in relation to the mesh. Results of this simulations are presented in Fig. 13. For the quadtree mesh, the RSD value is twice lower than for the uniform mesh when the angle changes and even three times lower when the position changes. The mean relative standard deviation (MRSD) for all repeated simulations was compared for the uniform and non-uniform mesh ( Table 2).

Conclusions
Due to the need for nonlinear image reconstruction in ECT, in which a forward problem is recalculated many times, a fast solver for electric field distribution is desirable. In response to this expectation, the implementation of the FVM method in non-uniform structural mesh was elaborated.
The level of mesh refinement should be selected to ensure the expected quality of calculations at a minimized computational cost. The method for determining the level of mesh refinement based on the observation of differences in capacitance values for pairs of adjacent electrodes in a cylindrical probe was proposed in this work.
Experiments showed that using this innovative method the effectiveness of mesh refinement can be assessed and the number of mesh elements can be reduced significantly which affects the speed of computation. Numerical simulations were carried out using uniform and non-uniform refined meshes to assess quality of numerical modeling of ECT sensor. For a uniform mesh, where small elements of numerical model are not described accurately enough, the capacitance values simulated for pairs of adjacent electrodes are subject to a significant error. Experiments showed that a moderate level of mesh refinement (realized using small number of additional elements relative to the  The refined mesh allows for more accurate calculations with the same computational complexity. The obtained time of forward problem solving allows for a nonlinear reconstruction of dynamic ECT images in an acceptable time. Funding Excellence Initiative-Research University grants by the Ministry of Science and Higher Education (PL).

Conflicts of interest Not applicable
Code availability Yes.

Data Availability Yes.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless  indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.