Numerical modelling of rock materials with polygonal finite elements

This article presents some preliminary results on numerical modeling of rock materials with polygonal finite elements. A method to describe the rock microstructure based on Voronoi diagrams, representing the rock grain texture, is sketched. In this method, the minerals constituting the rock are represented as Voronoi cells which themselves are polygonal finite elements. A three-point bending problem under plane stress linear elasticity condition is solved in order to compare the performance of polygonal elements to ordinary finite elements. Moreover, it is demonstrated by solving the stress state in uni-axial compression that the heterogeneity described with the present method results in short-range tensile stresses which could initiate mode-I cracks.


Introduction
Polygonal finite elements have been drawing increasing attention during the last 15 years. In comparison to the standard finite elements, these elements offer greater flexibility in meshing arbitrary geometries, better accuracy in the numerical solution, better description of certain materials, and less locking-prone behavior under volume-preserving deformation [1]. With these elements, some numerical problems, such as checkerboarding in topology optimization, can be alleviated [2]. Polygonal discrete elements based on Voronoi diagrams is already an established method in geomechanics (see, for example, UDEC manual [3]). On the other hand, the disadvantages of the polygonal finite element method include less sparse system matrix and the need for a higher order numerical integration quadrature to achieve high-accuracy.
In the present paper, some preliminary results on a project of modelling rock materials with polygonal finite elements based on the Voronoi tessellation are presented. Namely, a method to describe the rock microstructure based on the Voronoi tessellation is sketched. At this preliminary stage of developments, only 2D linear elasticity problems are considered as numerical examples.

Polygonal finite elements
The standard isoparametric mapping from a reference element to the physical element can be used with the polygonal finite elements as well [1], as illustrated in Figure 1b. Here, the work by Talischi et al. [2] based on the Wachspress rational interpolation functions, providing C 0 -continuity, is followed. For a reference n-gon (see Figure 1b for hexagon), the Wachspress shape function at node i is where A(a, b, c) denotes the signed area of triangle a, b, c (see see Figure 1a). The numerical integration is carried out by sub-dividing the reference polygon into triangles and using a three integration points scheme for each triangle (resulting 3n integration points for each n-gon), as illustrated in Figure 1b. The variational (weak) formulation of the polygonal finite elements is fully standard, see [1] for details.

Rock microstructure description
Polycrystalline materials such as rocks can be naturally described with non-regular meshes of polygonal elements. Here, the minerals constituting the rock are represented by random clusters of Voronoi cells which themselves are polygonal finite elements. For this end, the PolyMesher Matlab code by Talischi et al. [4] is employed. This code generates 2D Voronoi diagrams (tessellations) consisting of centroidal Voronoi cells, which are then used as polygonal finite element meshes in the present work. The rock material is assumed to consist of tree main minerals (quartz, feldspar, biotite) which have their respective material properties. Thereby, the rock heterogeneity can be described. This method is illustrated in Figure 3a and b where two different realizations of the numerical rock sample with 2000 grains (Voronoi cells/polygon elements) are generated. In the first sample, the generating seeds are random yielding a realistic rock microstructure while in the second the polygons in the mesh are centroidal, i.e. the generating seeds are at the centroids of the polygons (CVT mesh), leading to a more regular cell distribution.

Numerical examples
The three-point bending problem of a linear elastic (plane-stress) beam is solved first in order to compare the performance of the polygonal elements to that of standard CST    Figure 2e shows that the 4-node rectangle and the Wachspress polygonal element perform better than CST element with coarse meshes (as expected).
Finally, the stress fields in uniaxial compression of the numerical rock specimen in Figure 3a are solved. The loading is such that the nominal stress for a homogeneous specimen is 130 MPa. The mineral elastic properties are as follows: Quartz (33 %) E = 80 GPa, ν = 0.17; Feldspar (50 %) E = 60 GPa, ν = 0.29; Biotite (17 %) E = 20 GPa, The results for the stress fields in both co-ordinate directions in Figure 3b and c display the high heterogeneity of the numerical specimen. Moreover, the transverse stress (x -direction) at many locations exceeds 20 MPa so that these short-range tensile stresses would generate mode-I cracks parallel to the loading direction.

Conclusions
Polygonal finite elements method was employed in this paper for numerical description of rock materials. The rock microstructure was modelled by Voronoi diagrams with the constituting minerals being represented as random clusters of Voronoi cells (polygonal finite elements). The numerical tests demonstrated that the performance of even highly non-structured meshes is acceptable for present purposes. The uniaxial compression test under a nominal stress of 130 MPa showed that the heterogeneity of the numerical specimen causes short-range tensile stresses that exceed the tensile strength of most rocks. Therefore, it can be concluded that the present approached has a good potential to be further developed to include a fracture model. This avenue will be pursued in the future studies of the approach.