Elastic Constants and Homogenized Moduli of Monoclinic Structures Based on Density Functional Theory

Abstract

Elastic constants and homogenized properties of two monoclinic structures (gypsum and tobermorite) were investigated by first-principles method. The gypsum (chemical formula of CaSO4•2H2O) is an evaporite mineral and a kind of hydration product of anhydrite. Besides, the 11 Å tobermorite model (chemical formula: Ca4Si6O14(OH)4·2H2O) as an initial configuration of C-S-H structure is commonly used. Elastic constants are calculated based on density functional theory (DFT), which can also contribute to provide information for investigating the stability, stiffness, brittleness, ductility, and anisotropy of gypsum and tobermorite polycrystals. In addition, based on elastic constants (13 independent constants) of the monoclinic gypsum crystal, the elastic properties of polycrystals are obtained. The bulk modulus B, shear modulus G, Young’s modulus E, and Poisson’s ration ν are derived. Therefore, it is fairly meaningful to study the elastic constants to understand the physical, chemical, and mechanical properties of two monoclinic structures. Elastic constants can be used as the measure criterion of the resistance of a crystal to an externally applied stress. The calculated parameters are all in excellent agreement with reference.

1. Introduction

The density functional theory (DFT) is commonly used to study the crystal structure, lattice energy, the equation of state, the electronic bandgap, and vibration spectra properties [1]. Based on the kinetic energy density functional of Thomas [2] and the exchange-correlation effects of Dirac [3], DFT has been greatly developed by Kohn and Sham (KS) [4], who have established the fundamental approximation theorem on the functional status to describe real systems by electronic structure calculations. The eigenvalues of KS equations have no physical meaning, and the ionization energy is in the opposite state direction [5]. Moreover, one proposed approach is to introduce the eigenstates to calculate multi-body (many-body calculation) on the basis of Monte Carlo calculations [6] and perturbation theory [7]. The calculation of elastic constants is preceded by full geometry optimization and the stress tensor calculation of a number of distorted structures at the atomic scale. Polycrystalline structure constituted by a single crystal structure contains a variety of information (e.g., orientation) and the properties of a single crystal, such as anisotropy. Within the mechanics of typical crystals structures, the transition from the micro- to the meso-scale (homogenization) and vice versa (localization) can be estimated. Homogenization is an idealized description of a statistical distribution inside the actual heterogeneous material. Once the continuity model is admitted, the concept of homogeneity is deduced from it [8]. For quasi brittle materials, Zhu et al. [9] have formulated the anisotropic model in the framework of Eshelby-based homogenization methods. X-rays diffraction measurement is one of the stress assay test methods in physics field, of which the stress is actually determined by the strain [10]. Diffraction-based stress analysis depends critically on the use of the correct diffraction elastic constants [11]. X-ray method to test the material stress and to obtain elastic constants [12] is commonly based on the Reuss model [13]. Elasticity of single crystal and mechanical properties of polycrystalline material have been closely integrated. Various calculations methods are compared to determine homogenized moduli of the polycrystalline material composed of a single crystal, for example, the certain stress of Reuss model [13] and the certain strain of Voigt model [14].

DFT as a first-principles theory and a solid band theory in quantum mechanics has own a great success in linking physical properties and molecular structure, the calculation with exact accuracy but for low computational efficiency for macromolecular structure, which can be used to calculate elastic constants of anisotropic crystals, the monoclinic gypsum, and tobermorite crystals, for example. The chemical formula of gypsum is CaSO₄•2H₂O, which is an evaporite mineral and a kind of hydration product of anhydrite (chemical formula: CaSO₄). Moreover, the 11 Å tobermorite model (chemical formula: Ca₄Si₆O₁₄(OH)₄·2H₂O) as an initial configuration of C-S-H structure is commonly used. Since Young’s modulus parameters of gypsum and C-S-H are important to the multi-scale model [15], elastic constants of the gypsum crystal are investigated. The crystal is monoclinic, with 13 independent constants. For the homogenization of elastic deformation, especially for polycrystalline structures, the traditional Reuss-Voigt-Hill method is used to calculate the elastic moduli of monoclinic structures. Based on the ab initio plane-wave pseudopotential density functional theory method mentioned earlier, we focus on the monoclinic crystals to estimate their homogenized elastic moduli.

2. Theoretical calculation by density functional theory (DFT)

Despite the above advantages of DFT, however, the resolution of a system by Kohn-Sham equations involves difficulties due to an infinite number of electrons. These electrons maybe changed under an effective potential generated by an infinite number of cores or ions.

2.1. Equation of the theoretical approximate solution

From a microscopic point of view, Schrödinger equation describing a periodic crystal system composed of atomic nuclei n in mutual interaction and electron spin σ_i is positioned

and

respectively.

E1

Hamiltonian, in simple cases, consists of five terms: the kinetic energy of the electrons and nuclei, and the various interactions between them.

E2

The possible analytical representation and resolution of such a problem become a difficult task due to the limited memory of the computer tools. However, it is possible to reformulate the problem using appropriate theorems and approximations.

The fundamental principle approaches of mean field theory, in particular the DFT, are that any properties of an interacting particle system can be considered as a functional density in the ground state of the system n₀(r). Besides, the scalar function of the position n₀(r) essentially determines the wave functions of the system at the ground state and the excited states. Electronic and mechanical properties of a periodic crystal refer to solid state physics, quantum mechanics, and crystallography.

The crystalline ion movement of the electron is as

and assumes that the electron mobility (φ) does not depend on the speed nuclei but on their positions.

According to the Born-Oppenheimer or adiabatic approximation [16], the dynamics of the system (electrons and nuclei) is described. The electrons are assumed to react instantly to ionic motion. In electronic coordinates, the nucleus positions are considered as immobile external parameters.

E3

E4

where the last term of the Hamiltonian is constant and has been introduced in order to preserve the neutrality of the system and avoid the divergence of the eigenvalues. Clean the ground state of the system for fixed nuclear positions, total energy is given by the formula:

E5

This energy has a surface in the space coordinates that is said to be ionic Born-Oppenheimer surface. The ions move according to the effective potential energy, including Coulomb repulsion and the anchoring effect of the electron, which are as follows:

E6

E7

The dissociation degrees of freedom of electrons from those of nucleons, obtained through the adiabatic approximation, are very important, because if the electrons must be treated by quantum mechanics, degrees of freedom of ions in most cases are processed in a conventional manner.

This theorem/approach of Hohenberg and Kohn tries to make an exact DFT theory for many-body systems. This formulation applies to any system of mutually interacting particles in an external potential

, where the Hamiltonian is written as.

E8

DFT and its founding principle are summarized in two theorems, first introduced by Hohenberg and Kohn [17], which refer to the set of potential

and the density minimizing of Eq. (5).

The total energy of the ground state of a system for interacting electrons is functional (unknown) of the single electron density

$$E_{\mathit{HK}}\left[n\right]=T\left[n\right]+E_{\mathrm{int}}\left[n\right]+\int d^3{\mathit{rV}}_{\mathrm{ex}}\left(r\right)+E_{\mathit{nn}}\left(\overrightarrow{R}\right)\cdot F_{\mathit{HK}}\left[n\right]+\int d^3{\mathit{rV}}_{\mathrm{ex}}\left(r\right)+E_{\mathit{nn}}\left(\overrightarrow{R}\right)$$E9

As a result, the density n₀(r) minimizing the energy associated with the Hamiltonian (9) is obtained and used to evaluate the energy of the ground state of the system.

The principle established in the second theorem of Hohenberg and Kohn specifies that the density that minimizes the energy is the energy of the ground state

$$E^{\mathrm{BO}}\left[\overrightarrow{R}\right]=minE\left(\overrightarrow{R},n\left(\overrightarrow{r}\right)\right)$$E10

Because the ground state is concerned, it is possible to replace the wave system function by the electron charge density, which therefore becomes the fundamental quantity of the problem. In principle, the problem boils down to minimize the total energy of the system in accordance with the variations in the density governed by the constraint on the number of particles

. In this stage, the DFT can reformulate the problem rather than solve an uncertain functional F_HK(n).

2.2. The approximation approach of Kohn-Sham

The approach of Kohn-Sham system substitutes the interacting particles, which obeys the Hamiltonian in Eq. (3), by a less complex system easily solved. This approach assumes that the density in the ground state of the system is equal to that in some systems composed of non-interacting particles. This involves independent particle equations for the non-interacting system, gathering all the terms complicated and difficult to assess, in a functional exchange-correlation $E_{\mathit{XC}}\left(n\right)$.

E11

T is the kinetic energy of a system of particles (electrons) independently (non-interacting) embedded in an effective potential which is no other than the real system,

E12

The Hartree energy or energy of interaction is associated with the Coulomb interaction of the self-defined electron density.

E13

$$n\left(r\right)=\sum _{i=1}^{N_e}{\left|{\phi }_i,\left(r\right)\right|}^2$$E14

Solving the auxiliary Kohn and Sham system for the ground state can be seen as a minimization problem while respecting the density n(r). Apart from orbital function T_S, all other terms depend on the density. Therefore, it is possible to vary the functions of the wave and to derive the variational equation:

$$\frac{\delta E_{\mathit{KS}}}{{\mathit{\delta \phi }}_i^{\ast }\left(r\right)}=\frac{\delta T_S}{{\mathit{\delta \phi }}_i^{\ast }\left(r\right)}+\left[\frac{\delta E_{\mathit{ex}}}{\mathit{\delta n}\left(r\right)}+\frac{\delta E_{\mathit{\text{Hartree}}}}{\mathit{\delta n}\left(r\right)}+\frac{\delta E_{\mathit{xc}}}{\mathit{\delta n}\left(r\right)}\right]\frac{\mathit{\delta n}\left(r\right)}{{\mathit{\delta \phi }}_i^{\ast }\left(r\right)}=0$$E15

With the constraint of orthonormalization $⟨{\phi }_i\left|{\phi }_j\right.⟩={\delta }_{i,j}$, this implies the form of Kohn-Sham for Schrödinger equations:

$$\left({{\displaystyle \stackrel{\wedge }{H}}}_{\mathit{KS}}-{\epsilon }_i\right){\phi }_i\left(r\right)=0$$E16

${\mathit{\epsilon }}_i$ represents the eigenvalues, and ${{\displaystyle \stackrel{\wedge }{H}}}_{\mathit{KS}}$ is the effective Hamiltonian H,

$${{\displaystyle \stackrel{\wedge }{H}}}_{\mathit{KS}}\left(r\right)=-\frac{1}{2}{\nabla }^2+V_{\mathit{KS}}\left(r\right)$$E17

$$V_{\mathit{KS}}\left(r\right)=V_{\mathit{ex}}\left(r\right)+\frac{\delta E_{\mathit{\text{Hartree}}}}{\mathit{\delta n}\left(r\right)}+\frac{\delta E_{\mathit{xc}}}{\mathit{\delta n}\left(r\right)}$$E18

Eqs. (16)–(18) are known equations of Kohn-Sham, the density n(r) and the resulting total energy E_KS. These equations are independent of any approximation on the functional E_XC(n), resolution provides the exact values of the density and the energy of the ground state of the interacting system, provided that E_XC(n) is exactly known. The latter can be described in terms of Hohenberg Kohn function in Eq. (8)

$$E_{\mathit{xc}}\left[n\right]=F_{\mathit{HK}}\left[n\right]-\left(T_S\left[n\right]+E_{\mathit{\text{Hartree}}}\left[n\right]\right)$$E19

or more precisely,

$$E_{\mathit{xc}}\left[n\right]=⟨\stackrel{\wedge }{T}⟩-T_S\left[n\right]+⟨{\displaystyle \stackrel{\wedge }{V_{\mathrm{int}}}}⟩-E_{\mathit{\text{Hartree}}}\left[n\right]$$E20

This energy is related to potential exchange-correlation $V_{\mathit{xc}}=\frac{\partial E_{\mathit{xc}}}{\partial n\left(r\right)}$.

For the exchange-correlation functional, the only ambiguity in the approach of Kohn and Sham (KS) is the exchange-correlation term. It is subject to functional approximations of local or near local order of density that said energy E_KS can be written as

$$E_{\mathit{xc}}\left[n\right]=\int n\left(r\right){\epsilon }_{\mathit{xc}}\left(\left[n\right],r\right)d^{3}r$$E21

where ε_xc([n], r) is the exchange-correlation energy per electron at point r, it depends on n(r) in the vicinity of r. These approximations have made enormous progress in the field.

1. The approximation of the local density (LDA)

   The use of the local density approximation (LDA) in which the exchange-correlation energy $E_{\mathit{xc}}^{\mathit{LDA}}\left[n\right]$ is another integral over all space, assuming that ε_xc([n], r) is the exchange-correlation energy per particle of a homogeneous electron gas of density n

   $$E_{\mathit{xc}}^{\mathit{LDA}}\left[n\right]=\int n\left(r\right){\epsilon }_{\mathit{xc}}^{\mathit{bom}}\left[n\left(r\right)\right]d^{3}r=\int n\left(r\right)\left\{{\epsilon }_x^{\mathit{bom}}\left[n\left(r\right)\right]+{\epsilon }_c^{\mathit{bom}}\left[n\left(r\right)\right]\right\}{d}^{3}r$$E22

   The exchange term $E_x^{\mathit{bom}}\left[n\left(r\right)\right]$ can be expressed analytically, while the correlation term is computed accurately using the Monte Carlo by Ceperley Alder [18] and then set in different shapes [19].

   This approximation has been particularly checked to deal with non-homogeneous systems.

2. The generalized gradient approximation (GGA)

   The generalized gradient approximation (GGA) involves the local density approximation providing a substantial improvement and better adaptation to the systems. This approximation is equal to the exchange-correlation term only as a function of the density. A first approach (GEA) was introduced by Kohn and Sham then used by the authors of Herman et al. [20].

   This notion of GGA is the choice of functions, which allows us a better adaptation to wide variations so as to maintain the desired properties. The energy is written in its general form [21]:

   $$E_{\mathit{xc}}^{\mathit{GGA}}\left[n\right]=\int n\left(r\right){\epsilon }_{\mathit{xc}}\left[n,\left|{\nabla }_n\right|,\dots \right]d^{3}r=\int n\left(r\right){\epsilon }_x^{\mathit{bom}}\left(n\right)F_{\mathit{xc}}\left[n,\left|{\nabla }_n\right|,\dots \right]d^{3}r$$E23

   where ${\epsilon }_x^{\mathit{bom}}$ is the exchange energy of an unpolarized density n(r) system. There are many forms of F_XC, the most used are those introduced by Becke [22] and Perdew [23, 24].

2.3. Parameters of Bloch theorem and Brillouin zone

Different states of the Schrödinger equation for an independent particle in a system. By Kohn and Sham equations, responding to eigenvalue equation is as:

$${{\displaystyle \stackrel{\wedge }{H}}}_{\mathit{eff}}\left(r\right){\psi }_i\left(r\right)=\left[-\frac{{\mathrm{\hslash }}^2}{2m}{\nabla }^2+V_{\mathit{eff}}\left(r\right)\right]{\psi }_i\left(r\right)={\epsilon }_i{\psi }_i\left(r\right)$$E24

where the electrons are immersed in an effective potential

.

The effective potential has the periodicity of the crystal and may be expressed using Fourier series, in a periodic system:

$$V_{\mathit{eff}}\left(r\right)=\sum _m{V}_{\mathit{eff}}\left(G_m\right)exp\left({\mathit{iG}}_{m}r\right)$$E25

G_m is the reciprocal lattice vector:

$$V_{\mathit{eff}}\left(G\right)=\frac{1}{{\Omega }_{\mathit{\text{sell}}}}{\int }_{{\Omega }_{\mathit{\text{sell}}}}{V}_{\mathit{eff}}\left(r\right)exp\left(-{\mathit{iG}}_{m}r\right)\mathit{dr}$$E26

where ${\Omega }_{\mathit{\text{sell}}}$ is the volume of the original mesh.

As the translational symmetry, it is that states are orthogonal and conditioned by the limits of the crystal (infinite volume). In this case, the Eigen functions of KS are governed by the Bloch theorem: they have two quantum numbers: the wave vector k in the Brillouin zone (BZ) and the band index i, and this can be expressed by a product of a plane wave exp.(ik, r) and a periodic function:

$$\begin{array}{l}{\psi }_{i,k}\left(r\right)=exp\left(\mathit{ik},r\right)u_{i,k}\left(r\right)\\ u_{i,k}\left(r+R\right)=u_{i,k}\left(r\right)\\ R=\sum n_i{a}_i,n_i=1\dots N_i\end{array}$$E27

where R is the vector of direct space defined by a_i with $i\in \left\{1,2,3\right\}$ and N_i is the number of primitive cells in each direction ($N_i\to \infty$ in the case of perfect crystal).

Solving Eq. (24) is equivalent to increase the periodic function $u_{i,k}\left(r\right)$, in a database-dependent functions points $k:\left\{{\varphi }_j^k\left(r\right)\left|j=1\dots N_{\mathit{bas}}\left(k\right)\right.\right\}$:

$$u_{i,k}\left(r\right)={\sum }_j{C}_i^j{\varphi }_j^k\left(r\right)$$E28

where ${\varphi }_j^k$ is the wave function developed in a space of infinite dimensions; this means that j should be in principle infinite. But, in practice, we work with a limited set of basic functions, which imply that the description of ${\varphi }_j^k$ will approximate. That the selected database simply solves the system:

$$\begin{array}{l}\sum _{m\text{'}}{H}_{m,m\text{'}}\left(k\right)C_{i,m\text{'}}\left(k\right)={\epsilon }_i\left(k\right)C_{i,m}\left(k\right)\\ H_{m,m\text{'}}\left(k\right)=⟨{\phi }_{m,k}^j\left|{{\displaystyle \stackrel{\wedge }{H}}}_{\mathit{eff}}\right|{\phi }_{m,k}^j⟩\end{array}$$E29

where each point is a set of k eigenstates, the label having i = 1, 2, … obtained by diagonalization of the Hamiltonian in Eq. (29).

It is necessary to integrate the points k in the Brillouin zone. For a function f_i(k) where i defines the band index, the average value is

$$\overline{f_i}=\frac{1}{N_k}\sum _k{f}_i\left(k\right)\to \frac{{\Omega }_{\mathit{\text{cell}}}}{{\left(2\pi \right)}^d}{\int }_{\mathit{BZ}}{f}_i\left(k\right)\mathit{dk}$$E30

${\Omega }_{\mathit{\text{cell}}}$ is the cell volume of the original mesh in the real space and ${\left(2\pi \right)}^d/{\Omega }_{\mathit{\text{cell}}}$ of the cell volume of the Brillouin zone are determined using a sampling points k. Several election procedures exist for these points. Particularly those of Baldereschi [25], Chadi and Kohen [26], and Monkhorst and Pack [27] are the most frequently used.

3. Elastic constants and homogenized moduli of monoclinic structure

According to the crystal theory [28], any crystal lattice system contains six independent variables, namely the cell side length a, b, and c; unit cell angle α, β, and γ. Generally, the crystal under a certain deformation, temperature, and pressure can be described by the corresponding six-dimensional deformation tensor. The temperature and pressure will cause cell-deformed configuration tensor as

$$X^{P,T}=\left[\begin{array}{lll}{a}_x^{P,T}& a_y^{P,T}& a_z^{P,T}\\ b_x^{P,T}& b_y^{P,T}& b_z^{P,T}\\ c_x^{P,T}& c_y^{P,T}& c_z^{P,T}\end{array}\right]=\left[\begin{array}{lll}{a}_x^{0,0}& a_y^{0,0}& a_z^{0,0}\\ b_x^{0,0}& b_y^{0,0}& b_z^{0,0}\\ c_x^{0,0}& c_y^{0,0}& c_z^{0,0}\end{array}\right]\left(\left[\begin{array}{lll}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]+\left[\begin{array}{lll}{\alpha }_1\left(P,T\right)& {\alpha }_6\left(P,T\right)& {\alpha }_5\left(P,T\right)\\ {\alpha }_6\left(P,T\right)& {\alpha }_2\left(P,T\right)& {\alpha }_4\left(P,T\right)\\ {\alpha }_5\left(P,T\right)& {\alpha }_4\left(P,T\right)& {\alpha }_3\left(P,T\right)\end{array}\right]\right)$$E31

where $\left[\begin{array}{lll}{a}_x^{P,T}& a_y^{P,T}& a_z^{P,T}\\ b_x^{P,T}& b_y^{P,T}& b_z^{P,T}\\ c_x^{P,T}& c_y^{P,T}& c_z^{P,T}\end{array}\right]$ and $\left[\begin{array}{lll}{\alpha }_1\left(P,T\right)& {\alpha }_6\left(P,T\right)& {\alpha }_5\left(P,T\right)\\ {\alpha }_6\left(P,T\right)& {\alpha }_2\left(P,T\right)& {\alpha }_4\left(P,T\right)\\ {\alpha }_5\left(P,T\right)& {\alpha }_4\left(P,T\right)& {\alpha }_3\left(P,T\right)\end{array}\right]$ separately represent the cell configuration tensor and the deformation tensor at temperature T (K) under the pressure P (GPa).

3.1. Calculation of elastic constants for single crystal structure

A multi-particle electronic structure satisfies the Schrödinger equation. As in [29], Kohn-Sham equation as an approximation to simplify Schrödinger equation is described. For crystal composed by vibrator with the vibration frequency w_i, the total Helmholtz free energy is

$$F=E+\sum _i{F}_i^{\mathit{th}}=U+\sum _i\frac{1}{2}\mathrm{\hslash }{w}_i+k_{B}T\sum _i\mathit{\text{In}}\left(1-e^{-\frac{\mathrm{\hslash }{w}_i}{k_{B}T}}\right)$$E32

Helmholtz free energy can be calculated for all the thermodynamic quantities. DFT-QHA (quasi-harmonic approximation) is a precise calculation method to calculate thermodynamic properties of solid materials elastic constants and Debye temperature with the accurate predictions.

According to the theory of elasticity, under the isothermal strain, the elastic modulus of Helmholtz free energy can be described by the form of the Taylor expansion, of which the coefficients of the polynomial are the elastic coefficient:

$${\rho }_{0}F\left({\eta }_{\mathit{ij}},T\right)={\rho }_{0}F\left({\eta }_{\mathit{ij}},T\right)+\frac{1}{2}\sum _{\mathit{\text{ijkl}}}{c}_{\mathit{\text{ijkl}}}^T{\eta }_{\mathit{ij}}{\eta }_{\mathit{kl}}+\dots +\frac{1}{n!}\sum _{\mathit{\text{ijkl}}\dots }{c}_{\mathit{\text{ijkl}}\dots }^T{\eta }_{\mathit{ij}}{\eta }_{\mathit{kl}}\dots$$E33

where ${\eta }_{\mathit{ij}}$,${\eta }_{\mathit{kl}}$_, and ${\eta }_{\mathit{mn}}$ are the coefficients of Lagrange deformation tensor, $c_{\mathit{\text{ijkl}}}^T$is the isothermal first-order elastic coefficients, and $F\left({\eta }_{\mathit{ij}},T\right)$ is the Helmholtz free energy.

The components of the stress tensor can be extracted by ${\sigma }_i=\sum _{j=1}^6{c}_{\mathit{ij}}{\epsilon }_j$, after the applied strain, the total energy variation of the system can be expressed as

$$\mathit{\Delta E}=\frac{V}{2}\sum _{i=1}^6\sum _{j=1}^6{c}_{\mathit{ij}}{e}_i{e}_j$$E34

The second-order elastic coefficients can be obtained by the coefficient of the second-order Taylor expansion of Helmholtz free energy with the strain,

$$c_{\mathit{\text{ijkl}}}^T={\rho }_0\frac{{\partial }^{2}F\left({\eta }_{\mathit{ij}},T\right)}{\partial {\eta }_{\mathit{ij}}\partial {\eta }_{\mathit{kl}}}$$E35

Here, strain and thermodynamics deformation are symmetric. There is only six independent deformation tensor in the nine-dimensional deformation tensor. LCEC is a second-order linear combination of independent elastic coefficients corresponding to Helmholtz free energy coefficient under some deformation mode [30, 31]. For all directions under monoclinic crystals, if a strain is added, the corresponding simultaneous equations can be solved to determine all elastic coefficients.

3.2. The energy-volume relationship of the monoclinic crystal

Deformation tensors to calculate independent C_ij constants of monoclinic crystal are listed in Table 1.

Deformation tensor	ΔE-V relation of LCEC	LCEC
$e=\left(\delta ,\delta ,0,0,0,0\right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{11}}{2}+c_{12}+\frac{c_{22}}{2}\right){\delta }^2$	c11+ c22 + 2c12
$e=\left(0,\delta ,\delta ,0,0,0\right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{22}}{2}+c_{23}+\frac{c_{33}}{2}\right){\delta }^2$	c22 + c33 + 2c23
$e=\left(\delta ,0,\delta ,0,0,0\right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{11}}{2}+c_{13}+\frac{c_{33}}{2}\right){\delta }^2$	c11 + c33 + 2c13
$e=\left(0,0,0,\delta ,\delta ,0\right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{44}}{2}+c_{45}+\frac{c_{55}}{2}\right){\delta }^2$	c44 + c55 + 2c45
$e=\left(\delta ,0,0,0,0,\delta \right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{11}}{2}+c_{16}+\frac{c_{66}}{2}\right){\delta }^2$	c11 + c66 + 2c16
$e=\left(0,\delta ,0,0,0,\delta \right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{22}}{2}+c_{26}+\frac{c_{66}}{2}\right){\delta }^2$	c22 + c66 + 2c26
$e=\left(0,0,\delta ,0,0,\delta \right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{33}}{2}+c_{36}+\frac{c_{66}}{2}\right){\delta }^2$	c33 + c66 + 2c36
$e=\left(\delta ,\delta ,\delta ,0,0,0\right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{11}}{2}+\frac{c_{22}}{2}+\frac{c_{33}}{2}+c_{12}+c_{13}+c_{23}\right){\delta }^2$	c11 + c22 + c33 + 2c12 + 2c13 + 2c23
$e=\left(\delta ,-\delta ,{\delta }^2/\left(1-{\delta }^2\right),0,0,0\right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{11}}{2}-c_{12}+\frac{c_{22}}{2}\right){\delta }^2$	c11 + c22 − 2c12
$e=\left(\delta ,{\delta }^2/\left(1-{\delta }^2\right),-\delta ,0,0,0\right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{11}}{2}-c_{13}+\frac{c_{33}}{2}\right){\delta }^2$	c11 + c33 − 2c13
$e=\left({\delta }^2/\left(1-{\delta }^2\right),\delta ,-\delta ,0,0,0\right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{22}}{2}-c_{23}+\frac{c_{33}}{2}\right){\delta }^2$	c22 + c33 − 2c23
$e=\left({\delta }^2/\left(1-{\delta }^2\right),0,0,2\delta ,0,0\right)$	$\frac{\mathit{\Delta E}}{V_0}=2c_{44}{\delta }^2$	4c44
$e=\left(\delta ,0,0,0,2\delta ,0\right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{11}}{2}+2c_{15}+2c_{55}\right){\delta }^2$	c11 + 4c55 + 4c15
$e=\left(\delta ,0,0,0,-2\delta ,0\right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{11}}{2}-2c_{15}+2c_{55}\right){\delta }^2$	c11 + 4c55 − 4c15
$e=\left(0,\delta ,0,0,2\delta ,0\right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{22}}{2}+2c_{25}+2c_{55}\right){\delta }^2$	c22 + 4c55 + 4c25
$e=\left(0,0,\delta ,0,2\delta ,0\right)$	$\frac{\mathit{\Delta E}}{V_0}=\left(\frac{c_{33}}{2}+2c_{35}+2c_{55}\right){\delta }^2$	c33 + 4c55 + 4c35

Table 1.

Deformation tensors to calculate independent elastic constants of monoclinic crystal [30, 31].

For monoclinic crystal, elastic constants include C₁₁, C₂₂, C₃₃, C_l2, C₁₃, C₂₃, C₄₄, C₅₅, C₆₆ C_l5, C₂₅, C₃₅, and C₄₆; the strain energy-volume relation and elastic moduli of monoclinic symmetry based on E-V method can be obtained. The calculated E-δ points are fitted to second-order polynomials E(V, δ). For all strains, different strain forms δ are taken to calculate the total energies E for the strained crystal structure. By applying a series of δ strain amplitude, the independent elastic constants of monoclinic crystal by these simultaneous ΔE-δ equations can be obtained.

3.3. Homogenization of monoclinic polycrystals by RVH estimation

Stress-strain relation in an orthotropic monoclinic crystal can be defined by the independent elastic stiffness parameters [32]:

$$\left[\begin{array}{l}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{12}\\ {\sigma }_{13}\\ {\sigma }_{23}\end{array}\right]=\left[\begin{array}{llllll}{c}_{11}& c_{12}& c_{13}& 0& c_{15}& 0\\ c_{12}& c_{22}& c_{23}& 0& c_{25}& 0\\ c_{13}& c_{23}& c_{33}& 0& c_{35}& 0\\ 0& 0& 0& c_{44}& 0& c_{46}\\ c_{15}& c_{25}& c_{35}& 0& c_{55}& 0\\ 0& 0& 0& c_{46}& 0& c_{66}\end{array}\right]\left[\begin{array}{l}{\epsilon }_{11}\\ {\epsilon }_{22}\\ {\epsilon }_{33}\\ {\gamma }_{12}\\ {\gamma }_{13}\\ {\gamma }_{23}\end{array}\right]$$E36

where σ represents the normal stress and shear stress in each direction (unit: nN/nm²); ε and γ are the normal strain and shear strain in each direction, respectively.

The homogenized elastic properties of polycrystals can be calculated, of which elastic moduli and Poisson’s ratio can be obtained by calculating Voigt and Reuss bounds and averaging term as [32]

$$G_V=\frac{1}{15}\left[c_{11}+c_{22}+c_{33}+3\left(c_{44}+c_{55}+c_{66}\right)-\left(c_{12}+c_{13}+c_{23}\right)\right]$$E37

$$\begin{array}{l}{G}_R=15\{4[\left(c_{33}{c}_{55}-c_{35}^2\right)\left(c_{11}+c_{22}+c_{12}\right)+\left(c_{23}{c}_{55}-c_{25}{c}_{35}\right)\left(c_{11}-c_{12}-c_{23}\right)\hfill \\ +\left(c_{13}{c}_{35}-c_{15}{c}_{33}\right)\left(c_{15}+c_{25}\right)+(c_{13}{c}_{55}-c_{15}{c}_{35})\cdot \left(c_{22}-c_{12}-c_{23}-c_{13}\right)\hfill \\ +\left(c_{13}{c}_{25}-c_{15}{c}_{23}\right)\left(c_{15}-c_{25}\right)+f]/\Omega {+3\left[g/\Omega +\left(c_{44}+c_{66}\right)/\left(c_{44}{c}_{66}-c_{46}^2\right)\right]\}}^{-1}\hfill \end{array}$$E38

$$B_V=\left[c_{11}+c_{22}+c_{33}+2\left(c_{12}+c_{13}+c_{23}\right)\right]/9$$E39

$$\begin{array}{l}{B}_R=\Omega [\left(c_{33}{c}_{55}-c_{35}^2\right)\left(c_{11}+c_{22}-2c_{12}\right)+\left(c_{23}{c}_{55}-c_{25}{c}_{35}\right)\left(2c_{12}-2c_{11}-c_{23}\right)\hfill \\ +\left(c_{13}{c}_{35}-c_{15}{c}_{33}\right)\cdot \left(c_{15}-2c_{25}\right)+\left(c_{13}{c}_{55}-c_{15}{c}_{35}\right)\left(2c_{12}+2c_{23}-c_{13}-2c_{22}\right)\hfill \\ {+2\left(c_{13}{c}_{25}-c_{15}{c}_{23}\right)\left(c_{25}-c_{15}\right)+f]}^{-1}\hfill \end{array}$$E40

$$f=c_{11}\left(c_{22}{c}_{55}-c_{25}^2\right)-c_{12}\left(c_{12}{c}_{55}-c_{15}{c}_{25}\right)+c_{15}\left(c_{12}{c}_{25}-c_{15}{c}_{22}\right)+c_{25}\left(c_{23}{c}_{35}-c_{25}{c}_{33}\right)$$E41

$$g=c_{11}{c}_{22}{c}_{33}-c_{11}{c}_{23}^2-c_{22}{c}_{13}^2-c_{33}{c}_{12}^2+2c_{12}{c}_{13}{c}_{23}$$E42

$$\begin{array}{l}\mathrm{\Omega }=2\left.\left[c_{15}{c}_{25}\left(c_{33}{c}_{12}-c_{13}{c}_{23}\right)\right.+c_{15}{c}_{35}\left(c_{22}{c}_{13}-c_{12}{c}_{23}\right)+c_{25}{c}_{35}\left(c_{11}{c}_{23}-c_{12}{c}_{13}\right)\right]\\ \left.-\left[c_{15}^2\left(c_{22}{c}_{33}-c_{23}^2\right)\right.+c_{25}^2\left(c_{11}{c}_{33}-c_{13}^2\right)+c_{35}^2\left(c_{11}{c}_{22}-c_{12}^2\right)\right]+{\mathit{gc}}_{55}\end{array}$$E43

For monoclinic crystal structure, elastic constants include C₁₁, C₂₂, C₃₃, C_l2, C₁₃, C₂₃, C₄₄, C₅₅, C₆₆ C_l5, C₂₅, C₃₅, and C₄₆. The criteria for mechanical stability are given by Wu [32]:

$$c_{\mathit{ij}}>0\left(i=1,2,3,4,5,6\right)$$E44

$$\left(c_{44}{c}_{66}-c_{46}^2\right)>0$$E45

$$\left(c_{33}{c}_{55}-c_{35}^2\right)>0$$E46

$$\left(c_{22}+c_{33}-2c_{23}\right)>0$$E47

$$\left[c_{11}+c_{22}+c_{33}+2\left(c_{12}+c_{13}+c_{23}\right)\right]>0$$E48

$$\left[c_{22}\left(c_{33}{c}_{55}-c_{35}^2\right)+2c_{23}{c}_{25}{c}_{35}-c_{23}^2{c}_{55}-c_{25}^2{c}_{33}\right]>0$$E49

$$\begin{array}{l}\left\{2\left.\left[c_{15}{c}_{25}\left(c_{33}{c}_{12}-c_{13}{c}_{23}\right)\right.+c_{15}{c}_{35}\left(c_{22}{c}_{13}-c_{12}{c}_{23}\right)+c_{25}{c}_{35}\left(c_{11}{c}_{23}-c_{12}{c}_{13}\right)\right]\right.\\ \left.\left.-\left[c_{15}^2\left(c_{22}{c}_{33}-c_{23}^2\right)\right.+c_{25}^2\left(c_{11}{c}_{33}-c_{13}^2\right)+c_{35}^2\left(c_{11}{c}_{22}-c_{12}^2\right)\right]+{\mathit{gc}}_{55}\right\}>0\end{array}$$E50

Young’s modulus and Poisson’s ratio can be rewritten based on the Voigt-Reuss-Hill approximation [33]. In terms of the Voigt-Reuss-Hill approximations [34], M_H = (1/2)(M_R + M_V), M refers to B or G. Thus, Young’s modulus E and Possion’s ratio μ are obtained as

$$E=\frac{9\mathit{BG}}{3B+G}=\frac{9\left(B_V/2+B_R/2\right)\left(G_V/2+G_R/2\right)}{3\left(B_V/2+B_R/2\right)+\left(G_V/2+G_R/2\right)}$$E51

$$\mu =\frac{3B-2G}{2\left(3B+G\right)}=\frac{3\left(B_V/2+B_R/2\right)-2\left(G_V/2+G_R/2\right)}{6\left(B_V/2+B_R/2\right)+2\left(G_V/2+G_R/2\right)}$$E52

Then, Voigt-Reuss-Hill average [32] will be determined, and Young’s modulus can be calculated.

4. Modeling and homogenized elastic moduli of gypsum structure

4.1. Nanoscale modeling of monoclinic crystals

4.1.1. Nanoscale modeling of monoclinic gypsum crystal

The gypsum morphology is monoclinic, and the initial lattice is as a = 5.677Å, b = 15.207Å, c = 6.528Å, α = β = 90°, and γ = 118.49°, its structure is monoclinic with space group I 2/a [35].

In Figure 1, the gypsum crystal can be summarized as follows: (1) the two hydrogen atoms of water molecules formed weak hydrogen bonds with the O atoms of Ca and S polyhedra; (2) a stacking sequence of CaO₈ and SO₄ chains in the (010) plane alternates with water layers along the b-axis; and (3) in (010) plane, the sulfate tetrahedra and CaO₈ polyhedra alternate to form edge-sharing chains along [100] and zigzag chains along [001] direction [36] (Table 2).

Figure 1.

Modeling of gypsum crystal. (a) Gypsum structure [36] along [001]; (b) the real cell; (c) in x-direction; (d) in y-direction; and (e) in z-direction.

Atom	x	y	z	Occupancy rate	Uiso or Ueq
Ca	0.5000	0.0786	0.2500	1.00	1.00
S	0.0000	0.0787	0.7500	1.00	1.00
O1	−0.0384	0.1326	0.5512	1.00	1.00
O2	0.2429	0.0215	0.8347	1.00	1.00
Ow	0.3784	0.1825	0.4554	1.00	1.00
H1	0.2504	0.1615	0.5009	1.00	1.00
H2	0.4022	0.2435	0.4900	1.00	1.00

Table 2.

Atomic coordinates and displacement parameters of gypsum [36].

4.1.2. Nanoscale modeling of monoclinic 11 Å tobermorite crystal

Hamid model [37] as the 11 Å tobermorite (formula: Ca₄Si₆O₁₄(OH)₄·2H₂O) as an initial configuration is commonly used. The morphology is monoclinic, and the initial lattice is [37]: a = 6.69 Å, b = 7.39 Å, c = 22.779 Å, α = β = 90°, and γ = 123.49°, space group P21. Modeling of 11 Å tobermorite is shown in Figure 2.

Figure 2.

Modeling of 11 Å tobermorite crystal. Silicate chains, calcium octahedral, and oxygen atoms are shown as yellow tetrahedra, green spheres, and red spheres. (a) 11 Å Tob monoclinic crystal; (b) in x-direction; (c) in y-direction; and (d) in z-direction.

In Figure 2(a), the 11 Å tobermorite crystal can be summarized as follows: (1) the structure is basically a layered structure. (2) The central part is a Ca-O sheet (with an empirical formula: CaO₂, of which the oxygen in CaO₂ also includes that of the silicate tetrahedron part). (3) Silicate chains envelope the Ca-O sheet on both sides. (4) Ca²⁺ and H₂O are filled between individual layers to balance the charges. The infinite layers of calcium polyhedra are parallel to (001), with tetrahedral chains of wollastonite-type along b and the composite layers stacked along c and connected through the formation of double tetrahedral chains [38]. Atomic coordinates and displacement parameters are seen in Table 3.

Atomic species	X	Y	Z	Occupancy rate	Uiso or Ueq	Atomic species	X	Y	Z	Occupancy rate	Uiso or Ueq
Si1	0.7710	0.3830	0.1578	1	0.031	O8	0.7690	0.8430	0.0953	1	0.027
Si2	0.9250	0.7500	0.0721	1	0.030	O9	0.5370	0.7980	0.1968	1	0.036
Si3	0.7720	0.9620	0.1596	1	0.015	O10	0.0040	0.0420	0.2008	1	0.034
O1	0.7740	0.4950	0.0932	1	0.039	O11	0.4330	0.2230	−0.0250	0.5	0.072
O2	0.7620	0.1690	0.1305	1	0.019	O12	0.9490	0.2560	0.0000	1	0.080
O3	0.0020	0.5270	0.2000	1	0.032	O13	0.4300	0.7700	−0.0220	0.5	0.090
O4	0.5360	0.3040	0.1926	1	0.035	Cal	0.2770	0.4257	0.2083	1	0.024
O5	0.9100	0.7470	0.0000	1	0.034	Ca2	0.7630	0.9160	0.2951	1	0.027
O6	0.2020	0.8870	0.0942	1	0.053	Ca3	0.5620	0.0640	0.0450	0.25	0.038
O7	0.2890	0.4360	0.0940	1	0.076	—	—	—	—	—	—

Table 3.

Atomic coordinates and displacement parameters of 11 Å tobermorite [38].

4.2. Initial conditions and elastic constants of monoclinic crystals

4.2.1. Initial conditions and elastic constants of gypsum

The initial conditions are as follows: the pressure region of 0–1 GPa is used. Besides, a plane-wave basis set and ultrasoft pseudopotentials using GGA are used with a plane-wave cutoff energy of 400 eV. Brillouin zone is 6 × 6 × 4. Self-consistent convergence of the total energy per atom is chosen as 10⁻⁴ eV. Elastic constants of monoclinic gypsum crystal under 0–1.0 GPa are shown in Figure 3.

Figure 3.

Gypsum monoclinic crystal under pressure 0–1.0 GPa by DFT. (a) Relative change of a, b, c, and V and (b) elastic constants.

From Figure 3, elastic constants at 0 GPa are given as c₁₁ = 82.464 GPa, c₁₂ = 34.751 GPa, c₁₃ = 33.643 GPa, c₁₅ = −1.987 GPa, c₂₂ = 63.046 GPa, c₂₃ = 34.920 GPa, c₂₅ = −8.071 GPa, c₃₃ = 57.549 GPa, c₃₅ = −3.054 GPa, c₄₄ = 20.863 GPa, c₄₆ = −4.688 GPa, c₅₅ = 28.062 GPa, and c₆₆ = 28.556 GPa. It is found that the oxygen atom of the water molecule did not change its position or occupancy under pressure conditions. A simple pressure increase at an ambient temperature cannot induce dehydration because of the unchange of water molecular in the gypsum structure within pressure range [36].

Elastic constants of gypsum crystal model based on DFT are calculated, and parameters are detailed in Table 4.

P	C11	C12	C13	C15	C22	C23	C25	C33	C35	C44	C46	C55	C66
10^–4[36^]	—	—	—	—	—	—	—	—	—	—	—	—	—
0.0	82.46	34.75	33.64	−1.99	63.05	34.92	−8.07	57.55	−3.05	20.86	−4.69	28.06	28.56
0.1	79.82	32.64	29.2	1.8	71.04	29.61	−7.54	61.88	−3.22	20.13	−3.06	26.19	27.7
0.2	82.93	37.75	34.59	1.09	63.62	32.42	−7.35	50.64	−4.37	21.32	−1.1	25.8	17.8
0.3	82.82	39.77	32.81	0.17	65.64	29.61	−7.23	57.31	−4.45	26.43	−5.57	23.17	23.39
0.4	84.47	38.6	32.25	1.27	69.03	32.31	−8.51	53.41	−2.21	20.8	−2.03	28.41	22.34
0.5	75.84	43.68	29.39	0.57	68.7	33.18	−8.36	56.08	−2.52	29.7	−3.88	27.35	22.4
0.6	74.22	43.11	28.77	2.22	69.52	28.87	−7.81	53.19	−2.68	28.97	−1.49	23.24	15.53
0.7	88.37	41.74	32.85	2.25	70.09	32.28	−9.14	55.48	−4.28	24.66	−2.76	27.25	22.58
0.8	88.53	39.65	35.29	2.96	73.28	33.84	−8.02	62.22	−3.73	24.73	−3.44	26.37	24.39
0.9	88.7	45.09	37.97	4.54	66.78	36.02	−10.4	61.98	−1.2	25.15	−4.92	28.93	26.82
1.0	90.12	39.79	34.63	2.7	75.99	34.92	−8.74	68.31	−3.46	26.32	−5.83	28.15	30.19

Table 4.

Elastic coefficient C_ij (GPa) of gypsum by DFT.

4.2.2. Initial conditions and elastic constants of tobermorite

Initial conditions of tobermorite are quite the same with that of gypsum crystal. Elastic constants of 11 Å tobermorite crystal under 0–1.0 GPa are shown in Figure 4. Elastic constants are shown in Table 5.

Figure 4.

11 Å tobermorite monoclinic crystal under pressure 0–1.0 GPa by DFT. (a) Relative change of a, b, c, and V and (b) elastic constants.

P/GPa	C11	C12	C13	C15	C22	C23	C25	C33	C35	C44	C46	C55	C66
SHA^[39^]	102.65	41.68	27.70	1.25	125.05	18.83	−4.10	83.80	−3.38	22.90	−11.93	23.25	50.20
0.0	106.63	50.37	41.09	−3.50	131.67	22.78	−0.78	71.45	−0.83	26.03	−0.02	27.61	45.26
0.1	118.37	45.40	35.91	−3.52	129.18	17.19	0.11	67.84	−0.55	32.51	3.90	32.74	40.07
0.2	109.13	45.84	35.63	−3.22	136.79	23.05	0.03	82.75	0.06	28.88	1.21	22.40	45.69
0.3	115.53	46.36	40.17	−4.46	142.59	27.65	−0.04	95.03	0.02	31.08	0.49	32.38	50.57
0.4	102.65	35.38	38.73	−6.32	123.43	18.11	−1.92	74.28	0.05	18.14	−0.83	22.38	40.44
0.5	100.08	42.58	36.10	−4.52	137.56	21.68	−0.26	90.87	0.36	29.92	−0.46	29.66	51.82
0.6	97.87	44.09	28.76	−5.85	162.17	25.77	0.19	93.71	−0.14	24.89	−1.20	26.63	40.26
0.7	108.73	48.60	34.07	−4.55	147.09	26.78	−0.14	92.64	−0.01	21.56	2.06	44.25	41.23
0.8	122.87	55.30	40.62	−4.05	155.75	29.54	−0.25	103.3	−0.57	24.90	0.72	33.31	42.67
0.9	120.77	44.19	45.41	−4.82	139.59	13.68	0.09	88.25	−0.19	27.18	−0.35	26.85	53.22
1.0	127.01	41.78	45.00	−4.47	143.72	23.65	−0.02	98.30	−0.12	29.98	0.71	32.08	45.68

Table 5.

Elastic coefficient C_ij (GPa) of 11 Å tobermorite by DFT.

A comparisonal results of Shahsavari [39] are provided. Elastic constants at 0 GPa are as follows: c₁₁ = 106.63 GPa, c₁₂ = 50.37 GPa, c₁₃ = 41.09 GPa, c₁₅ = −3.50 GPa, c₂₂ = 131.67 GPa, c₂₃ = 22.78 GPa, c₂₅ = −0.78 GPa, c₃₃ = 71.45 GPa, c₃₅ = −0.83 GPa, c₄₄ = 26.03 GPa, c₄₆ = −0.02 GPa, c₅₅ = 27.61 GPa, and c₆₆ = 45.26 GPa. Thus, elastic modulus can be homogenized to compare with the results of LD C-S-H phase in nano-indentation test by Vandamme and Ulm [40].

4.3. Homogenized elastic moduli of typical monoclinic structures

4.3.1. Elastic modulus of monoclinic gypsum structure

Based on elastic constants, the elastic moduli of gypsum at 0 GPa are verified and averaged in Figure 5.

Figure 5.

Elastic moduli of gypsum crystal under pressure 0–1.0 GPa.

As gypsum shows anisotropic compressibility along three crystallographic axes with b > c > a below 5 GPa [44], the pressure region of 0–1.0 GPa is used to verify whether the performance of model under low pressure is stable. Mechanical moduli of gypsum polycrystalline are listed in Table 6.

Pressure (GPa)	Gv (GPa)	Bv (GPa)	Gr (GPa)	Br (GPa)	B (GPa)	G (GPa)	E (GPa)	μ
Reference [43]	26.5333	39.2556	24.8077	39.2381	25.6705	39.2469	63.2265	0.2315
0.0	22.1459	45.5208	19.7054	43.8224	20.9257	44.6716	54.2985	0.2974
0.1	22.8896	43.9624	21.7569	42.9250	22.3233	43.4437	57.1766	0.2806
0.2	19.1472	45.1906	17.4395	41.9064	18.2934	43.5485	48.1394	0.3158
0.3	21.5041	45.5718	19.3501	43.6675	20.4271	44.6197	53.1678	0.3014
0.4	21.2263	45.9146	19.5324	43.2463	20.3794	44.5805	53.0537	0.3017
0.5	22.1809	45.9026	19.4980	43.7566	20.8395	44.8296	54.1306	0.2988
0.6	19.9615	44.2709	17.7554	41.5818	18.8585	42.9264	49.3486	0.3084
0.7	22.0356	47.5189	20.1833	44.0623	21.1095	45.7906	54.8931	0.3002
0.8	22.7817	49.0659	21.3745	47.0681	22.0781	48.0670	57.4399	0.3008
0.9	22.7399	50.6242	19.4094	48.1542	21.0747	49.3892	55.3510	0.3132
1.0	25.2707	50.3430	23.4785	48.9327	24.3746	49.6379	62.8382	0.2890

Table 6.

Mechanical moduli of gypsum polycrystalline by different methods.

As an acoustic method [41] and mechanical properties [42] have been investigated, according to elastic constants of gypsum crystal [43], elastic moduli by experiment can be calculated, as shown in Table 6. Elastic moduli are as follows: G_v = 22.146 GPa, G_r = 19.705 GPa, B_v = 45.521 GPa, B_r = 43.822 GPa, B = 44.672 GPa, G = 20.926 GPa, E = 54.299 GPa, and μ = 0.2974. These results are close to the plane-strain value of Young’s modulus by reference [44] E = 50 GPa, μ = 0.45. By comparison of gypsum crystal and CH crystal, axial moduli of gypsum in x, y, and z directions are 57.75, 37.22, and 34.91 GPa, while axial moduli of Ca(OH)₂ in x, y, and z directions are 93.75, 93.75, and 42.39 GPa, showing that gypsum crystal is much less anisotropic than hydrogen-bonded layered Ca(OH)2 structure [42].

4.3.2. Elastic modulus of monoclinic tobermorite structure

Based on elastic constants of 11 Å tobermorite crystal using GGA calculation method by DFT, bulk modulus B and shear modulus G are separately calculated by Eqs. (37)–(50) (Figure 6).

Figure 6.

Elastic moduli of 11 Å tobermorite crystal under pressure 0–1.0 GPa.

Elastic moduli at 0 GPa are verified and averaged as G_v = 32.815 GPa, B_v = 59.803 GPa, G_r = 29.908 GPa, B_r = 54.276 GPa, E = 79.512 GPa, and μ = 0.268. Young’s modulus is about 79.512 GPa by Reuss-Voigt-Hill estimation, which is close to the simulation result of 89 GPa [45] by Pellenq and result of 78.939 GPa [39] by Shahsavari. Mechanical moduli by different methods are listed in Table 7.

Pressure (GPa)	Bv (GPa)	Br (GPa)	Gv (GPa)	Gr (GPa)	B (GPa)	G (GPa)	E (GPa)	μ
Reference [191]	54.2133	51.6976	34.1560	28.9168	52.9555	31.5364	78.9391	0.2516
0.0	59.8066	54.2778	32.8140	29.9063	57.0399	31.3615	79.5121	0.2677
0.1	56.9306	50.1535	35.5205	33.4468	53.5438	34.4843	85.1689	0.2349
0.2	59.7454	56.4236	34.3364	31.4164	58.0846	32.8763	82.9743	0.2619
0.3	64.6138	62.4671	38.7385	36.7029	63.5388	37.7202	94.4669	0.2522
0.4	53.8665	50.9999	30.0678	26.1841	52.4333	28.1257	71.5786	0.2725
0.5	58.8039	56.9635	37.4898	34.7400	57.8831	36.1145	89.6903	0.2417
0.6	61.2236	57.0694	35.3648	32.2837	59.1442	33.8242	85.2259	0.2598
0.7	63.0391	60.0785	37.3432	34.0050	61.5598	35.6744	89.6966	0.2572
0.8	70.3213	67.4998	37.2739	34.9015	68.9022	36.0849	92.1653	0.2771
0.9	61.6837	57.8325	37.8067	33.5164	59.7599	35.6606	89.2325	0.2511
1.0	65.5442	63.4872	38.6855	36.6712	64.5147	37.7292	94.7226	0.2553

Table 7.

Mechanical moduli of 11Å tobermorite polycrystalline by different methods.

However, these values considering the ordered Si-chain at a long range are far away from the nano-indentation experiment performed on the C-S-H phase [40, 46]. It confirms the absence of order at a long range in this phase and that the up-scaling to polycrystals cannot be done with the tobermorite model. Modeling of C-S-H structure with disordered Si chain should be fairly considered.

5. Conclusions

Elastic constants of gypsum and tobermorite structures under a certain pressure region are calculated by DFT method, which has a certain value for both application and reference. Results are as follows:

1. For monoclinic gypsum and tobermorite crystals, elastic coefficients are obtained in 0–1-GPa pressure range to verify the reliability of the model by comparing other literatures.

2. Elastic constants of gypsum single crystal at 0 GPa are given as follows: c₁₁ = 82.464 GPa, c₁₂ = 34.751 GPa, c₁₃ = 33.643 GPa, c₁₅ = −1.987 GPa, c₂₂ = 63.046 GPa, c₂₃ = 34.920 GPa, c₂₅ = −8.071 GPa, c₃₃ = 57.549 GPa, c₃₅ = −3.054 GPa, c₄₄ = 20.863 GPa, c₄₆ = −4.688 GPa, c₅₅ = 28.062 GPa, and c₆₆ = 28.556 GPa.

3. Elastic constants of 11Å tobermorite single crystal at 0 GPa are as follows: c₁₁ = 106.63 GPa, c₁₂ = 50.37 GPa, c₁₃ = 41.09 GPa, c₁₅ = −3.50 GPa, c₂₂ = 131.67 GPa, c₂₃ = 22.78 GPa, c₂₅ = −0.78 GPa, c₃₃ = 71.45 GPa, c₃₅ = −0.83 GPa, c₄₄ = 26.03 GPa, c₄₆ = −0.02 GPa, c₅₅ = 27.61 GPa, and c₆₆ = 45.26 GPa.

4. Young’s modulus of gypsum is about 54.299 GPa. Elastic moduli at 0 GPa are as follows: G_v = 22.146 GPa, G_r = 19.705 GPa, B_v = 45.521 GPa, B_r = 43.822 GPa, E = 54.299 GPa, and μ = 0.297.

5. Young’s modulus of 11Å tobermorite is about 79.512 GPa. Elastic moduli at 0 GPa are as follows: G_v = 32.815 GPa, B_v = 59.803 GPa, G_r = 29.908 GPa, B_r = 54.276 GPa, E = 79.512 GPa, and μ = 0.268.

Structural, elastic properties of monoclinic crystals are investigated, and C_ij determination is given by DFT method. Reuss-Voigt-Hill estimation has been used for polycrystal structures and can be seen as an intermediate step in the homogenization of elastic properties.

Acknowledgments

The authors greatly acknowledge the financial support for this work provided by the China Scholarship Council (CSC) and the support of start-up foundation of Xi’an Shiyou University. Thanks to Qiufeng Wang for her proofreading.