In this paper we describe the method of adaptive mesh refinement, based on the estimation of eigenvalues of discretization matrix. This estimation based on Gershgorin Circle Theorem. This method can be used for unstructured meshes in two-dimensional problems as well as and in three-dimensional. The implementation of the grid adaptation algorithm was made within OpenFOAM open source library of continuum mechanics. This library consists of a set of modules for computational needs, modules written in C++. We give two numerical examples, which show the effectiveness of the proposed method of mesh adaptation.
volume as well as for the whole computational domain and for any number of control volumes.
One of the advantages of the finite volume method over finite difference methods is that it does not require a structured mesh. Finite volume methods are especially powerful on coarse nonuniform grids and in calculations where the mesh moves to track interfaces or shocks.
In general case it is impossible to know a priori how to design an optimal mesh, i.e. mesh with minimal number of cells still satisfying the defined tolerance of the computational error. For transient problems, where the flux is unsteady and the points of interest can reposition during the simulation, uniform mesh becomes ineffective and would need to be very fine to satisfy the error tolerance throughout the whole simulation. To solve this problem a scheme where the mesh self-adapts its structure upon some criteria can be used.
The refinement based commonly on the physical quantity gradient field (see [3] and [4]). Our algorithm based on discretization matrix eigenvalues estimation, which mark cells to be refined or likewise mark cells to be unrefined. This approach has a more effective use of cells and thereby lower computational cost.
Suppose that after equations have been discretized by FVM we obtain the system of algebraic equations, which can be written in matrix form as Ax = b where A – square × discretization matrix, − × 1 column vector of unknown variable, − × 1 right-hand column vector. The problem consists in finding vector x whose elements are the values of physical quantity in the cell centers.
It is known that the accuracy of the solution of the problem is largely connected
with conditionality of the matrix A. Perform the following transformation with the
expression (
Note that error of the k-th and (k+1)-th iteration respectively can be written as
follows:
Subtract (
Expression (
(M) 1
Thus we see that the eigenvalues of the matrix M = I - A depend on the matrix A elements and play an important role in achieving the required accuracy. In particular, they show whether the problem under these conditions converge, converge to the correct solution and how fast.
On other hand, it is known that the elements of matrix A and the right-hand side b are functions of the mesh spacing. Consequently, there is a relation between conditionality and mesh discretization step. Identification of this relations using approximated equations is not easy task. Moreover, in closed source software, this details are hidden. Therefore we construct procedures of direct analysis of matrix M to identify the relations of conditionality and solution accuracy with mesh and adaptive mesh refinement.
It is known that the solution of eigenvalues problem is significant and has high computational cost. Moreover, if the case is ill-conditioned, then small eigenvalues can be calculated with low accuracy. Therefore, we have sufficient interest to use simple estimates of eigenvalues, estimates which calculated by the matrix elements insensitive to its conditioning.
Gershgorin circles method [4] is a well-known method for the localization of eigenvalues. According to Gershgorin’s Theorem every eigenvalues satisfies: Mii Mij ,
i j where i 1, 2, ..., n .
Let di Mij . Then the set
ji
is called the ith Gershgorin disc of the matrix M. This disc
is the interior plus the boundary of a circle. The circle has a radius di and is centered
at (the real part of Mii , the imaginary part of Mii ).
(
Calculating of estimates of the eigenvalues based on equation (6) is the computationally simple task. However, sufficiently strong upper and lower estimations are possible only for diagonally dominant matrices: Mii Mij .
i j If the matrix is not diagonally dominant, then lower and upper estimations of eigenvalues are undefined.
We offer to predict conditionality and mesh quality associated with conditionality through the right boundary of the Gershgorin circle: Fi mii ij mij . (7)
This boundary can be easily calculated. The maximum value, which calculated among the boundaries for all Gershgorin circles, defines an upper bound for the maximum eigenvalue of the matrix. Because of change of mesh all eigenvalues “shift” together with appropriate Gershgorin circles. There is reason that changes in the mesh, which lead to increase in (7), may lead to increase small eigenvalues.
Based on these assumption we have the mesh adaptation algorithm, which based on the analysis of scalar field i 1, 2, ..., n , which formed by calculating the values of (7) for all mesh nodes. Then mesh adaptation was made based on this scalar field. The normalization of field F is performed before each iteration:
Fi normalised
Fi max(Fi )
Thus the values of the field F are within the semi-interval (0;1] and it let us to set in OpenFOAM to refine mesh for cells with close to 1 and unrefined cells with Fi around 0. Values of field F decrease while refining mesh and increase while coursing mesh.
In order to illustrate feasibility and effectiveness of the algorithm, we use following two examples [5]. Both examples use modified laplacianFoam OpenFOAM solver called laplacianFoamF. New solver has ability to work with dynamic mesh, which allows to adapt mesh.
We compare AMR based on temperature scalar field T, i.e. based on temperature gradient minimization (AMR T) and AMR based on described above scalar field F (8) – (AMR F). For every example we tried to get about the same amount of cells.
Initial geometry of first example, as shown on Fig.1, thin square plate. The length along the x-axis and z-axis is 100 meters, along y-axis is 1 meter. (8) has temperature T2 273C . On all other surfaces the temperature gradient is set to 0, that walls do not conduct heat (adiabatic walls).
During testing, we found that AMR based on scalar field F detects too large cells, but does not take into account the boundary conditions (see Fig. 4). This occurs due to the fact that the boundary conditions are contained in right-hand column vector b of Eq.1, but almost no effect on discretization matrix M.
AMR based on temperature gradient field T vice versa takes into account the boundary conditions, but skips “bad” cells.
Therefore for comparison was made third “hybrid” variant, which include first five iterations of AMR F and after that it continues with AMR T. For comparison, residual plots of three cases are given (see Fig. 6).
As can be seen from the figure, AMR FT has the best residual, AMR F and AMR T slightly worse.
In the proposed method mesh refinement based on discretization matrix conditioning. As can be seen from above two examples – the method does not provide a significant performance increase compared to AMR, based on temperature gradient minimization, but use AMR T not always convenient or possible. Such situations are possible in complex problems with dynamic geometry, multiphase flows, etc. Our proposed method allows to choose more suitable AMR settings than in case of AMR Т.
In This work was supported by the Ministry of Education and Science of the Russian Federation (project №"2.2335.2014/K").