Open access peer-reviewed chapter

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

By Jia Fu

Submitted: June 14th 2017Reviewed: November 8th 2017Published: December 20th 2017

DOI: 10.5772/intechopen.72301

## 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.

### Keywords

• DFT calculation
• single crystal
• nano scale
• elastic constants
• homogenized moduli

## 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 CaSO4•2H2O, which is an evaporite mineral and a kind of hydration product of anhydrite (chemical formula: CaSO4). Moreover, the 11 Å tobermorite model (chemical formula: Ca4Si6O14(OH)4·2H2O) 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 positionedandrespectively.

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 n0(r). Besides, the scalar function of the position n0(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 asand 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 potentialand 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

EHKn=Tn+Eintn+d3rVexr+EnnRFHKn+d3rVexr+EnnRE9

As a result, the density n0(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

EBOR=minERnrE10

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 FHK(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 EXCn.

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
nr=i=1Neφir2E14

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 TS, all other terms depend on the density. Therefore, it is possible to vary the functions of the wave and to derive the variational equation:

δEKSδφir=δTSδφir+δEexδnr+δEHartreeδnr+δExcδnrδnrδφir=0E15

With the constraint of orthonormalization φiφj=δi,j, this implies the form of Kohn-Sham for Schrödinger equations:

HKSεiφir=0E16

εirepresents the eigenvalues, and HKSis the effective Hamiltonian H,

HKSr=122+VKSrE17
VKSr=Vexr+δEHartreeδnr+δExcδnrE18

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

Excn=FHKnTSn+EHartreenE19

or more precisely,

Excn=TTSn+VintEHartreenE20

This energy is related to potential exchange-correlation Vxc=Excnr.

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 EKS can be written as

Excn=nrεxcnrd3rE21

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 ExcLDAnis 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

ExcLDAn=nrεxcbomnrd3r=nrεxbomnr+εcbomnrd3rE22

The exchange term Exbomnrcan 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]:

ExcGGAn=nrεxcnnd3r=nrεxbomnFxcnnd3rE23

where εxbomis the exchange energy of an unpolarized density n(r) system. There are many forms of FXC, 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:

Heffrψir=22m2+Veffrψir=εiψirE24

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:

Veffr=mVeffGmexpiGmrE25

Gm is the reciprocal lattice vector:

VeffG=1ΩsellΩsellVeffrexpiGmrdrE26

where Ωsellis 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:

ψi,kr=expikrui,krui,kr+R=ui,krR=niai,ni=1NiE27

where R is the vector of direct space defined by ai with i123and Ni is the number of primitive cells in each direction (Niin the case of perfect crystal).

Solving Eq. (24) is equivalent to increase the periodic function ui,kr, in a database-dependent functions points k:ϕjkrj=1Nbask:

ui,kr=jCijϕjkrE28

where ϕjkis 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 ϕjkwill approximate. That the selected database simply solves the system:

m'Hm,m'kCi,m'k=εikCi,mkHm,m'k=φm,kjHeffφm,kjE29

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 fi(k) where i defines the band index, the average value is

fi¯=1NkkfikΩcell2πdBZfikdkE30

Ωcellis the cell volume of the original mesh in the real space and 2πd/Ωcellof 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

XP,T=axP,TayP,TazP,TbxP,TbyP,TbzP,TcxP,TcyP,TczP,T=ax0,0ay0,0az0,0bx0,0by0,0bz0,0cx0,0cy0,0cz0,0100010001+α1PTα6PTα5PTα6PTα2PTα4PTα5PTα4PTα3PTE31

where axP,TayP,TazP,TbxP,TbyP,TbzP,TcxP,TcyP,TczP,Tand α1PTα6PTα5PTα6PTα2PTα4PTα5PTα4PTα3PTseparately 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 wi, the total Helmholtz free energy is

F=E+iFith=U+i12wi+kBTiIn1ewikBTE32

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:

ρ0FηijT=ρ0FηijT+12ijklcijklTηijηkl++1n!ijklcijklTηijηklE33

where ηij,ηkl, and ηmnare the coefficients of Lagrange deformation tensor, cijklTis the isothermal first-order elastic coefficients, and FηijTis the Helmholtz free energy.

The components of the stress tensor can be extracted by σi=j=16cijεj, after the applied strain, the total energy variation of the system can be expressed as

ΔE=V2i=16j=16cijeiejE34

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

cijklT=ρ02FηijTηijηklE35

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 Cij constants of monoclinic crystal are listed in Table 1.

Deformation tensorΔE-V relation of LCECLCEC
e=δδ0000ΔEV0=c112+c12+c222δ2c11+ c22 + 2c12
e=0δδ000ΔEV0=c222+c23+c332δ2c22 + c33 + 2c23
e=δ0δ000ΔEV0=c112+c13+c332δ2c11 + c33 + 2c13
e=000δδ0ΔEV0=c442+c45+c552δ2c44 + c55 + 2c45
e=δ0000δΔEV0=c112+c16+c662δ2c11 + c66 + 2c16
e=0δ000δΔEV0=c222+c26+c662δ2c22 + c66 + 2c26
e=00δ00δΔEV0=c332+c36+c662δ2c33 + c66 + 2c36
e=δδδ000ΔEV0=c112+c222+c332+c12+c13+c23δ2c11 + c22 + c33 + 2c12 + 2c13 + 2c23
e=δδδ2/1δ2000ΔEV0=c112c12+c222δ2c11 + c22 − 2c12
e=δδ2/1δ2δ000ΔEV0=c112c13+c332δ2c11 + c33 − 2c13
e=δ2/1δ2δδ000ΔEV0=c222c23+c332δ2c22 + c33 − 2c23
e=δ2/1δ2002δ00ΔEV0=2c44δ24c44
e=δ0002δ0ΔEV0=c112+2c15+2c55δ2c11 + 4c55 + 4c15
e=δ0002δ0ΔEV0=c1122c15+2c55δ2c11 + 4c55 − 4c15
e=0δ002δ0ΔEV0=c222+2c25+2c55δ2c22 + 4c55 + 4c25
e=00δ02δ0ΔEV0=c332+2c35+2c55δ2c33 + 4c55 + 4c35

### Table 1.

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

For monoclinic crystal, elastic constants include C11, C22, C33, Cl2, C13, C23, C44, C55, C66 Cl5, C25, C35, and C46; 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]:

σ11σ22σ33σ12σ13σ23=c11c12c130c150c12c22c230c250c13c23c330c350000c440c46c15c25c350c550000c460c66ε11ε22ε33γ12γ13γ23E36

where σ represents the normal stress and shear stress in each direction (unit: nN/nm2); ε 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]

GV=115c11+c22+c33+3c44+c55+c66c12+c13+c23E37
GR=15{4[(c33c55c352)(c11+c22+c12)+(c23c55c25c35)(c11c12c23)+(c13c35c15c33)(c15+c25)+(c13c55c15c35)(c22c12c23c13)+(c13c25c15c23)(c15c25)+f]/Ω+3[g/Ω+(c44+c66)/(c44c66c462)]}1E38
BV=c11+c22+c33+2c12+c13+c23/9E39
BR=Ω[(c33c55c352)(c11+c222c12)+(c23c55c25c35)(2c122c11c23)+(c13c35c15c33)(c152c25)+(c13c55c15c35)(2c12+2c23c132c22)+2(c13c25c15c23)(c25c15)+f]1E40
f=c11c22c55c252c12c12c55c15c25+c15c12c25c15c22+c25c23c35c25c33E41
g=c11c22c33c11c232c22c132c33c122+2c12c13c23E42
Ω=2c15c25c33c12c13c23+c15c35c22c13c12c23+c25c35c11c23c12c13c152c22c33c232+c252c11c33c132+c352c11c22c122+gc55E43

For monoclinic crystal structure, elastic constants include C11, C22, C33, Cl2, C13, C23, C44, C55, C66 Cl5, C25, C35, and C46. The criteria for mechanical stability are given by Wu [32]:

cij>0i=123456E44
c44c66c462>0E45
c33c55c352>0E46
c22+c332c23>0E47
c11+c22+c33+2c12+c13+c23>0E48
c22c33c55c352+2c23c25c35c232c55c252c33>0E49
2c15c25c33c12c13c23+c15c35c22c13c12c23+c25c35c11c23c12c13c152c22c33c232+c252c11c33c132+c352c11c22c122+gc55>0E50

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], MH = (1/2)(MR + MV), M refers to B or G. Thus, Young’s modulus E and Possion’s ratio μ are obtained as

E=9BG3B+G=9BV/2+BR/2GV/2+GR/23BV/2+BR/2+GV/2+GR/2E51
μ=3B2G23B+G=3BV/2+BR/22GV/2+GR/26BV/2+BR/2+2GV/2+GR/2E52

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 CaO8 and SO4 chains in the (010) plane alternates with water layers along the b-axis; and (3) in (010) plane, the sulfate tetrahedra and CaO8 polyhedra alternate to form edge-sharing chains along [100] and zigzag chains along [001] direction [36] (Table 2).

AtomxyzOccupancy rateUiso or Ueq
Ca0.50000.07860.25001.001.00
S0.00000.07870.75001.001.00
O1−0.03840.13260.55121.001.00
O20.24290.02150.83471.001.00
Ow0.37840.18250.45541.001.00
H10.25040.16150.50091.001.00
H20.40220.24350.49001.001.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: Ca4Si6O14(OH)4·2H2O) 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.

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: CaO2, of which the oxygen in CaO2 also includes that of the silicate tetrahedron part). (3) Silicate chains envelope the Ca-O sheet on both sides. (4) Ca2+ and H2O 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 speciesXYZOccupancy rateUiso or UeqAtomic speciesXYZOccupancy rateUiso or Ueq
Si10.77100.38300.157810.031O80.76900.84300.095310.027
Si20.92500.75000.072110.030O90.53700.79800.196810.036
Si30.77200.96200.159610.015O100.00400.04200.200810.034
O10.77400.49500.093210.039O110.43300.2230−0.02500.50.072
O20.76200.16900.130510.019O120.94900.25600.000010.080
O30.00200.52700.200010.032O130.43000.7700−0.02200.50.090
O40.53600.30400.192610.035Cal0.27700.42570.208310.024
O50.91000.74700.000010.034Ca20.76300.91600.295110.027
O60.20200.88700.094210.053Ca30.56200.06400.04500.250.038
O70.28900.43600.094010.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−4 eV. Elastic constants of monoclinic gypsum crystal under 0–1.0 GPa are shown in Figure 3.

From Figure 3, elastic constants at 0 GPa are given as c11 = 82.464 GPa, c12 = 34.751 GPa, c13 = 33.643 GPa, c15 = −1.987 GPa, c22 = 63.046 GPa, c23 = 34.920 GPa, c25 = −8.071 GPa, c33 = 57.549 GPa, c35 = −3.054 GPa, c44 = 20.863 GPa, c46 = −4.688 GPa, c55 = 28.062 GPa, and c66 = 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.

PC11C12C13C15C22C23C25C33C35C44C46C55C66
10–4[36]
0.082.4634.7533.64−1.9963.0534.92−8.0757.55−3.0520.86−4.6928.0628.56
0.179.8232.6429.21.871.0429.61−7.5461.88−3.2220.13−3.0626.1927.7
0.282.9337.7534.591.0963.6232.42−7.3550.64−4.3721.32−1.125.817.8
0.382.8239.7732.810.1765.6429.61−7.2357.31−4.4526.43−5.5723.1723.39
0.484.4738.632.251.2769.0332.31−8.5153.41−2.2120.8−2.0328.4122.34
0.575.8443.6829.390.5768.733.18−8.3656.08−2.5229.7−3.8827.3522.4
0.674.2243.1128.772.2269.5228.87−7.8153.19−2.6828.97−1.4923.2415.53
0.788.3741.7432.852.2570.0932.28−9.1455.48−4.2824.66−2.7627.2522.58
0.888.5339.6535.292.9673.2833.84−8.0262.22−3.7324.73−3.4426.3724.39
0.988.745.0937.974.5466.7836.02−10.461.98−1.225.15−4.9228.9326.82
1.090.1239.7934.632.775.9934.92−8.7468.31−3.4626.32−5.8328.1530.19

### Table 4.

Elastic coefficient Cij (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.

P/GPaC11C12C13C15C22C23C25C33C35C44C46C55C66
SHA[39]102.6541.6827.701.25125.0518.83−4.1083.80−3.3822.90−11.9323.2550.20
0.0106.6350.3741.09−3.50131.6722.78−0.7871.45−0.8326.03−0.0227.6145.26
0.1118.3745.4035.91−3.52129.1817.190.1167.84−0.5532.513.9032.7440.07
0.2109.1345.8435.63−3.22136.7923.050.0382.750.0628.881.2122.4045.69
0.3115.5346.3640.17−4.46142.5927.65−0.0495.030.0231.080.4932.3850.57
0.4102.6535.3838.73−6.32123.4318.11−1.9274.280.0518.14−0.8322.3840.44
0.5100.0842.5836.10−4.52137.5621.68−0.2690.870.3629.92−0.4629.6651.82
0.697.8744.0928.76−5.85162.1725.770.1993.71−0.1424.89−1.2026.6340.26
0.7108.7348.6034.07−4.55147.0926.78−0.1492.64−0.0121.562.0644.2541.23
0.8122.8755.3040.62−4.05155.7529.54−0.25103.3−0.5724.900.7233.3142.67
0.9120.7744.1945.41−4.82139.5913.680.0988.25−0.1927.18−0.3526.8553.22
1.0127.0141.7845.00−4.47143.7223.65−0.0298.30−0.1229.980.7132.0845.68

### Table 5.

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

A comparisonal results of Shahsavari [39] are provided. Elastic constants at 0 GPa are as follows: c11 = 106.63 GPa, c12 = 50.37 GPa, c13 = 41.09 GPa, c15 = −3.50 GPa, c22 = 131.67 GPa, c23 = 22.78 GPa, c25 = −0.78 GPa, c33 = 71.45 GPa, c35 = −0.83 GPa, c44 = 26.03 GPa, c46 = −0.02 GPa, c55 = 27.61 GPa, and c66 = 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.

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.533339.255624.807739.238125.670539.246963.22650.2315
0.022.145945.520819.705443.822420.925744.671654.29850.2974
0.122.889643.962421.756942.925022.323343.443757.17660.2806
0.219.147245.190617.439541.906418.293443.548548.13940.3158
0.321.504145.571819.350143.667520.427144.619753.16780.3014
0.421.226345.914619.532443.246320.379444.580553.05370.3017
0.522.180945.902619.498043.756620.839544.829654.13060.2988
0.619.961544.270917.755441.581818.858542.926449.34860.3084
0.722.035647.518920.183344.062321.109545.790654.89310.3002
0.822.781749.065921.374547.068122.078148.067057.43990.3008
0.922.739950.624219.409448.154221.074749.389255.35100.3132
1.025.270750.343023.478548.932724.374649.637962.83820.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: Gv = 22.146 GPa, Gr = 19.705 GPa, Bv = 45.521 GPa, Br = 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)2 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).

Elastic moduli at 0 GPa are verified and averaged as Gv = 32.815 GPa, Bv = 59.803 GPa, Gr = 29.908 GPa, Br = 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.213351.697634.156028.916852.955531.536478.93910.2516
0.059.806654.277832.814029.906357.039931.361579.51210.2677
0.156.930650.153535.520533.446853.543834.484385.16890.2349
0.259.745456.423634.336431.416458.084632.876382.97430.2619
0.364.613862.467138.738536.702963.538837.720294.46690.2522
0.453.866550.999930.067826.184152.433328.125771.57860.2725
0.558.803956.963537.489834.740057.883136.114589.69030.2417
0.661.223657.069435.364832.283759.144233.824285.22590.2598
0.763.039160.078537.343234.005061.559835.674489.69660.2572
0.870.321367.499837.273934.901568.902236.084992.16530.2771
0.961.683757.832537.806733.516459.759935.660689.23250.2511
1.065.544263.487238.685536.671264.514737.729294.72260.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: c11 = 82.464 GPa, c12 = 34.751 GPa, c13 = 33.643 GPa, c15 = −1.987 GPa, c22 = 63.046 GPa, c23 = 34.920 GPa, c25 = −8.071 GPa, c33 = 57.549 GPa, c35 = −3.054 GPa, c44 = 20.863 GPa, c46 = −4.688 GPa, c55 = 28.062 GPa, and c66 = 28.556 GPa.

3. Elastic constants of 11Å tobermorite single crystal at 0 GPa are as follows: c11 = 106.63 GPa, c12 = 50.37 GPa, c13 = 41.09 GPa, c15 = −3.50 GPa, c22 = 131.67 GPa, c23 = 22.78 GPa, c25 = −0.78 GPa, c33 = 71.45 GPa, c35 = −0.83 GPa, c44 = 26.03 GPa, c46 = −0.02 GPa, c55 = 27.61 GPa, and c66 = 45.26 GPa.

4. Young’s modulus of gypsum is about 54.299 GPa. Elastic moduli at 0 GPa are as follows: Gv = 22.146 GPa, Gr = 19.705 GPa, Bv = 45.521 GPa, Br = 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: Gv = 32.815 GPa, Bv = 59.803 GPa, Gr = 29.908 GPa, Br = 54.276 GPa, E = 79.512 GPa, and μ = 0.268.

Structural, elastic properties of monoclinic crystals are investigated, and Cij 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.

chapter PDF
Citations in RIS format
Citations in bibtex format

## More

© 2017 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## How to cite and reference

### Cite this chapter Copy to clipboard

Jia Fu (December 20th 2017). Elastic Constants and Homogenized Moduli of Monoclinic Structures Based on Density Functional Theory, Density Functional Calculations - Recent Progresses of Theory and Application, Gang Yang, IntechOpen, DOI: 10.5772/intechopen.72301. Available from:

### Related Content

Next chapter

#### Application of Density Functional Theory in Soil Science

By Jiena Yun, Qian Wang, Chang Zhu and Gang Yang

First chapter

#### Numerical Solution of Linear Ordinary Differential Equations in Quantum Chemistry by Spectral Method

By Masoud Saravi and Seyedeh-Razieh Mirrajei

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

View all Books