Research on Mechanical Properties of High-Pressure Anhydrite Based on First Principles

: This article focuses on the elucidation of a three-dimensional model of the structure of anhydrite crystal (CaSO 4 ). The structure parameters of anhydrite crystal were obtained by means of first principles after structure optimization at 0~120 MPa. In comparison with previous experimental and theoretical calculation values, the results we obtained are strikingly similar to the previous data. The elastic constants and physical parameters of anhydrite crystal were also studied by the first-principles method. Based on this, we further studied the Young's modulus and Poisson's ratio of anhydrite crystal, the anisotropy factor, the speed of sound, the minimum thermal conductivity and the hardness of the material. It was shown that the bulk modulus and Poisson's ratio of anhydrite crystal rose slowly with increasing pressure. The anisotropy characteristics of the Young's modulus and shear modulus of anhydrite crystal were consistent under various pressure levels, while the difference in the anisotropy characteristics of the bulk modulus appeared. The acoustic velocities of anhydrite crystal tended to be stable with increasing pressure. The minimum thermal conductivity remained relatively unchanged with increasing pressure. However, the material hardness declined gradually with increasing pressure.


Introduction
The development and utilization of underground resources has always been an important aspect of human energy development. With the development of science and technology and the shortage of energy, the development depth of underground space has improved every year [1]. Therefore, the research and development of deep-underground technology has important strategic and military significance. The difference in the geological and environmental conditions between the deep underground and Earth's surface is very large. The mechanical properties of deep rock mass are a critical factor of deep-underground projects, and anhydrite rock is one of the major components of rock mass, meaning that the mechanical properties of anhydrite make a great contribution to the mechanical properties of rock mass [2]. Therefore, the study of the mechanical properties of anhydrite rock is of great significance.
Mohamed Chikhaoui experimentally tested the hydraulic permeability, isotherm, shear properties and elastic parameters of gypsum soil [3]. Marius manufactured a new type of high-performance composite using gypsum material and tested its mechanical properties, such as tensile strength [4]. Austa studied the chemical mechanism for low-salinity enhanced oil-recovery effects in carbonates by using a preserved core material containing dissolvable anhydrite [5]. Scholars have also undertaken a large amount of research work that investigates the physical and mechanical properties of anhydrite crystals. Sytle M. Antao studied the spatial structure of anhydrite crystals and obtained the lattice constants of anhydrite crystals [6]. Ghazal Azimi et al. studied the mechanism and kinetics of transformation, rehydration and crystallization between gypsum and anhydrite [7][8][9]. Aquilano put forward the theoretical equilibrium and growth shape of anhydrite crystals by means of the theory of Hartman and Perdok and compared them with the morphology of natural samples [10].
With the development of quantum theory and quantum mechanics, many scholars have used first principles to theoretically study the physical and chemical properties of crystals. E. Narsimha Rao used first principles to study the mechanical and thermodynamic properties of fluorocarbon materials and obtained physical property parameters such as the crystal vibration mode, free energy and specific heat capacity [11]. Etienne Balan used first principles to add sulfate to different CaCO3 samples to study the stability of its crystal structure [12]. With the development of high-speed computing, some software companies have developed software to calculate first principles, which has accelerated the progress in material research. At present, the popular programs on the market for first-principles calculation are VASP and CASTEP. Weck studied the physical properties of anhydrite crystal using VASP software, such as the heat capacity, entropy and free energy [13]. Benmakhlouf and Walker used CASTEP software to study the energy band and thermodynamic properties of MgSO4 and carbonate [14,15]. In the early research, most scholars mainly used the first-principles method to study the physical and chemical properties of crystals, such as free energy, entropy and enthalpy. Later, Zhao et al. studied the elastic constant of CaCO3 using first principles and obtained the shear modulus, bulk modulus, Poisson's ratio and other mechanical parameters of CaCO3 crystal [16]. DV Korabel'nikov also studied the physical properties of anhydrite rock using first principles, such as its structure parameter, elastic constants and vibrational energy states [17]. However, the specific research conditions-such as the pressure-were not specified in the research paper, nor is the specific elastic constant.
The static pressure of deep rock mass is about 120 MPa, and, in this paper, CASTEP software is mainly used to study mechanical properties such as the structural stability and elastic constants of anhydrite crystals at 0~120 MPa. Therefore, this work plays an important fundamental role in the study of the mechanical properties of deep rock mass.

Research Methods
A three-dimensional model of anhydrite (CaSO4) crystal was first established. Because of the orthogonality of anhydrite crystals, there are nine independent elastic constants, which can be obtained through the following formula [18,19]: where  is the reduced Planck constant, m is the electronic mass, eff V is the electronic potential, substituting Equation(2c) into Equation(2b). The same is done for the iterative calculation. The iterative calculation is finally stopped by the control of the self-consistent field function.

Specific Parameter Setting
All calculations were based on density functional theory, and the calculation processes were performed using the CASTEP software package. The electronic potential was treated with ultra-soft pseudopotential approximation. The exchange-correlation potential energy was treated with PW91 of the generalized gradient approximation. The cut-off energy of the plane-wave was set at 500 eV. A convergence criterion of 10 −6 eV/atom was used for the electronic self-consistency loop. The sampling unit of k points in the Brillouin zone was 2 × 2 × 2. The optimization of the structure was calculated at six pressure points between 0~120 MPa. The values of the convergence thresholds in each set are given in Table 1.  Figure 1 is a 3D drawing of unit cell of anhydrite (CaSO4) and shows the spatial distribution of anhydrite crystal atoms under a pressure of 0 MPa. Table 2 shows the fractional coordinates of each atom. There are 24 atoms in the unit cell.

Properties of Unit Cell Under 0 MPa Pressure
The lattice constants were measured by Antao using synchrotron high-resolution powder X-ray diffraction [20]. Hawthorne also reported the lattice parameters of anhydrite crystals through experiments [13]. In addition, Weck and Korabel'nikov analyzed the anhydrite crystal lattice constant theoretically [13,17]. Table 3 shows the comparison of the lattice parameters between the data used in this paper and previous data. It is concluded that the lattice parameters used in this paper are in good agreement with the former experimental and theoretical values. It also suggests that the structural model of anhydrite crystal that is established in this paper is reliable.  According to the mechanical stability condition of the orthogonal crystal system [20]: It is known that anhydrite rock is stable at 0 MPa.

Relationship of Crystal Structure Parameters with Pressure Variation
Lattice parameters and densities of anhydrite crystals measured under different pressures are presented in Table 4. It is clear from the table that the lattice parameters decline slowly with increasing pressure. It shows the compressible property of anhydrite rock. The parameters b/a and c/a tend to stabilize at about 1.126 and 1.13, respectively. It follows from the above information that the compressible property of anhydrite crystals is almost identical in all directions. The crystal density increases slowly with the increasing pressure. There is a minimum density of 2.855 g/cm 3 at 0 MPa. In order to more clearly describe the changing relationship between the lattice constant and pressure, the lattice parameters and volume change ratio of the lattice under pressure are represented by a graph (as shown in Figure 2). The a0, b0 and c0 are the lattice parameters of anhydrite crystals under pressure of 0 MPa. It can be seen from the graph that the ratio of the lattice constant in each direction decreases with the increase of pressure.  A polynomial was used to fit the lattice parameter curve, and the functional expressions of the lattice parameters in the region [0,120] MPa were obtained as follows: The above equation was used to draw the curve of the fitting function, as shown in Figure 3. It can be seen from the figure that all the parameters have good fitting degrees. Therefore, the relationship between pressure and lattice parameters can be approximately described by this expression in the region 0~120 MPa.  Table 5 shows the elastic constants of anhydrite crystals under different pressures. The elastic constants fluctuate and tend to increase slightly with the increase of pressure. As can be seen from  (3) under a pressure of 0~120 MPa. The bulk modulus B reflects the volume change of the crystal under pressure, the shear modulus G reflects the degree of shear deformation of the crystal, and the Young's modulus E reflects the deformation capacity of the crystal in the axial tensile load.

Pressure-Dependent Elastic Constants and Mechanical Properties
Once the elastic constant of the crystal is determined, the bulk modulus B and shear modulus G of the crystal can be obtained by Voigt-Reuss-Hill approximation. The bulk modulus B and shear modulus G of the crystal can be calculated by the following formula [21]:       2  12  22  11  33  11  23  13  12  23  22  13  23  12 13 where B is the bulk modulus, G is the shear modulus, BR and GR are, respectively, the bulk modulus and shear modulus of Reuss, and BV and GV are the bulk modulus and shear modulus of Voigt, respectively. In addition, Young's modulus and Poisson's ratio can also be obtained by the following formula: Figure 4 shows the variation of modulus under different pressures. The figure shows that the bulk modulus rises slowly with the increase of pressure. The shear modulus and Young's modulus of anhydrite crystals oscillate with the increase of pressure. In summary, the value of the three parameters is stable under the pressure of 0~120 MPa. In 1954, Pugh proposed the use of the ratio of the bulk modulus to shear modulus to evaluate the toughness/brittleness of materials [22]. If the B/G value is large, the toughness of the material is good (i.e., the value is greater than 1.75 B/G, meaning the material shows ductile behavior); otherwise, it exhibits brittle deformation behavior. Figure 5a shows the variation in the ratio of the bulk modulus to shear modulus with the increase of pressure. As can be seen from the diagram, the value is always greater than 1.75 B/G and tends to grow with the increase of pressure. This illustrates the toughness of anhydrite crystals under conditions of low pressure. In addition, Frantsevich uses Poisson's ratios to reflect the brittle behavior of materials [23]. If υ < 0.26, this indicates that the material is a brittle material; otherwise, it is a ductile material. Figure 5b is the Poisson's ratio of the anhydrite crystal change with pressure. As can be seen in the picture, the Poisson's ratio of anhydrite tends to grow with the increase of pressure, and the Poisson's ratios are always greater than 0.26 under a pressure in the range of 0~120 MPa. This shows that anhydrite crystals exhibit ductile behavior in this region. The Cauchy pressure of a crystal is expressed as the difference between two given elastic constants [24]. Figure 6 shows the variation of the three Cauchy pressures. The higher the Cauchy pressure, the greater the toughness of the material; otherwise, the material is brittle. If the Cauchy pressure is positive, the material is characterized by ductile behavior; otherwise, the material exhibits brittle behavior. As can be seen from the graph, the Cauchy pressure is negative in direction c, and in the other two directions, the Cauchy pressure is positive. Therefore, the toughness of anhydrite in direction b is inferior compared to in the other two directions. Hardness is a material's resistance to elasto-plastic deformation under external load, which is related to the elastic constant, plasticity, strain, brittleness, strength and other characteristics of the material. Theoretically, the hardness (H) of the material can be obtained by the following equation [25]: In order to visually show the relationship between hardness and pressure of anhydrite crystal, the hardness as a function of pressure was plotted, as shown in Figure 7. It is well known that all crystals are anisotropic. Therefore, the study of the anisotropic properties of crystals is of great significance to crystal physics and engineering science. The compression anisotropy factors, shear anisotropy factors and general anisotropy factors can be obtained by the following formula [26,27]: Figure 8 is the variation of anisotropic characteristic parameters for anhydrite crystals under different pressures. The anisotropy factor is zero for isotropic crystals. Regarding the difference between the anisotropy factor and zero reflection degree of the anisotropic elasticity of crystals, Figure 8a shows that the compression anisotropy factor of anhydrite tends to increase slightly with increasing pressure. Figure 8b,c shows that the shear anisotropy factor and general anisotropy factor increase slightly with the increase of pressure. Therefore, the compression and shear anisotropic characteristics of anhydrite crystal have the same tendency with increasing pressure. In addition, the directivity of the Young's modulus reflects the elastic anisotropic characteristics of the crystal. In materials science, elastic anisotropy is the directional dependence of physical properties of materials. Most materials exhibit elastic anisotropy. The bulk modulus, Young's modulus, shear modulus and Poisson's ratio are directional.
The bulk modulus of uniaxial compression along directions a and c is obtained by the following:  13  23  33  12  2  12  22  11   22  13  23  12  11  23  13  12   12 1 As can be seen from Figure 9, Ba and Bc have exhibit similar trends as a function of pressure, showing a gradual increasing trend with the increase of pressure. The fluctuation of data may be caused by the convergence accuracy of the calculation program.   (15) where B is the bulk modulus, G is the shear modulus, E is Young's modulus, and S is the flexible tensor.

Pressure-Dependent Wave Velocity and Thermal Conductivity
The wave velocities are related to the symmetry of the crystal structure and have directionality.
Then, the average elastic wave velocity is Figure 11 shows the relationship between the longitudinal wave velocity, transverse wave velocity and average wave velocity of polycrystalline anhydrite and pressure. It can be seen from the figure that the polycrystalline wave velocity fluctuates slightly under different pressures and tends to stabilize between 0 MPa and 120 MPa. To study the directivity of the crystal wave velocity, the slowness surface equation is taken as follows: where l is the direction cosine vector of the wave vector, and v is the velocity of the wave.
Then, the matrix element of the Christoffel tensor of an orthogonal crystal system is [29] 2  3  33   2  2  44   2  1  55  33   3  2  44  23  32  23   2  3  44   2  2  22   2  1  66  22   3  1  55  13  31  13   2  1  66  12  21  12   2  3  55   2  2  66   2  1  11 t v , respectively, represent the first and second transverse wave velocities, and  is the mass density of the crystal. Then, the minimum thermal conductivity can be obtained by the following formula [30]: where v n is the density of atoms per volume. Because the total thermal conductivity is already treated as the summation of one longitudinal and two transverse acoustic branches, the equation might be suitable to study the anisotropic thermal conductivities of the crystal. Figure 12 shows

Conclusions
The structure and anisotropic elastic properties of anhydrite crystals were studied by the first-principles method at a pressure in the range of 0~120 MPa. The optimized structure constants of anhydrite crystal remain strikingly similar to the previous experimental values and theoretical calculation values reported by other scholars, and the elastic tensor of anhydrite crystals was obtained under different pressures. The mechanical stability was ascertained by using elastic constants at each stage. The results show that the elastic constants meet the stability conditions effectively under all pressures. By studying the Young's modulus, anisotropy factor, sound velocity and other properties of anhydrite, it was found that the value of the bulk modulus and the anisotropy factor increase slowly with increasing pressure from Figures 4 and 8, respectively. The anisotropy characteristics of the Young's modulus and shear modulus of anhydrite were consistent under different pressures, and the anisotropic characteristics of the bulk modulus were different from Figure 10. The sound velocity of anhydrite crystal is essentially unchanged under a pressure of 0~120 MPa. The minimum thermal conductivity also remains relatively unchanged with the increase of pressure. The hardness of anhydrite crystal declines with increasing pressure, as shown in Figure 7. In this paper, the results of this study carry important implications regarding the understanding of dynamic characteristics of deep rock mass.
Author Contributions: X.Z. and S.Y. conceived and designed the calculations; X.Z. and L.L. performed the calculations; G.H. and W.Z. analyzed the data; Z.L. contributed analysis tools; X.Z., S.Y. and Y.X. wrote the paper. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflicts of interest.