1. Introduction
Polycrystalline materials are widely used as an important material in many fields such as construction, medical, and aerospace. Recent studies have demonstrated that the deformation, damage, and even fracture of polycrystalline materials at the macroscopic level depend on the microstructural characteristics of the materials, so it is important to understand the deformation mechanism of polycrystalline materials at the grain level to study the mechanical properties of these materials in the macroscopic field. The microscopic characteristics of polycrystalline materials include the single crystal behavior, the distribution of grain sizes and grain orientations, and the polycrystal morphology, etc. Therefore, the development of optimization techniques for polycrystalline material design requires a topological structure model with appropriate description of the statistical, topological, and physical characteristics of the microstructure.
Modeling techniques for polycrystalline material microstructures have developed rapidly in recent decades, either by material testing methods based on one-to-one experimental measurements, or by computer simulation methods based on statistical laws and algorithms.
With the development of materials testing technology, experimental measurement techniques have been used to measure real microstructural topological information of materials (e.g., grain size distribution, first neighborhood number, and grain morphology). Ultrasonic excitation-based nondestructive detection techniques that rely on the reflections of high frequency acoustic waves are not very effective for imaging grains in metals. Subsequently, the application of X-ray computed tomography technology solved this problem. Simonovski et al. established a microstructure model of stainless steel based on X-ray diffraction contrast tomography data [
1]. This method achieved the visualization of polycrystalline material microscopic modeling. However, complemented techniques (e.g., electron backscatter diffraction (EBSD)) must be used to identify the distribution of grain orientation.
Introduced more recently, the 3D real tissue modeling method based on serial sectioning techniques is one of the most prospective methods due to its high-fidelity structure of the materials. This technique builds 3D microstructure models in 3D space by stacking the two-dimensional microstructure information obtained following serial sectioning. This technique still needs an EBSD system to obtain the grain orientation [
2]. Generally speaking, all these modeling methods based on experiments require high research funding and must deal with the huge amounts of data generated during the experiments, which require complex post-processing techniques. Furthermore, the transformation process between 2D and 3D models generates non-unique topological relationships due to the overlap of grains.
Due to the aforementioned difficulties in experimental methods and the difficulty in generating large-scale grains, computer modeling and simulation techniques have gained researchers’ attention with their advantages in investigating the microstructure of polycrystalline materials. Furthermore, they bring great convenience to the establishment of finite element models in the subsequent simulation of deformation and damage of crystalline materials.
Initially, researchers investigated the micro-mechanics of polycrystalline materials with numerical simulations based on regular two-dimensional shapes, including Baczmanski [
3], who proposed thousands of square grains to simulate the residual stresses in steel in plastic deformation, followed by Ortiz [
4], who used ortho-hexagonal grains for residual stress simulations and concluded that at least two hundred grains were required to make the calculations more accurate. These regular shape grain models are simple to establish and easy to mesh. However, during the crystallization processes, because of the mutual resistance between adjoining grains, the grains spontaneously form polyhedral irregular grains. The above square or ortho-hexagonal grain models do not match with the actual grain shapes and cannot reflect the grain irregularity and the grain inhomogeneity deformation of inner materials. Among the numerical techniques developed to generate representative models of polycrystalline microstructures, Voronoi tessellation techniques are generally considered to offer an excellent compromise between representativeness and simplicity of formulation. Therefore, the modeling methods based on the Voronoi diagram have been widely used in the field of polycrystal modeling in recent years.
With the improvement of computing power, polycrystal modeling methods based on the Voronoi diagram were developed, and polycrystal models were developed from simple 2D polygons to 3D polyhedral models that are more consistent with the actual grain morphology. One of the most representative 3D polycrystal modeling methods was proposed by Quey et al., which is a modeling method for generating 3D random models of large-scale grain polycrystals [
5]. This method uses the open-source software Neper to generate the topological information of the grain model and reconstruct the polycrystal model by topological information in reverse topology in CAE software. Neper runs on any Unix-like system, which requires a large number of commands to configure the system and compile the software, and the process is troublesome and time-consuming for the researchers.
In recent years, Liang et al. generated the basic Voronoi cells information with MATLAB, stored the cells’ information as data files according to a certain order, and imported the data files into ABAQUS through the Python interface to establish polycrystalline microstructure models [
6]. This method implements the Quey.R method in Windows, reducing the researchers’ usage requirements. A large amount of grain information data is generated through the Voronoi diagram function in MATLAB, and then the polycrystalline microstructure can be modeled by reverse topological reconstruction using grain information (such as vertices, lines, surfaces, bodies) in the ABAQUS package. It remains an arduous task to establish polycrystalline models.
There are several variations of the Voronoi tessellation algorithm available in the open literature; Poisson–Voronoi, Hardcore–Voronoi, and Laguerre–Voronoi formulations are widely used to statistically represent the real polycrystalline microstructures. Compared with the first two variants that have only limited control over the shape and size of the cells, the Laguerre–Voronoi formulation imposes constraints on the initial state [
5], thus allowing the formulation of the tessellation to be modified more deeply. The Laguerre–Voronoi formulation can be used to model a wider range of grain structures (e.g., metals [
7], foams [
8], and granular matter [
9]).
In addition, it is universally acknowledged that grain boundaries (GBs) govern many properties of polycrystalline materials. With the decrease of grain size, the volume ratio of grain boundary structures increases with more significant effects on macroscopic material properties, especially in the nanomaterials [
10,
11]. However, the above-mentioned research works model grain boundaries as ideal interfaces without thickness that does not consider the 3D solid structure of grain boundaries.
In this paper, a new implementation of the Voronoi diagram in Laguerre geometry is presented for the generation of numerical models of polycrystalline microstructures. This method directly models the grains in CAD software by repeating cutting of the corresponding cell with several cutting planes, which are created by the Voronoi diagram. The 3D grain boundary structure can be constructed in the microstructure models. There is no need for cumbersome data processing and reverse topology reconstruction. Moreover, the resulting models can be directly imported into various computer aided engineering (CAE) software (e.g., ABAQUS, ANSYS, COMSOL) for the simulation and analysis of polycrystalline materials. This paper is organized as follows: In
Section 2, the procedure of modeling polycrystalline microstructures is described. In
Section 3, the finite element homogenization method based on the concept of representative volume element (RVE) is implemented to evaluate the effective elastic properties of the
ceramic material. In
Section 4, statistical analyses of grain face and grain size distribution are performed with the models, and the macroscopic elastic properties of polycrystalline ceramic materials are simulated to verify the capability of the presented method.
2. Polycrystal Modeling
Voronoi tessellations and its variants provide an analytical formulation to reproduce the non-regularity of polycrystalline morphologies. According to Boots’s hypothesis on crystals [
12],
All crystalline nuclei have the same weight, appear at the same time, and remain fixed in the same location during the growth process.
The growth is uniform and isotropic with a constant rate.
Grain growth in a direction stops when two grain boundaries contact each other.
In the aggregate there are no voids in between the grains or grain overlapping.
Among all the numerical techniques for generating representative models of polycrystalline microstructures, the Voronoi-based tessellation techniques are generally considered to offer an excellent compromise between representativeness and simplicity of formulation. Additionally, the straight edges and the planar faces of the grains are advantageous for further spatial discretization in finite element models.
2.1. Voronoi Tessellation
The Voronoi diagram is essentially a spatial partitioning structure. In three-dimensional space
, given a finite set of
points
(hereinafter called nuclei), the space is divided into a series of regions using the nucleus. Each region is the set of nearest points to the nucleus, and in this way a boundary is set for the region to achieve spatial partitioning. The set is defined as follows:
where
is any point in the
n-dimensional space
, and
denotes the distance from point
to
.
Voronoi tessellations have been applied to spatially discretize models in a variety of fields such as astronomy, materials science, and biology. It is obvious that Voronoi tessellation has obvious advantages in simulating the grain growth results of polycrystalline materials, such as metals or ceramics, and establishing the finite element model of polycrystalline materials. This method partitions the space to form a 3D structure by randomly creating a nucleus in the 3D space. The method is simple, stores less data, and has mature algorithms through recent development, which makes it easier to achieve the establishment of polycrystal microstructure models.
With the development of Voronoi tessellation in the field of polycrystalline microstructure modeling, the limitations of this algorithm are exposed: the grain sizes of polycrystalline models being normally distributed and not necessarily compatible with real polycrystalline materials, and the singularity and non-tunable nature of Voronoi tessellation in polycrystal modeling. In particular, they do not represent a large range of real grain sizes and the presence of large grains within the microstructure [
13]. Therefore, it is imperative to design a core algorithm that can satisfy the requirements of polycrystal modeling.
2.2. Laguerre–Voronoi Tessellation
The Laguerre–Voronoi diagram is studied on the basis of the Voronoi diagram. The Laguerre–Voronoi diagram divides the 3D space into N regions by defining Laguerre distances, which are composed of the nearest points to the nucleus Laguerre distance. Laguerre distance is defined as follows:
The Laguerre–Voronoi diagram is defined by the following equation:
where
and
are the radius and center, respectively, of any sphere
in space
;
denotes the Laguerre distance from point
to circle
.
Thus, the N-dimensional space is divided into n Laguerre–Voronoi regions and the corresponding boundaries, which constitute the Laguerre–Voronoi diagram. Compared with the Voronoi model, the grain size of the Laguerre–Voronoi model is controlled by the Laguerre distance, which changes the distribution of grain size so that the grain size is more consistent with the real distribution of grains.
From the mathematical description of the Laguerre–Voronoi tessellation, the overlap between the power circles can be avoided by imposing a non-overlapping condition on the spheres representing the weights of the nuclei. The Random Close Packing of Spheres (RCPS) model is usually adopted as the conditioning method [
14,
15]. This conditioning method, adopted in this paper, generates dense three-dimensional packing of spheres with an arbitrary radii distribution.
Moreover, polycrystalline microstructure models with 3D grain boundaries have also been established based on the Laguerre–Voronoi tessellation formulation. As shown in
Figure 1,
are the radii of the power circles corresponding to the adjacent nucleus,
is the distance between nuclei,
is the thickness of the grain boundary, and the point
is the location of the cutting plane between two nuclei. According to the Laguerre–Voronoi diagram, the distance ratio from the point
to the nucleus
is denoted as
, and the tangential distances from the point
to the two spheres are equal. Thus, the position of the cutting plane can be calculated by the derivation of the equation as follows:
deduce that
It is evident that when the grain boundary thickness is set to zero where the grain boundary is considered as ideal interfaces without thickness, the distance from the cutting plane to the nuclei is , and if the thickness of the grain boundary is set to , this value is equal to . In addition, when the RCPS model has the same radius of each power circle, the polycrystalline model is converted to Voronoi tessellation.
In this paper, a new polycrystalline microstructure modeling method was implemented. The flow chart is shown in
Figure 2. Firstly, spherical particles corresponding to the number of grains are randomly generated in the RVE region, and the positions and radius sizes of these particles are adjusted so that they meet the stacking requirements [
16,
17], as shown in
Figure 3. The seed points are established using the spherical center coordinates of the spherical particles. Secondly, a 3D Delaunay triangulation is created from the seed points. According to the data structure of vertices and edges in the Delaunay triangulation network, the adjacent seed points around each seed point which will achieve an effective cut on its cells can be efficiently determined. Thirdly, the corresponding cube is built by taking a seed point as the center and establishing the cutting planes between this seed point and its adjacent seed points. The grain model is obtained by cutting the corresponding cube with all the established cutting planes. Finally, Loop the above steps to obtain the grain model of each seed point, and a Boolean operation is conducted between each grain with the original cube to obtain the 3D grain boundary model; the program terminates when the final polycrystalline model is obtained (
Figure 4).
The procedure can be presented as follows:
The input and output of algorithm are as follows:
Input: The coordinates of the seed points , radius sizes of spheres and thickness of the grain boundary .
Output: Achieving the establishment of polycrystalline material models with solid grain boundary model.
- 1.
Create initial cube and seed points
- 2.
Delaunay triangulation of seed points and determine the effective cutting seed points around each seed point
- 3.
While the number of loops less than of seed points:
- 3.1.
Select a seed point to create a seed point cube
- 3.2.
Create cutting planes
If the cutting plane intersects with the grain model:
Cutting the seed cube with cutting plane
Else:
Finding next cutting plane
If there is no cutting plane
End of if
Go back to (3)
End While
- 4.
Boolean operation between each grain and the original cube to get the thickness of h grain boundary model
- 5.
Return Final Model
In particular, because of the presented formulation being implemented in the CAD package, it is possible to tessellate arbitrary desired solid models with Voronoi cells where the seed points are assigned randomly, as shown in
Figure 5.
4. Results and Discussion
In this section, statistical analysis and various finite element models of different polycrystalline models are used to illustrate the capability of the presented method for producing statistically representative microstructures.
Firstly, statistical analysis of grain face and grain size distribution are performed with the polycrystalline model. Compared with the Voronoi tessellation, the grain face and grain size distributions are obviously lognormal distribution.
Then, the polycrystalline model is assigned with linear elastic anisotropy properties, and the effective mechanical parameters are calculated by finite element homogenization method. The numerical results are compared with experimental measurement to validate the representativeness of the polycrystalline models.
4.1. Microstructure Analysis
4.1.1. The Distribution of Grain Size
The grain size distribution in real polycrystalline materials has been suggested to be lognormal in the previous studies [
22,
23,
24]. Its probability density function is described as
where
and
are standard deviation and arithmetic mean of ln
x, respectively.
For the purposes of comparison, two sets of polycrystalline microstructure models are generated separately in RVE based on the Voronoi and Laguerre–Voronoi diagram. To explore the effect of the number of grains, each group consists 125 and 500 grains. In order to have statistically valuable number results, ten different models—for a given number of grains (e.g., 125 and 500 grains)—were generated.
During the statistical analysis, the grain size distributions (
Figure 9) show that the grain size of the Voronoi tessellation-based polycrystalline model has a Gaussian normal distribution,
and the Laguerre–Voronoi tessellation-based one has a log-normal distribution:
where the
and
are standard deviation and arithmetic mean, respectively.
With the increasing of the number of grains, the grain size distributions of the polycrystalline models seem to become more stable, and the fitting degree is higher. Compared with the Voronoi tessellation-based polycrystalline models, the Laguerre–Voronoi tessellation-based ones allowing a wider range of grain structures to be modelled.
4.1.2. The Distribution of Face and Edge Number
Statistical analysis for the number of grain faces and edges was also performed for a polycrystalline model with 500 grains, and the statistical results are shown in
Figure 10. The number of grains faces in polycrystalline models is concentrated between 10 and 15, and the number of edges concentrates between 25 and 40, which is consistent with the actual grain microscopic morphology [
2,
5]. Furthermore, as with the distribution of grain size, the distribution functions of the number of grain faces and edges are fitted to the normal and log-normal distribution functions, respectively.
In addition, the 3D grain boundary can also be constructed in the polycrystalline model by the presented method in this paper. The distributions of grain size, and the number of grains faces and edges is shown in
Figure 11 for a grain boundary thickness
h of 0.1 and a number of grains of 200. The results show that the distributions of grain size, and the number of grain faces and edges of polycrystalline model grains that have solid grain boundaries and established by the Laguerre–Voronoi method conform to the log-normal distribution function shown in Equation (18). Moreover, the grain topology information in
Figure 11 is consistent with the experimentally measured grain information in a statistically regular manner [
25].
4.2. Predicting the Effective Elastic Modulus of Al2O3 Ceramic Materials
For predicting the elastic modulus of ceramic materials, linear elastic deformation is considered, there no material plasticity, and there are no defects. The effects of grain size, grain orientation, and grain boundary in the polycrystalline model to predict the effective elastic modulus of ceramic materials are investigated in this section.
4.2.1. Effect of Grain Size
An RVE of a polycrystalline model usually requires enough number of grains to convey consistent bulk properties. In order to investigate the adequate number of grains to construct an RVE of polycrystalline ceramic materials, different polycrystalline microstructure models with different numbers of grains are generated. Benedetti [
20] found that polycrystalline finite anisotropic materials (e.g., nickel, copper, gold) with cubic grain structures can be calculated with 20 grains for elastic properties within 10% error under specific boundary conditions. The
ceramic materials have finite anisotropy with linear elastic properties, so the minimum RVE structure is defined as consisting of 25 grains in this paper.
In this section, finite element models with 25, 50, 100, 200, and 500 grains are established to predict the effective elastic modulus of
ceramic materials (
Figure 12), and 10 different topologies are established for the given number of grains to reduce the influence of random errors on the calculation results.
According to references [
26,
27,
28], the elastic modulus interval of
material is obtained, as shown in
Table 1. The polycrystalline model in this paper is an ideal model and does not consider defects such as pores and microcracks. Therefore, the maximum value of 410 in the interval is selected as a reference, that is,
, to evaluate the accuracy.
Table 2 shows the effective elastic modulus of
ceramic materials predicted with these four sets of RVE models.
,
, and
are the average values of elastic modulus in
x,
y, and
z directions for 10 different topologies, respectively, and the relative error
is calculated using Equation (19). It can be concluded that the effective elastic moduli are all distributed within the same ranges of values measured experimentally and reported in the literature [
21]. With the increase of the number of grains, the errors of calculation results gradually decrease. Considering the balance between computational cost and accuracy, for a grain number larger than 100, the effect of single grain spatial distribution on the stress–strain response or overall elastic modulus can be ignored. The deformation of the RVE model of polycrystalline material with 100 grains under the six kinds of boundary conditions in this paper is shown in
Figure 13.
4.2.2. Effect of Grain Orientation
The grain orientation is determined by the order of atomic arrangement, and the atoms are randomly arranged during grain growth, so the polycrystalline materials are usually intergranular heterogeneous in microstructure, and the grain orientation can be simulated by the random factor method. However, the bulk properties of polycrystalline materials are isotropic. According to
Section 4.2.1, the adequate number of grains is taken as 100, and the effect of grain orientation on the effective elastic modulus is investigated in this section by assigning different grain orientations to the same finite element models.
In this section, the grain orientations are generated by the random factor method with randomly distributed orientations for each grain, and five different sets of grain orientations are assigned to the identical polycrystalline finite element models. Then, the effective elastic modulus of each model in three different directions is obtained, and the calculation results are shown in
Figure 14.
When the number of grains is 100, the effective moduli of elasticity for different grain orientations are in the range of 407–412 GPa and the average modulus of elasticity is within 1% error from the value obtained according to reference [
26,
27,
28] (
). It can be concluded that the grain orientation constructed by the random factor method can be used to simulate the randomness of grain orientation in the polycrystalline materials.
4.2.3. Effect of Grain Boundary
The difference in grain orientation leads to the formation of grain boundaries between adjacent grains during the grain growth process, and the existence of grain boundaries is the fundamental reason for the discontinuity of material properties, so it is often regarded as a weak point of material structure in the study of material structure. Nowadays, due to the small size of grain boundaries and the complex chemical composition, structure, and mechanical environment, there is no direct and accurate quantitative analysis method and experimental data for the thickness of grain boundaries and various material parameters. Therefore, the material properties of 3D grain boundaries are simplified as isotropic in this paper, and the elastic modulus of 3D grain boundaries are considered in three cases:
- 1.
The elastic modulus of the grain boundaries is 50% of the single crystal of the ceramic material.
- 2.
The elastic modulus of the grain boundaries is consistent with the elastic modulus of the ceramic material.
- 3.
The elastic modulus of the grain boundaries is 150% of the single crystal of the ceramic material.
For the three cases, the effective elastic modulus at different grain sizes was given as shown in
Figure 15. In the figure,
is the elastic modulus of grain boundaries, and
is the single crystal elastic modulus of
. When the value of
is 1.5 where the grain boundaries are considered to hinder grain deformation in the actual polycrystalline material, the overall macroscopic elastic modulus is strengthened as the average grain size decreases (the number of grains in the RVE model increases) where the density of grain boundaries increases, which is in accordance with the Hall–Petch equation. When the value of
is 0.5, with the decrease of grain size, the volume fraction of grain boundaries increases, and the macroscopic elastic modulus becomes smaller with grain refinement to achieve the toughening effect. It is concluded that when predicting the mechanical properties of polycrystalline materials, the presence of solid grain boundaries in the polycrystalline model results in the calculations being more consistent with the mechanical properties of the actual materials.
5. Conclusions
A new implementation of the Voronoi diagram in Laguerre geometry is presented in this paper for the generation of numerical models of polycrystalline microstructures, where the size and shape of the grains can be controlled, and the 3D grain boundaries can be modeled with a specified thickness. This method directly models the grains in CAD software. There is no need for cumbersome data processing and reverse topology reconstruction. Moreover, the resulting models can be directly imported into various finite element packages (e.g., ABAQUS, ANSYS, COMSOL, etc.) for the simulation and analysis of polycrystalline materials.
In order to verify the capability of the presented method, the statistical analysis of grain size, grain face, and grain edge distribution are performed with the established polycrystalline models, and the macroscopic elastic properties of polycrystalline ceramic materials are simulated. It is shown that the grain size, grain face, and grain edge distribution can be fitted to a lognormal distribution, compared with the normal distribution in Voronoi-based tessellation methods.
The effective elastic modulus of polycrystalline ceramic materials is predicted by the RVE-based finite element homogenization method. The effects of grain size, grain orientation, and grain boundary in the polycrystalline model materials are investigated. The predicted effective elastic moduli are compared with the experimental measurement, and the validity of the proposed method is verified.