The elastic constant of the crystal reflects the degree of crystal’s response to external forces, and is an important parameter that determines the elastic properties of the material. In this paper, the First Principle method based on density functional theory is used to study the elastic properties of mineral crystals.
3.3.1 Simulation Process and Geometric Optimization
The calculation uses Materials Studio software. Input the unit cell parameters to establish the crystal parameters of calcite and dolomite, build the model shown in Fig. 3, and then select the CASTEP module for geometric optimization. The parameters are as follows: select Geometric Optimization in the Task, select Fine for the calculation accuracy, the system energy converges to1.0×10−5eV/atom, the average atomic force converges to 0.03 eV/nm, the tolerance deviation converges to 0.001 Å, and the stress deviation converge to 0.05 GPa, select BFGS algorithm. In addition, the exchange correlation function, cutoff energy, and K point parameter selection are tested, and the results are shown in Fig. 5.
The lattice parameters and energy of calcite and dolomite will change along with the simulation parameters. Thermodynamics shows that when the total energy is the smallest, the crystal structure is the most stable. According to the calculation, the calcite and dolomite are selected for geometric optimization under the conditions of exchange correlation function GGA/PW91, cutoff energy of 380eV, K point 3×3×1, and the optimized calcite lattice parameters are: a=b=5.05229 Å, c=17.2129 Å, the dolomite lattice parameters are: a=b=4.87871 Å, c=16.2701 Å.
3.3.2 Simulation Results and Analysis of Crystal Elastic Properties
The elastic constants and other elastic properties of crystals can give us a deeper understanding of the micro-mechanical properties of mineral crystals. The molecular simulation technology to calculate the elastic parameters of the crystal is based on the "stress-strain" method, that is, the elastic constant of the crystal is obtained by the stress change caused by the different strains. For trigonal crystals, such as calcite and dolomite, the "stress-strain" relationship can be expressed in matrix form, as shown in Eq. (1).
\(\left(\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{23}\\ {\sigma }_{31}\\ {\sigma }_{12}\end{array}\right)=\left(\begin{array}{cccccc}{C}_{11}& {C}_{12}& {C}_{13}& {C}_{14}& & \\ {C}_{12}& {C}_{11}& {C}_{13}& -{C}_{14}& & \\ {C}_{13}& {C}_{13}& {C}_{33}& & & \\ {C}_{14}& -{C}_{14}& & {C}_{44}& & \\ & & & & {C}_{44}& {C}_{14}\\ & & & & {C}_{14}& \frac{{C}_{11}-{C}_{12}}{2}\end{array}\right)\left(\begin{array}{c}{e}_{11}\\ {e}_{22}\\ {e}_{33}\\ {e}_{23}\\ {e}_{31}\\ {e}_{12}\end{array}\right) (\) 1)
Materials Studio software is used to calculate the elastic constants of optimized calcite and dolomite crystals. To calculate all the elastic constants, the crystal structure needs to be optimized in six strain modes, and the maximum strain is set to 0.003. To improve the precision of elastic constants, the total energy convergence standard is increased to 2.0×10-6e V/atom. Other convergence standards and K point settings are consistent with the previous geometric optimization. Calculate the elastic modulus C11, C12, C13, C14, C33 and C44 of calcite and dolomite. According to the elastic constant of the crystal, the Voigt-Reuss-Hill theory (Hill 1952; Pham 2003) can be applied to further calculate the crystal elastic properties of calcite and dolomite. Among them, the calculation formulas for the bulk modulus, shear modulus and Young's (elastic) modulus of the trigonal crystal are shown in Eq. (2).
$${B}_{V}=\frac{1}{9}\left(2{C}_{11}+2{C}_{12}+4{C}_{13}+{C}_{33}\right)$$
$${G}_{V}=\frac{1}{30}\left(7{C}_{11}-5{C}_{12}-4{C}_{13}+2{C}_{33}+12{C}_{44}\right)$$
$${B}_{R}=\frac{\left({C}_{11}+2{C}_{12}\right){C}_{33}-2{C}_{13}^{2}}{{C}_{11}+{C}_{12}-4{C}_{13}+2{C}_{33}}$$
$${G}_{R}=\frac{15}{2}\left[\frac{2{C}_{11}+2{C}_{12}+4{C}_{13}+{C}_{33}}{{C}_{33}\left({C}_{11}+{C}_{12}\right)-2{C}_{13}^{2}}+\frac{3{C}_{11}-3{C}_{12}+6{C}_{44}}{{C}_{44}\left({C}_{11}-{C}_{12}\right)-2{C}_{14}^{2}}\right]$$
$$B=\frac{1}{2}\left({B}_{V}+{B}_{R}\right)$$
$$G=\frac{1}{2}\left({G}_{V}+{G}_{R}\right)$$
The elastic parameters of calcite and dolomite crystals are shown in Table 4 and Table 5.
Table 4
Elastic parameters of calcite crystal (unit: GPa)
elastic parameters
|
present work
|
calculated
|
experiment
|
Junhua (2009)
|
Ayoub (2011)
|
Hearmon(1979)
|
Chienchih (2001)
|
C11
|
147.20
|
146.82
|
152.3
|
144.0
|
149.4
|
C12
|
23.68
|
47.87
|
57.05
|
53.9
|
57.9
|
C13
|
53.29
|
46.05
|
54.83
|
51.1
|
53.5
|
C14
|
15.54
|
-16.81
|
17.14
|
-20.5
|
-20.2
|
C33
|
85.16
|
91.76
|
87.77
|
84.0
|
85.2
|
C44
|
38.3
|
32.52
|
36.00
|
33.5
|
34.1
|
B
|
70.20
|
|
|
|
|
G
|
39.74
|
|
|
|
|
E
|
100.29
|
|
|
|
|
Table 5
Elasticity parameters of dolomite crystal (unit: GPa)
elastic parameters
|
present work
|
calculated
|
experiment
|
Titiloye(1998)
|
Bakri(2011)
|
Humbert (1972)
|
Pofei (2006)
|
C11
|
200.97
|
201.6
|
196.6
|
205
|
204.1
|
C12
|
33.23
|
68.4
|
64.4
|
71
|
68.5
|
C13
|
53.38
|
58.2
|
54.71
|
57.4
|
45.8
|
C14
|
14.05
|
10.7
|
22.45
|
-19.5
|
20.6
|
C33
|
106.55
|
105.4
|
110.01
|
112.8
|
97.4
|
C44
|
43.62
|
41.7
|
41.57
|
39.8
|
39.1
|
B
|
84.98
|
|
|
|
|
G
|
53.88
|
|
|
|
|
E
|
133.44
|
|
|
|
|
According to the lattice dynamics theory (Born 1955), the elastic constant Cij of the crystal needs to meet the mechanical stability, and the conditions of the structural stability of the trigonal crystal is shown in Eq. (3).
$${C}_{11}-\left|{C}_{12}\right|>0$$
$$\begin{array}{c}\left({C}_{11}+{C}_{12}\right){C}_{33}-2{C}_{13}^{2}>0\end{array}$$
$$\left({C}_{11}-{C}_{12}\right){C}_{44}-2{C}_{14}^{2}>0$$
3
Substituting the calculated elastic constants of calcite and dolomite crystals into Eq. (3), the calculation results meet the conditions, indicating that the calculated elastic constants Cij of calcite and dolomite crystals meet the structural stability standard, that is, it is feasible to calculate the elastic constants of crystals with application of First Principle.
Normally, the elastic constants C11 and C33 indicate the resistance of the material to the strain along the axial-axis and the radial-axes (Hadi 2019). When C11>C33, the compressibility along the axial-axis is weaker than the radial-axes, otherwise the opposite. It can be seen from Table 4 and Table 5 that both calcite and dolomite have relatively weak compressibility along the axial-axis. The difference in value between C12 and C44 (C12-C44) is defined as the Cauchy pressure, which can be used to differentiate the brittleness and ductility of the material. If the Cauchy pressure value is positive, the material is ductile, otherwise the material is brittle (Du 2008). Table 4 and Table 5 show Cauchy pressures of calcite and dolomite crystals are 8.14GPa and 19.18GPa, indicating that both calcite and dolomite crystals are brittle. The calculation shows that the bulk modulus, shear modulus, and elastic modulus of dolomite are greater than that of calcite, and because the overall chemical bond strength of dolomite crystal is greater than that of calcite, we speculate that there is a correlation between the overall chemical bond strength and the elastic modulus of the crystal. The stronger the chemical bond strength of the crystal, the greater the elastic modulus of it.