Atomistic Simulation of Anistropic Crystal Structures at Nanoscale

Overview of Nanostructure Modeling and Atomistic Simulation Methods

Abstract

The development of molecular modeling tools to describe and predict the mechanical properties of structural materials reveals an undeniable practical importance. It is now well recognized that such an objective can be achieved through the linking of the structure of materials at the nanoscale or with their performances. At the nanometer scale, anisotropic materials exhibit differences due to the directional arrangements of atomic structure, the force between atoms, such as van der Waals force, typical Coulomb force, and other forces. Nanoscale modeling and mechanical properties by using the density functional theory (DFT), a so-called atomic finite element method (AFEM), and the classical molecular dynamics (MD) method are especially concerned according to the modeling requirement of different crystal structures investigated. Omitting the description of these structures’ importance, elastic moduli are separately calculated by either homogenization or curve fitting of the linear portion of the stress-strain curve by the corresponding numerical simulations. This chapter is committed to introduce the modeling and simulation for calculating mechanical properties (Young’s modulus especially) of typical anisotropic crystal structures using three methods (DFT, AFEM, and MD) mentioned above. It is therefore asked to connect to the nanoscale modeling and continuous pattern of behavior by identifying the relevant output data at small scales and bringing the necessary information to higher scales.

Keywords: nanoscale, crystal structures, anisotropic elasticity, DFT, AFEM, molecular dynamics simulation, nanoindentation, mechanical properties

1. Introduction

Nanotechnology is considered “high-tech” since it requires deep knowledge of the system considered at the nanoscale in order to intelligently design. Nanomaterials refer to one dimension in the nanometer range (nanosized particles, atomic clusters etc.) or a basic unit at least among a three-dimensional space (nanotubes, ultrathin films, multilayers, superlattices, etc.). Nanoscience and technology turns up in the late 1980; its basic meaning is to understand and transform the nature in the nanometer size range (10⁻¹⁰–10⁻⁷ m), directly through the manipulation and arrangements of atoms. Material properties mainly include chemical properties, physical properties (density, viscosity, particle size, specific heat, thermal conductivity, linear expansion coefficient, etc.), and mechanical properties (elasticity, hardness, strength, plasticity, ductility, toughness, impact resistance, fatigue limit, etc.). At nanoscale, the crystalline material with the smallest unit holds anisotropy, meaning that crystal has different properties in different directions.

With the development of material science and technology, the scale of material has expanded from macro to micro, also from micro to nanoscale. Nanostructure scientific research at nanometer scale includes elaboration, observation, characterization, and analysis. Modeling of miniaturized materials has a lot of new properties, where characteristic space and time scales correspond to typical simulation methods. Modeling and applications of structures at typical space/time scales [1] are shown in Figure 1.

Figure 1.

Modeling and applications of structures at typical space/time scales [1]. (a) Some characteristic space and time scales and (b) characteristic times of typical simulation methods.

From Figure 1, for electronic structure, DFT-QHA method is commonly used, while for atomistic structure, Monte Carlo method or molecular dynamics is used. However, real materials are never perfectly isotropic. At nanoscale, the crystal structure reorientations are not a direct result of the applied stress but are a geometrical requirement. Bulk anisotropy due to crystal orientation is therefore induced by plastic strain and is only indirectly related to stress. Recrystallization during annealing usually changes the crystallographic texture but does not cause randomness [2]. When the grain size is reduced to the nanoscale, its Young’s modulus and hardness have great changes. The microscopic deformation mechanism is very complex for nanograins, which depends on the shape of the crystal orientation, surface effect, the substrate effect, the grain boundaries affect, etc. In some cases (e.g., composite materials, single crystal), the differences in properties for different directions are so large that one cannot assume isotropic behavior, which leads to another behavior—anisotropic.

Anisotropic behavior is widespread among the variety of materials, and it has a great meaning in practical engineering application in material science field. Anisotropy is the property of being directionally dependent, as opposed to isotropy, which implies identical properties in all directions. It can be defined as a difference, when measured along different axes, in a material’s physical or mechanical properties. From the plastic performance perspective, the anisotropy of plastic deformation of the material shows the performance of texture, ears, and other aspects. Mechanical working produces preferred orientations or crystallographic textures. The primary cause of anisotropy of plastic properties is preferred orientation of grains. In contrast, anisotropy of fracture behavior is largely governed by mechanical alignment of inclusions, voids, and grain boundaries. However, from the elastic performance perspective of crystal structure, the mechanism of elastic anisotropy behavior is still not being recognized, especially for knowing the connection between atomic bonds and mechanical properties of anisotropic crystals structures at nanoscale. Therefore, a comprehensive understanding of the intrinsic nature and the anisotropy behavior at nanoscale has become an urgent task.

Essentially, material properties due to anisotropy are closely related to its internal structure and deformation conditions in the research field of materials science. The microstructure provides a link between processing (how a material is made) and properties (how a material behaves) [3]. It relates with the thermal, electrical, and mechanical properties of a crystal [3]. Alongside the experiment and theory, numerical simulation is a way additional access to the understanding of physical systems. Indeed, it can calculate experimentally measurable quantities and predict properties that are inaccessible in the laboratory or in the model to be validated. In general, both the chemical structure and the microstructure of a material control its properties, of which the chemical structure is relatively fixed and the microstructure depends strongly on how the material is made. However, from mesoscopic level to atomic scale, the drop of material scale makes the contained atoms of nanometer system greatly reduced; besides, macroscopic quasistationary continuous band disappeared, thus showing the energy level separation. The quantum size effect makes physical properties of nanosystem different from other conventional materials, leading to many novel features. At the nanometer scale, anisotropic materials exhibit differences due to these directional arrangements of atomic structure, the force between atoms, such as van der Waals force (caused by hydrogen bonding, ionic bond), typical Coulomb force (caused by charge), and other forces caused by covalent bond, atomic distortion, defects, hydrophilic, etc. Numerical simulation of anisotropic plastic behavior has been investigated by Yonggang [4] in Harvard University, while the numerical simulation and elastic properties of crystal structure at nanoscale have not been investigated systematically.

Although material performance is not the same in different directions in crystal, however, it has a strict symmetry. Anisotropic behavior of crystalline/structure can be reflected in the constitutive relation, where the elastic coefficient matrix corresponding to crystal system is certain. For different crystal/system, the number of cubic crystal has the least number of 3, hexagonal crystal 6, monoclinic crystal 13, and triclinic crystal with the most number of 21. Elasticity coefficient has important practical significance on the scientific basis and engineering applications of the material. Different crystal systems can be characterized exclusively by their symmetries. Each crystal has a certain level of crystal symmetry with corresponding different number of independent elastic coefficients. By obtaining elastic coefficient matrix, homogenized elastic properties (bulk, shear, and Young’s moduli) of polycrystals at larger scale can be determined, using the Voigt-Reuss-Hill estimation, for example. Because of complex issues mentioned above, for different anisotropic structure, selecting the appropriate modeling approach is particularly critical.

2. Nanostructural modeling and simulation methods at nanometer scale

With the development of the modeling and computer technology, the scale of material has expanded from macro to micro, also from micro to nanoscale. Nanomechanics aims to study fundamental mechanical properties of material structure at the nanoscale. Nanoscale modeling methods mainly include quantum molecular (QM), Monte Carlo (MC), molecular (structural) mechanics (MM) [5], molecular dynamics (MD) [6], nonlocal continuum theory [7], the ab initio calculation [8], tight-binding molecular dynamics (TBMD), and the density functional theory (DFT) [9]. Typical application, different space scale, the corresponding time scales, and simulation methods in computational materials science fields have been discussed by Raabe [1, 10]. Molecular simulation is also known as molecular modeling, which refers to theoretical methods and computational techniques to model or imitate the behavior of molecular. The model size at different scales has gradually decreased to continuum and mesoscopic level (material model diameter greater than 10⁻⁴ m), the mesolevel (about 10⁻⁶–10⁻⁴ m), the microlevel (about 10⁻⁷–10⁻⁶ m), and atomic/nanolevel (about 10⁻¹⁰–10⁻⁷ m). Figure 2 shows typical scales of material computational mechanics.

Figure 2.

Scales of material computational mechanics.

From Figure 2, scales in computational mechanics field can be divided into four scales: the macro (macroscale), mesoscopic (mesoscale), microscopic (microscale), and nano (nanoscale). Correspondingly, the continuum mechanics is commonly used to solve the macroscale structures, where finite element method can be used. For molecular dynamics for microscale, mesoscale, and nanoscale structures, for example, the discrete model can be an option at mesoscale. For quantum mechanics for both nanoscale and atomic scale structure, for example, density functional theory can be an effective option.

During the last decade, nanomechanics has emerged on the crossroads of classical mechanics, solid-state physics, statistical mechanics, materials science, and quantum chemistry. Numerical methods represent the most versatile computational method for the various engineering disciplines, and the scale of material modeling is gradually transited from bulk scale to nanoscale in Figure 3.

Figure 3.

Length/time scales using various methods in computational materials mechanics.

From Figure 3, the simulation approach for the structure in each domain is not all the same. The internal distance of each structure from quantum scale (nanoscale and atomic scale) to bulk scale covers the range of nanometer to micrometer and the time scale within the range of femtosecond to seconds.

Correspondingly, numerical modeling in computational material mechanics includes several analysis techniques such as ab initio methods, molecular dynamics, finite elements, boundary elements, distinct elements, and other numerical approaches that depend on the material/structure at different scales. Relationship between modeling methods and time/length scale is shown in Table 1.

Modeling methods	Length scales (m)	Time scales (s)
Quantum mechanics method	<10^−8	10^−15–10^−12
Molecular dynamics method	10^−8–10^−6	10^−12–10^−9
Micromechanics methods	10^−6–10^−4	10^−9–10^−3
Continuum method	>10^−3	>10^−3

Table 1.

Relationship between modeling methods and time/length scale.

Table 1 shows typical space, time scales, and simulation methods. These methods can separately link the relationship between modeling methods and real structure size. The fundamental characteristic of numerical methods is that a structure is discretized into small elements. Then the corresponding models that describe the individual elements and their interactions are constructed. The scale of material has expanded from macro to micro, also from micro to nanoscale.

2.1 Classical quantum mechanics and molecular dynamics methods

2.1.1 Brief introduction of quantum mechanics (QM) method

In the early 1900, the basics of quantum mechanics were gradually proposed by Schrödinger and others to study the quantum phenomena when the dimension of the system is in the size of atoms and molecules. Schrödinger equation [11] is proposed and it plays the same role as Newton’s equation in classical systems but exactly describes the wave function of a system that evolves with time. Then, some approximations [12], such as the first ab initio Hartree-Fock (HF) calculations on diatomic molecules, were carried out at the Massachusetts Institute of Technology in 1956 [13]; thus, the numerical resolution of the time-independent Schrödinger equation was calculated. These traditional methods are briefly introduced as follows.

First-principles calculations are also called ab initio method, which refers five fundamental constants of physics—mass of the electron, the electron charge, Planck’s constant, the speed of light, and Boltzmann constant. The calculation is to compute the system integrating of all the electronics solving the Schrodinger equation without considering any empirical parameters and thus can reasonably predict the status and nature of microscopic systems. Besides, it is a strong theoretical approach developed from quantum mechanics and quantum chemistry, mainly to study hundreds of atoms or periodic crystal material. First-principles calculations can be used not only to explain the physical nature of the experimental phenomena, relationship between microscopic electronic structure and the physical parameters, but also to design new materials and to theoretically predict many physical properties.

Density functional theory is a kind of first-principles molecular dynamics calculations [14], which is also called ab initio molecular dynamics. First-principles molecular dynamics make up shortcomings that the classical molecular dynamics need to know the atomic interaction potential. The atomic interaction potential is solved by DFT [15] and then substituted into the atomic motion equations to solve the trajectory of the particles of the system. The main difference between first-principles molecular dynamics and classical molecular dynamics is to calculate the atomic interaction force [14, 15]. Advantages of first-principles molecular dynamics are as follows: (1) quantum mechanics is used to describe the interaction between electrons and atomic structure and (2) it is able to calculate the effect of nonharmonic crystal thermal vibrations. The downside of DFT is that it is time-consuming [16].

2.1.2 Brief introduction of molecular dynamics (MD) method

The molecular dynamics (MD) simulation to describe the atomic behavior by classical models and equations was reported in 1956 [17]. Molecular dynamics is to calculate a set of molecular orbital phase space. The MD method uses a microclassic calculation, assuming that atomic motion can be described by Newton’s motion equations, which means that atomic motion is associated with a particular track. This assumption is feasible and available only when the movement of quantum effects can be ignored and the adiabatic approximation is strictly satisfied. Molecular dynamics is a deterministic method that offers the possibility of a microscopic description of a physical system in consideration of all the interactions involved. The main advantage of this method is that it gives the information on the evolution of the system over time by numerically solving Hamilton equations of motion, Lagrange, or Newton.

Before showing how organized a simulation by molecular dynamics, it is important to define some necessary terms (NVE, NVT, NPH, and NPT ensembles) in such a simulation. The MD method is a deterministic simulation technique for evolving systems to equilibrium by solving Newton’s laws numerically. Molecular dynamics method can be used to study the molecule substances with lots of atoms (even up to millions atoms). The implemented potential was tested to reproduce the basic physical parameters of bulk, defects, and surface. The bulk physical parameters, such as the cohesive energy, lattice constant, and bulk modulus, were obtained on the supercell structure with a certain boundary condition, which were equilibrated by MD method with the set of a time step.

2.2 The AFEM methodology based on molecular mechanics method

Molecular mechanics method aims to study fundamental mechanical properties at the atomic scale. It has emerged on the crossroads of classical mechanics, solid-state physics, statistical mechanics, materials science, and quantum chemistry. Length/time scales using various methods in computational materials mechanics are shown in Figure 5.

From Figure 4, we can see that there is a gap between time record and space record. The QM method, such as the density functional theory (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. Molecular dynamics (MD) has obvious advantages in simulating macromolecular structure, with very high computational efficiency but for much dependency of atomic potential parameters. Besides, MD method is limited into a range of structure at atomic scale level, and it cannot effectively solve continuum mechanics problems like practical engineering application using finite element method (FEM). In fact, the physicist Garrett Lisi in 2007 has proposed “a theory of everything” and has proposed an E₈ model [18]. Garrett Lisi believes that all known physics particles and elementary various forces can be accommodated at the E₈ mode where the plurality of interactions between particles is naturally appeared.

Figure 4.

Length/time scales using various methods in computational materials mechanics.

So how to establish a connection between quantum-mechanical (QM), molecular dynamics (MD), and continuum mechanics, so as to predict or calculate macroscopic properties at the atomic scale based on some kinds of continuum methods?

These possible modeling methods cause great interest and concern of the scholars. Quantum-mechanical (QM) or ab initio methods mainly aim to solve electronic Schrodinger equation for atoms and molecules, as well as various approximations of that solution. Electronic structure of each atom is accounted for. Major methods include density functional theory, Hartree-Fock, and tight binding. QM methods are very accurate but extremely expensive. Otherwise, without consideration of the complex electronic structure of each atom, the whole atom calculation methods can be realized by modeling the structure/system just as a type of ball (soft sphere), using a series of empirical interatomic potentials. These major methods (molecular mechanics, molecular dynamics, Monte Carlo) are less accurate but relatively inexpensive. Although MD method can simulate the microscopic deformation in a certain degree, some effects due to size are too small but still cannot be obtained because of computational limitations of molecular dynamics simulations, and continuum mechanics finite element method cannot effectively simulate the microstructure deformation behavior. Because of these deficiencies of molecular dynamics and continuum methods, multiscale method came into being, which can be used to simulate in larger sizes without ignoring the deformation at smaller scale.

In addition to the traditional simulation methods that we briefly introduced above, a so-called atomic finite element method (AFEM) [19] at the atomic scale (or nanoscale) has attracted our great interest. In fact, models of crystal structures at the nanometer scale often tend to exhibit anisotropy. Modeling of anisotropic crystal structure needs to include the atomic and force field information, which leads to anisotropy due to the consequence of differences in the directional arrangements of atomic structure. The force between atoms, such as van der Waals force, typical Coulomb force, and others, should be considered.

In 2011, an atomic finite element (AFEM) as a bridge between molecular dynamics and continuum mechanics is proposed, which makes the classical finite element method (FEM) to model at nanoscale.

Considering Coulomb interaction, the total steric potential energy U by Rappe is as [20]:

$$U=\sum U_r+\sum U_{\theta }+\sum U_{\varphi }+\sum U_{\omega }+\sum U_{\mathit{vdw}}+\sum U_{\mathit{\text{coul}}}\hspace{1em}\mathrm{(1)}$$

where U_r is bond stretching, U_θ is bond angle bending, U_ϕ is dihedral angle torsion, U_ω is improper torsion, U_vdw is nonbonded van der Waals interaction, and U_coul is Coulomb (electrostatic) interaction.

At this scale, the interatomic forces are considered binding van der Waals, electrostatic force, and the covalent chemical bond. However, its application is limited as well, including the simple binary alloy and covalent bonds [19]. The application of AFEM in the field of ionic crystal is blank. Description of main force field of molecular mechanics is shown in Figure 5.

Figure 5.

Description of main force field of molecular mechanics.

From Figure 5, covalent bonds can be characterized to be equivalent beam by the molecular structural mechanics, and nonlinear spring element can be used to characterize the coupled field of van der Waals force and Coulomb force between atoms or charged ions by using Lennard-Jones potential. For AFEM modeling, the equivalent beam’s definition of other covalent bond can be approximated solved by a linear bond energy ratio according to C▬C bond energy. This book tries to approach this goal, elastic constants under various anisotropic crystal directions and homogenized elastic moduli of polycrystal become the research priorities of subject. In addition, the conventional DFT methods and MD methods are also used to compare mechanical properties and the simulation results of research objects.

As a structural mechanics approach [5] has successfully modeled the frame-like structure of C▬C covalent bonds by equivalent beam elements and an AFEM methodology has solved the interatomic bonds modeled by nonlinear spring elements, modeling of nanostructure is realized by an equivalent beam, and nonlinear spring elements become possible. Such modeling tool may be used to further widen the application fields of nanoindentation simulation with consideration of atomic forces at nanoscale level.

3. Research background on the multiscale modeling of concrete platform

Concrete is the world’s most widely used man-made material with the history of more than 160 years. It is a complex mixture of numerous compounds, which has been used in the application in the antiseismic building, civil engineering, construction, marine engineering projects, etc. [21] Our research background is based on some typical crystal structure of cement pastes in the reinforcing concrete, which accounts for the platform MuMoCC (Multi-Scale Modeling of Computational Concrete).

3.1 Modeling scales of concrete and cement hydration products

Concrete is a highly heterogeneous composite construction material whose microstructure contains randomly dispersed features. Concrete is primarily composed of cement paste, sand, coarse aggregate, and admixtures. Each of these ingredients is heterogeneous themselves at a certain level and has different stiffness and strength. Therefore, the heterogeneity of concrete exists in a variety of length scales.

3.1.1 Different scales of the concrete modeling

All the types of cement pastes, in general, are truly complex materials, porous, multicomponent, and multiscalar. Moreover, as the hydration products have a very dissimilar structure and sizes, the cement paste itself has different hierarchical levels of organization at different scales. Figure 6 shows schematic representation of top-down five scale levels model for normal concrete [22, 23], which gives level classification and typical structures of cement-based composite material.

Figure 6.

Level classification and microstructures of cement-based composite material [22, 23].

From Figure 6, the typical structures can be divided into five levels, from the scale of pavement (10⁻¹ m) down to the C▬S▬H solid phase (10⁻¹⁰ m) [22]. Each length scale is a random composite, and the description of the cement paste structure will be presented for each scale. A multiscale level model has been defined based on morphology and length scale of the structural elements [24]. The most usual classification divides the multiscale structure into four size ranges: macroscale, mesoscale, microscale, and nanoscale.

3.1.2 Cement hydration process and typical hydrated structures

Hydration refers to the chemical reactions between cement and water; the hydration process causes cement paste to first set and then harden. Meanwhile, hydration products—some new solid phases—are formed by hydration. Cement hydration refers to a lot of complex mixtures of numerous compounds. The major compounds found in cement are tricalcium silicate (C₃S), dicalcium silicate (C₂S), tricalcium aluminate (C₃A), tetracalcium aluminoferrite (C₄AF), and gypsum. Some tricalcium silicate produces C▬S▬H gel and calcium hydroxide (CH) is formed at the early stage of hydration, while dicalcium silicate produces C▬S▬H gel and CH is formed at a later stage in the hydration process.

Figure 7 summarizes the order of the different hydrates and their relative quantities [25].

Figure 7.

Crystals of silicates during solidification process [25].

Figure 7 shows the degree of hydration, particularly at later times. Dashed curve of the porosity decreases in proportion with the rate of hydration of the cement paste, due to the progressive filling empty initial intergranular by various hydrates. The C▬S▬H gel is the most important product of hydration. Furthermore, it is found that after a few hours, ettringite increases, reaches a maximum, and then decreases in favor of monosulfoaluminate, calcium aluminates, and hydrated aluminates C₄(A,F)H₁₃. At the end of stage (after 28 days), the initial cement has hydrated, and the paste has undergone final set.

Typical microstructures and scanning electron microscope (SEM) of cement hydration products by Regourd [26] in 1975 is shown in Figure 8.

Figure 8.

Typical microstructures of cement hydration products [26].

In Figure 8, we can see the SEM image of C▬S▬H microstructures that are nanocrystalline compounds of nanometric particle aggregated with each other. Among the SEM of a clinker grain, the well-defined side and the pseudohexagonal contour of the upper face of the crystal of C₃S are noted. Aluminates are as a “stick” structure between the crystals C₃S calcium silicates and/or C₂S. Alkali sulfates are condensed on the surfaces of the silicate crystals during solidification process. The hydration of C₂S and C₃S forms calcium hydrosilicate (C▬S▬H), which is accompanied with calcium hydroxide (CH). The continuous and relatively rapid deposition of hydration products (primarily C▬S▬H gel and CH) into the capillary porosity occurs, which is the space originally occupied by the mix water, leading to a large decrease in the total pore volume and a concurrent increase in strength [26].

C▬S▬H is responsible for most of the cement paste properties, which is amorphous and heterogeneous in composition. It is well known that the C▬S▬H gel can be classified as inner product and outer product, of which the inner product of C▬S▬H gel grows inward the boundaries of the clinker grains, occupying the place of the anhydrous phases, while the outer product of C▬S▬H gel is formed outward the boundaries of the clinker grains, in the water filled space [27]. In general, the inner product has a more compact and more amorphous structure, while the outer product has been found to form bundles of fibers radiating from the cement grains [28]. The composition of both products seems to be similar, although some authors report slight differences in the Ca/Si ratio, which is higher for the inner product [29].

3.2 Background on reinforcing concrete and MuMoCC platform

3.2.1 Research background on concrete reinforced by nanotubes

Nowadays, the use of self-compacting concrete (SCC) is progressively changing the method of concrete placement on building sites [30]. Concrete is a heterogeneous material composed of a mixture of sand, aggregates embedded in a hardened cement paste. Cement is the binder component of concrete, the glue that holds the filler together to create a uniform, strong material. The cement notation uses letters that symbolize the main components expressed as oxides: CaO (C), SiO₂ (S), Al₂O₃ (A), Fe₂O₃ (F), H₂O (H), SO₃ (S), etc., forming crystal morphology. As a reinforcing material in cement matrix, carbon nanotubes have high stability and significantly improved its mechanical properties, impact resistance, antistatic properties, and wear resistance properties. Therefore, the development of a cement matrix of carbon nanotubes becomes scholars’ interest. In 2012, the conception that “cement nanotubes as a natural means for reinforcing concrete [31]” is proposed, shown in Figure 9.

Figure 9.

Main cement paste types and nanotube reinforcing concrete expectation [31].

From Figure 9, it has been concluded that nanotubes make calcium silicate hydrates an ideal cement-reinforcement paste, of which the stress deformation of portlandite nanotubes increases almost linearly up to a strain of 27% [31]. Mechanical properties of cement matrix reinforced by carbon nanotubes are super excellent; it is not only resistance to fatigue and creep but also dimensionally stable. Carbon nanotubes can be formed by rolling a cylindrical graphene and thus can be used as a reinforcing material of brittle materials such as cement.

3.2.2 Research background about MuMoCC platform

Modeling of nanostructure and atomic/ionic crystal and their elastic properties has become a hot field for many scholars [32, 33, 34, 35, 36, 37, 38, 39, 40, 41]. Another research background is that a hierarchical multiscale modeling of the behavior of cement-based materials and the multiscale modeling of computational concrete (MuMoCC) platform [42, 43] have been proposed; thus, behavior of hydrates (like C▬S▬H, CH, etc.) and admixtures (like calcite, CNT, etc.) are needed to be investigated in detail. Elastic properties (Young’s modulus, Poisson ratio, etc.), tensile strength, and postfailure tensile stress displacement are thus the unique data that are transferred from cement paste to mortar [43]. Besides, the elastic displacements are found in every pixel, and the average strain and stress is computed and averaged over the entire microstructure to give the effective elastic properties. Input data of these solid phases at the nanoscale are needed. Different scales in structural modeling of concrete using HYMOSTRUC model of multiscale structures by UT Delft [44] and MuMoCC platform [42, 43] are in Figure 10.

Figure 10.

Different scales in structural modeling of concrete. (a) HYMOSTRUC model of multiscale structures by UT Delft [44] and (b) the representative MuMoCC platform [42, 43].

The methodology has been successfully applied to determine mechanical as well as transport properties of cement-based materials, and also the evolution of these properties according to time, and their modification due to degradation phenomena such as leaching.

The prediction of the mechanical effective properties has been an active research area for the few decades. Numerical models seem to be a well-suited approach to describe the behavior of semianalytical material models, because there is no restriction on the geometry, on the material properties, on the number of phases in the composite, and on the size. Except for the experimental studies, either micro- or macromechanical methods are used to obtain the overall properties of material phases. In the macromechanical approach, the heterogeneous structure of the composite is replaced by a homogeneous medium with anisotropic properties. Micromechanical method provides overall behavior of the matrix-inclusion composites from known properties of their constituents (inclusion and matrix) through an analysis of representative volume element or a unit cell model. The advantage is not only the global properties of the composites but also various mechanisms such as damage initiation and propagation, through the analysis.

4. Objectives and characteristics of typical anisotropic crystals at nanoscale

Under the research background of the conception that “cement nanotubes as a natural means for reinforcing concrete” and the platform construction of MuMoCC (multiscale modeling of computational concrete), 10 kinds of crystals structures are, respectively, studied so as to achieve the calculated values of their mechanical properties (Young’s modulus in particular) and to perfect the crystal elastic theory of typical crystal structures simultaneously.

4.1 Objectives of typical anisotropic crystal structures at nanoscale

4.1.1 Nanoscale modeling on concrete reinforced by nanotubes

In 2016, Eftekhari and Mohammadi [45] have investigated the molecular structure of CNT-reinforced C▬S▬H phase, the C▬S▬H model is 3.925 nm × 3.626 nm × 4.768 nm, shown in Figure 11.

Figure 11.

The molecular structure of CNT-reinforced C▬S▬H phase [45]. (Atom types: red, oxygen; cyan, hydrogen of water molecules; green, calcium; yellow, silicon; blue, CNT) (diameter of 0.813 nm with length of 4.0 nm).

From Figure 11, the enhanced mechanical properties of C▬S▬H reinforced by embedding CNT in its molecular structure are simulated by MD method. It is indicated that the tensile strength of CNT-reinforced C▬S▬H is substantially enhanced along the direction of CNT as compared to the pure C▬S▬H. In addition, CNTs can severely intensify the “transversely isotropic” response of the CNT-reinforced C▬S▬H, which can be used for multiscale simulation of crack bridging at macroscale specimen of CNT-reinforced cement.

In all, from a cementitious composite [46] reinforced with carbon CNT and/or polypropylene microfibers, a dense C▬S▬H formation that appears to be tightly bonded to the SWCNT, producing reinforcing behavior [47]. In fact, nanotubes and graphene as their superior mechanical properties have been more and more confirmed to be the reinforcing materials to enhance the toughness of brittle materials [48].

4.1.2 Objectives of typical anisotropic crystal structures at nanoscale

The use of CNT as a reinforcing material is intended to move the reinforcing behavior from the macroscopic to the nanoscopic level, which has the advantage of extremely high strengths, Young’s modulus, elastic behavior, and advantageous thermal properties. Crystal elastoplastic behavior has been the long concern for scholars, with the introduction of nanoscale conception; the performances of anisotropic crystal or polycrystalline become very interesting. The investigation of the crystal elastic simulation has become an urgent issue to be resolved. Combining the introduction mentioned above and the actual research needs, Young’s modulus of typical crystal structures is especially concerned, and models of these considered structures are shown in Figure 12.

Figure 12.

Crystal modeling objects and corresponding real structures at nanoscale. (a) CNT and hexagonal graphene. (b) Cubic CaO structure. (c) Cubic MgO structure. (d) Hexagonal CH structure. (e) Hexagonal calcite structure. (f) Monoclinic 11 Å tobermorite. (g) Monoclinic 14 Å tobermorite. (h) Monoclinic gypsum structure. (i) Amorphous disordered C▬S▬H structure.

From Figure 12, depending on crystal types of typical structures and practical research needs, these crystals structures are divided into cubic crystals (CaO and MgO), hexagonal crystals (calcite and portlandite), monoclinic crystals (11 Å tobermorite and gypsum), monolithic disordered structures (C▬S▬H (I) with Ca/Si of 0.67 and C▬S▬H (II) with Ca/Si of 1.67), and other crystals structures (CNT, graphene).

This book mainly includes three types of the cubic, hexagonal, and monoclinic crystals, as well as the other crystal structures considered. Abandoning the lengthy discussion of the important description in these typical structures, we mainly focus on modeling and calculating the mechanical properties of these typical anisotropic crystal structures at nanoscale with a high efficiency and then determining homogenized elastic moduli of these structures with a larger scale. The numerical simulation to determine elastic modulus of these typical structures is a meaningful research.

4.2 Determination of crystal-independent elastic constant coefficients

In matrix format, the stress-strain relation shows the 36 (6 × 6) independent components of stiffness. Crystal elastic anisotropy can be expressed by a fourth-order tensor matrix of strain-stress relation; Hooke’s law thus can be written as a general symmetric 6 × 6 matrix in Voigt notation.

$$\left[\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\sigma }_4\\ {\sigma }_5\\ {\sigma }_6\end{array}\right]=\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& c_{14}& c_{15}& c_{16}\\ c_{21}& c_{22}& c_{23}& c_{24}& c_{25}& c_{26}\\ c_{31}& c_{32}& c_{33}& c_{34}& c_{35}& c_{36}\\ c_{41}& c_{42}& c_{43}& c_{44}& c_{45}& c_{46}\\ c_{51}& c_{52}& c_{53}& c_{54}& c_{55}& c_{56}\\ c_{61}& c_{62}& c_{63}& c_{64}& c_{65}& c_{66}\end{array}\right]\left[\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\\ {\epsilon }_3\\ {\epsilon }_4\\ {\epsilon }_5\\ {\epsilon }_6\end{array}\right]\hspace{1em}\mathrm{(2)}$$

Elastic constant matrix C_ijkl refers to the parameters that resist physical resilience when the material is subjected to the stress. Elastic constants reflect the macroscopic mechanical properties. Further reductions in the number of independent constants are possible by employing other symmetry considerations.

Elastic constants of single crystal can be determined by various methods, such as classical MD method calculation [49], the universal linear-independent coupling strains [50], first-principles methods [51, 52], and experimental measurements. For an aggregate crystalline structure, Voigt [53] proposed to find the average of stresses at all possible lattice orientations under a given strain, while Reuss [54] suggested to calculate the average of strains under a given stress instead. In addition to the experimental technique, there exist mainly two methods of calculating the elastic constants: (1) stress-strain curve relationship and (2) strain energy-strain relationship [55]. Elastic constants can be determined by the first-order derivative of the curve. Finally, the bulk and shear modulus of the crystalline phases can be estimated based on averaging methods.

The independent elastic coefficients corresponding to each crystal are listed in Table 2.

		Crystal type	Number of Cij	Independent elastic coefficient
		Isotropic	2	c11, c12
a = b = c	α = β = γ = 90°	Cubic	3	c11, c12, c44
a = b ≠ c	α = β = 90°, γ = 120°	Hexagonal	5	c11, c12, c13, c33, c44
a = b ≠ c	α = β = γ = 90°	Tetragonal	6	c11, c12, c13, c33, c44, c66
a ≠ b ≠ c	α = β = 90°, γ = 120	Trigonal	6	c11, c12, c13, c33, c14, c44
a ≠ b ≠ c	α = β = γ = 90°	Orthorhombic	9	c11, c22, c33, c12, c13, c23, c44, c55, c66
a ≠ b ≠ c	α = γ = 90°, β ≠ 90°	Monoclinic	13	c11, c22, c33, cl2, c13, c23, c44, c55, c66, cl6, c26, c36, c45
a ≠ b ≠ c	α ≠ γ ≠ β ≠ 90°	Triclinic	21	c11, c12, c13, c14, c15, c16, c22, c23, c24, c25, c26, c33, c34, c35, c36, c44, c45, c46, c55, c56, c66

Table 2.

Independent elastic constants for seven kinds of crystal symmetry.

As we can see in Table 2, the elastic properties of an isotropic material can be fully determined by calculating independent coefficients; thus, bulk modulus and shear modulus can be finally determined, using the relationship between the homogenized elastic moduli and elastic constants, and Young’s modulus, shear modulus, the speed of sound, and other physical properties of the material.

Simulations of typical structures to determine their elastic moduli based on AFEM, DFT, and MD are shown in Figure 13.

Figure 13.

Elastic moduli determination based on AFEM, DFT, and MD.

Independent elastic constants for typical crystals/structure considered are in Table 3. A kind of AFEM method in Chapter 2 to simulate the previously proposed models at the atomic scale is shown in Chapter 2.

Crystal type	Lattice relations	Cij number and elastic coefficient		Method	Object
Cubic	a = b = c; α = β = γ = 90°	3	c11, c12, c44	DFT	CaO, MgO
Hexagonal	a = b ≠ c; α = β = 90°, γ = 120°	5	c11, c12, c13, c33, c44	DFT	Calcite, portlandite
				AFEM	Portlandite, graphene
Monoclinic	a ≠ b ≠ c; α = γ = 90°, β ≠ 90°	13	c11, c22, c33, c12, c13, c23, c44, c55, c66, cl5, c25, c35, c46	DFT	Tobermorite, gypsum
				MD	Tobermorite
Others	—	—	—	AFEM	CNT
				MD	C▬S▬H with C/S = 1.67

Table 3.

Independent elastic constants for typical crystals/structure considered.

From Table 3, we can see that during the investigation of nanoscale modeling and mechanical properties, the DFT, AFEM, and MD method are especially concerned according to the modeling requirement of different crystal structures investigated. In summary, our research will provide further detailed parameters of our research objects and thus leads to a certain practical significance.

5. Mechanical properties calculation of typical structures for elastic properties

To get all 21 elements, only one of the six independent strains (ε₁, ε₂, ε₃, ε₄, ε₅, ε₆) nonzero each time. Using the average stress tensor,<σ_ij>, one can then define the effective elastic tensor by <σ_ij>=<C_ijkl>ε_ij, where ε_ij is the applied strain. Using this technique, all 21 elements can be evaluated. Methodology used for elastic constants evaluation is in Table 4.

Method	The method of measurement or calculation	Measured or calculated parameters
Experiment	Static, quasistatic method	Strain, stress
	Indentation method	Strain, stress
	Pulse-echo method	Velocity
	Resonance ultrasonic spectrum (RUS) method	Frequency
	X-ray or neutron method	Diffraction or scattering spectra
	Brillouin scattering method	Diffraction or scattering spectra
Calculation	Pending elasticity parameter method	Parameters from experiment
	Classical molecular dynamics	Energy, stress, strain
	Monte Carlo method	Energy, stress, strain
	First-principles molecular dynamics	Energy, stress, strain
	DFT-QHA method	Phonon density of states, Helmholtz free energy
(Present)	Atomic finite element method	Strain, stress

Table 4.

Methodology used for elastic constant evaluation.

In Table 4, the experimental technique is to build a microstructure using one of the above rules, assign elastic properties to the phases, and use the programs to compute the effective elastic moduli of the model.

Methodologies for elastic constant evaluation of considered structures are in Table 5.

Method	Measurement or calculation methods	Measured or calculated parameters
Experiment	Indentation method	Strain, stress
Classical numerical calculation	Classical molecular dynamics	Energy, stress, strain
	First-principles molecular dynamics	Energy, stress, strain
Method developed	AFEM	Stress, strain

Table 5.

Methodologies for the elastic constant evaluation of considered crystal structures.

From Table 5, it shows advantages and disadvantages of each modeling method. Besides, different typical structures are modeled using different modeling method to calculate mechanical properties. Researchers in materials science may be tasked with a number of challenging goals such as the development of new compounds or new molecular structure that may be faced with simply understanding and describing fundamental processes, explaining why one particular material is better than another. Modeling can address all of these challenges, provided that the method is fast, accurate, and works at the atomic scale.

5.1 AFEM method for C▬C and hexagonal structures by ABAQUS

5.1.1 AFEM methodology for C▬C and hexagonal structures

The AFEM/FEM linkage provides a seamless multiscale computation method for static analysis [56]. The AFEM element can be implemented into the ABAQUS finite element program with its USER-ELEMENT subroutine to solve several atomic-scale problems. To ensure that the AFEM/FEM multiscale computation method accurately represents the material behavior at both atomic and continuum scales, the continuum FEM elements should be based on the same interatomic potential as AFEM elements. The nonlinear spring element’s definition of ionic bond or metallic bond can be calculated by MATLAB according to proper potential functions in various corresponding equilibrium conditions and then be imported into ABAQUS. For AFEM modeling, the calculation of U_r, U_θ, U_φ, and U_ω of Eq. (1) has been investigated in the reference [5], where the molecule's chemical energy is equivalent to the beam strain energy. The coupled field of van der Waals force between atoms or charged ions may be considered by means of the Lennard-Jones potential. Nonlinear spring elements can represent U_vdw potential.

5.1.2 The atomic finite element method modeling and ABAQUS

As a so-called atomic finite element method (AFEM) as a bridge between molecular dynamics and continuum mechanics has been proposed based on molecular mechanics method, where main force fields can be divided into two parts—covalent bond interaction and nonbonded force field. The nonlinear ABAQUS is a commercial finite element analysis program, which has been in use to analyze the structure in mechanical/material industry, civil engineering, and other field industries. The software is capable to analyze the stress and strain buildup in a variety of problems, to design the AFEM multiunit/element models. For defining nonlinear spring behavior, we can define nonlinear spring behavior by giving pairs of force-relative displacement values, which should be given in ascending order of relative displacement and should be provided over a sufficiently wide range of relative displacement values so that the behavior is defined correctly. Nonlinear spring is used to define force-relative displacement relations [57]. For crystal structure before homogenization, its elastic moduli in three directions are different, so is there a link of theoretical expression formula between axial elastic moduli and elastic coefficient? From Hooke’s law, it can be inferred indirectly that axial elastic moduli are closely related to elastic coefficient. For mechanical properties, the anisotropy of the material shows the performance difference of axial moduli in three axial directions; formula between axial elastic moduli and elastic coefficient is derived in Chapter 2.

5.2 DFT method for typical structures by CASTEP and homogenization

Elastic constants can be calculated by CASTEP module using DFT method. For DFT method, the Hamiltonian and Kohn-Sham equation [58] as an approximation to simplification of Schrodinger equation can be used to solve the multiparticle systems. However, at a sufficiently fine scale, all materials are heterogeneous or can be said to be homogeneous, which introduces an important concept—homogenization.

5.2.1 The density functional theory calculation and CASTEP

The Hamiltonian and Kohn-Sham equations as an approximation are used for the simplification of Schrodinger equation to solve the multiparticle systems [59]. DFT includes not only the nucleus and electron kinetic energy term, but also interaction term of nucleus-nucleus and electron-electron [49]. CASTEP is an ab initio quantum mechanical program employing density functional theory (DFT) to simulate the properties of a wide range of solid materials at the atomic scale. CASTEP uses quantum mechanical calculations to study problems in chemicals and materials research, which is able to predict the structure of a material as well as many essential properties. In particular, it can predict electronic properties (e.g., band gaps and Schottky barriers), optical properties (e.g., phonon dispersion curves, polarizability, and dielectric constants), or physical properties such as elastic constants.

The computed elastic tensor of engineering material will have the symmetry of elastic constants matrix C_ijkl in general. One can analytically average this tensor, using the full C_ijkl form of the elastic moduli tensor, over the three Euler angles (α, β, γ) [60]. For solid material at macroscale, it can be seen that polycrystalline material is a collection of small grain, so elastic, plastic, and physical properties of the polycrystalline at macroscale are commonly thought to be relative with the mechanical properties of single crystals and polycrystalline themselves, including grain orientation and grain boundary structure.

First-principles calculations allow researchers to investigate the nature and origin of the electronic, optical, and structural properties of a system without the need for any experimental input. CASTEP is thus well suited to research problems in solid-state physics, materials science, chemistry, and chemical engineering where empirical models are lacking and experimental data may be sparse. In these areas, researchers can employ computer simulations to perform virtual experiments, leading to tremendous savings in costly experiments and shorter developmental cycles. Key procedure includes a transition state search algorithm that greatly facilitates determination of reaction profiles and energy barriers, essential to an understanding of kinetics. The full 6 × 6 tensor of the elastic constants can be predicted for a periodic structure of any symmetry. In the density functional theory part of the book, CASTEP is used to calculate the elastic constants of the single crystal structure, seen in Chapter 3.

5.2.2 Homogenized moduli by Y-parameter method based on elastic constants

Although single crystals generally exhibit anisotropic mechanical behaviors, macroscopic averaging assuming a randomly distribution of initial crystals leads to isotropic mechanical properties such as Young’s modulus. Polycrystalline material is composed of single crystal, the probability distribution depicting upward almost the same as in three-dimensional space, resulting in an isotropic behavior of mechanical properties at macro-scopic scale. Meanwhile, Voigt model based on the assumption that all grains have the same deformation strain is given, by which the upper bound of polycrystalline and elastic constitutive relation can be given, while the Reuss model is based on the hypothesis of the same stress for polycrystalline, where the lower bound of polycrystalline and another form of elastic constitutive relation can be given. However, the constitutive relation by Voigt and Reuss model is not always accurate, especially for the lower bound. Y-parameter of cubic structure was firstly introduced and presented by Zheng and Min [61, 62] in 2009. Based on the elastic constants of cubic crystal, the corresponding Y-parameter, as new evaluation criteria, has been proposed. Elasticity of single crystal and mechanical properties of polycrystalline material has been closely integrated. By comparing various calculations method to determine homogenized moduli of the polycrystalline material composed of a single crystal, for example, the certain stress of Reuss model [54] and the certain strain of Voigt model [53], the Y-parameter, in the theoretical calculation to forecast the elastic modulus of polycrystalline material, has a high consistency. So our work aims to complete this part of the theoretical analysis for polycrystalline, with comparison of elastic constants measured by references.

This interesting Y-parameter for cubic crystals can be used to determine the homogenized moduli of cubic structures. Moreover, the structural and the elastic properties of CH and calcite are investigated by ab initio plane-wave pseudopotential density functional theory method in Chapter 3, and then the credibility of Y parameters for determining elastic moduli of CH and calcite structures is proved. We will then describe these interesting Y parameters of cubic and hexagonal polycrystalline in detail in Chapter 4. In fact, due to the difficulty to accurately measure elastic constants of each crystal for polycrystalline, the fourth-order tensor connections between mean stress and mean strain called the elastic constitutive relation are impossible to be determined by numerical integration method. Then elastic moduli (Young’s, shear, and bulk moduli) can be obtained by means of Voigt and Reuss bounds, which are based on the C_ij calculated. As the complication of the monoclinic crystals (11 Å tobermorite and gypsum), elastic moduli of monoclinic crystals will be still calculated by the classical Voigt-Reuss-Hill estimation, in Chapter 4.

5.3 MD method and indentation simulation for C▬S▬H by LAMMPS

5.3.1 Elastic moduli by stress-strain curve slope method using MD

To simulate large-scale structures instead of these small crystals by DFT, the corresponding potentials were implemented in large-scale atomic/molecular massively parallel simulator (LAMMPS) code [63]. Note that LAMMPS is an open-source code and thus can be actively modified by users as well as the developers. At the end of each time interval, a new box length is computed. The elastic modulus of typical structures can be approximated and computed from the slope of the linear portion of the stress-strain curve.

In the molecular dynamic part of this book, LAMMPS (large-scale atomic molecular massively parallel simulator) is a classical MD code used in this study. LAMMPS is developed by Steve Plimpton and his group of Sandia National Laboratories, which can be used to calculate the tension, compression, and shear deformation processes of the large-scale supercell structures. An input file was generated for LAMMPS code to calculate the response of C▬S▬H structures to uniaxial loading, which leads to strain-stress data to calculate elastic moduli we aim to obtain in Chapter 5.

5.3.2 Elastic moduli by nanoindentation simulation

Mechanical properties of concrete as a composite material depend on microstructure and even depend on properties of crystal of each paste at nanoscales. Concrete is a complex heterogeneous material whose mechanical properties can vary substantially from point to point. Nanomechanics and nanotechnology are used to explore the composition, behavior, and nanomechanical properties of cement paste phases with submicron resolution. Nanoindentation technique bridges the gap between atomic force microscopy (AFM) and macroscale mechanical testing. The corresponding load-depth curve can be measured directly by nanoindentation simulation; thus, elastic modulus can be determined.

Therefore, nanoindentation technique as a useful tool can be used to better understand the submicron mechanical properties of cement pastes. Because of its small probe size, indentation tools are previously mainly used to measure local material properties in small, thin, and heterogeneous materials. Nanoindentation simulation becomes a hot spot in recent years; based on continuum mechanics, nanoindentation simulation can reveal atomic deformation mechanism under a indent force; the smallest unit with embedded atomic force is used in the indentation area of the multiscale model. Some methods or programs can be used partially to design and analyze multiscale model, which includes the finite element method, the boundary element method, finite difference method, and the discrete element method. ABAQUS can be used to program and simulate the nanoindentation process even with the nanoscale unit. The nanoindentation model with the unit of a nanometer size will be simulated and discussed in Chapter 6.

6. Conclusions

The development of modeling tools to describe and predict the mechanical properties of structural materials proves an undeniable practical importance. It is now well recognized that such an objective can be achieved through the linking of the structure of materials at the microscopic nanoscale or with their performance. Nanoscale modeling and mechanical properties by using the density functional theory, a so-called atomic finite element method, and the classical molecular dynamics methods are especially concerned in this work. Such tools will be more accurate and they can relate the structure of materials at the atomic or microscopic scale with their performance. In quantum mechanics field, when applied to real systems including several atoms and electrons, a proper solution of the stationary equation remains a challenge due to its complexity. However, the mean field theory (electrons move as independent particles in an effective potential generated by the ions and the other electrons) is also effective. Especially for molecular dynamics simulations, the model scale can be greatly extended.

In addition to the tradition of DFT methods and MD, a methodology for AFEM modeling at the atomic scale (or nanoscale) is explored since the macroscopic properties depend largely on the physicochemical properties of the interatomic bonds. This book is committed to model and to calculate mechanical properties (especially Young’s modulus) of the different typical anisotropic crystal structures by programming using three methods above. It therefore connects nanoscale modeling and continuous pattern of deformation behavior by identifying relevant parameters from small scales to larger scales. Typical crystal structures will be discussed systematically in this book to determine their elastic moduli.

In summary, such modeling tools will be more accurate and they can relate the structure of materials at the atomic or microscopic scale with their performance. It is not far for us to deeply investigate mechanical properties with the development of computer hardware level, the interatomic potential functions accurately described, new high-efficiency algorithms, etc.

Acknowledgements

The author acknowledges the financial support provided by the China Scholarship Council and the startup foundation of Xi’an Shiyou University. Thanks Qiufeng Wang for her proofreading, and very sincere thankness to prof. Fabrice Bernard and prof. Siham Kamali-Bernard for the their enthusiastic and responsible guideness.

The Advanced Atomic Finite Element Method: Modeling and Application

Abstract

An atomic finite element analysis is developed in this chapter. At atomic scale, the interatomic bonding forces of van der Waals and the covalent chemical bond are taken into account. In addition to the tradition atomistic simulation methods, the more specific objective of this chapter is to develop a modeling methodology at the atomic scale (or nanoscale) since the macroscopic properties largely depend on the physicochemical properties of the interatomic bonds. Above all, the methodology is applied to study the behavior of carbon nanotubes, whose development has experienced strong growth in recent years and can be used for quality mechanical reinforcement. These carbon nanotubes are formed by repeating zigzag carbon-carbon bonds. Development of atomic finite element method (AFEM) can be traced back to the homogenized elastic properties of various graphene structures (single-layer graphene sheet, zigzag single-walled carbon nanotubes, triple-layer graphene sheet). Moreover, the AFEM is developed to investigate the Young modulus of one of the main constituents of cement-based materials, portlandite (CH). Besides, a theoretical formula between axial modulus (Ex, Ey, and Ez) and elastic constants of hexagonal crystal is derived by crystal elastic theory. Modeling of single crystal can be traced back to the homogenized elastic properties of polycrystals.

Keywords: Portlandite, AFEM, DFT, anisotropy, Young’s modulus

1. Introduction

Composite material reinforced by carbon nanotubes (CNT) has a high strength, mechanical shock, thermal shock resistance, and fracture toughness, thus has a potential application with greatly improved properties. The “ene aluminum alloy”-graphene-reinforced aluminum matrix nanocomposites are firstly discovered for high properties in recent years. Carbon nanotube (CNT) and graphene have established a natural link; graphene was firstly discovered and separated by Geim [1] and the influential study shows that CNT can be formed by rolling a cylindrical graphene [2]. CNT has caused a widespread concern in civil engineering application. Due to disadvantages of lower tensile strength and lower impact strength of cement-based material, CNT as a reinforcing matrix material has higher mechanical stability properties and impact resistance. The difficulty in seeking the analytical solution of CNT Young’s modulus is existing due to choosing various wall thicknesses, zigzag or armchairs, D/L ratio, indirect measurements, dimensions, and modeling methods, etc. Young’s modulus of single-walled carbon nanotubes (SWCNT) is about 1.5–5.0 TPa [3] for individual CNT, 1.033–1.042 TPa [4] with thickness 0.34 nm by using nanoscale continuum mechanics method (CM), 5.5 TPa for SWCNT based on MD method with thickness 0.066 nm by Yakobson et al. [5], 4.88 TPa for SWCNT based on CM by Vodenitcharova and Zhang [6], 5.1 TPa for CNT based on the electronic energy band theory (EEBT) by Zhou et al. [7] with considering the total energy of all electrons. There is no agreement of Young’s modulus of CNT and graphene in the analytical method up to nowadays; the disagreement of Young’s modulus of CNT and graphene is actually due to the restriction of both the theory and experiment. Especially for graphene, atomic van der Waals forces between the layers cannot be ignored. As the measurement of a single CNT and graphene has some difficulty especially for a certain type (zigzag type or air-chairs type), it therefore causes our interest to present the exact forms of modeling at nanoscale.

It is common that graphene and nanotube structures sharing the same C▬C bond have a wide potential application. CNT has similar properties of fullerene molecules; its high strength has caused a widespread concern in various application fields. Graphene is considered to be an atomic crystal with flat polycyclic aromatic hydrocarbons (the structure is very stable). It can be rolled into a barrel-shaped carbon nanotube [2], showing the similar properties of nanotube and fullerene molecule. However, the theoretical mechanical deformation analysis is so complicated that the nonlinear analytical method is needed where van der Waals force should be considered and the corresponding nonlinear transformation is needed.

For the value of Young’s modulus of carbon nanotubes (CNT), after the separation of graphene by Geim [1], Young’s modulus is 1.033–1.042 TPa [8] with thickness 0.34 nm by using nanoscale continuum mechanics method, then 1.5–5.0 TPa for CNT by Treacy et al. [3], 5.5 TPa for SWCNT based on MD method with thickness 0.066 nm by Yakobson et al. [5], 4.88 TPa for SWCNT based on continuous mechanics method by Vodenitcharova and Zhang [6], 5.1 TPa for CNT based on the electronic energy band theory by Zhou et al. [7] with considering the total energy of all electrons. Researches of CNT are mainly listed in Table 1.

Year	Author	Method	t (nm)	E (TPa)	E.t
1996	Yakobson et al. [5]	MD	0.066	5.5	0.3630
1999	Goze et al. [9]	Tight binding	0.34	1.26	0.4284
2000	Popov et al. [10]	FCM	0.33	1.02	0.3368
2000	Zhou et.al [7]	Tight-binding MD	0.074	5.1	0.3774
2001	Kudin et al. [11]	Ab initio	0.089	3.859	0.3450
2002	Tu and Ou-Yang [12]	LDA	0.075	4.7	0.3525
2003	Li and Chou [8]	Structure mechanics	0.34	1.02	0.3463
2003	Chang and Gao [13]	Molecular mechanics	0.36	0.997	0.3592
2004	Sears and Batra [14]	Molecular mechanics	0.134	2.52	0.3377
2004	Pantano et al. [15]	Continuum shell	0.075	4.84	0.3630
2005	Zhang et al. [16]	C-B rule	0.335	1.08	0.3618
2006	To [17]	FEM	0.34	1.024	0.3482
2007	Chandraseker and Mukherjee [18]	Ab initio	0.335	0.99	0.3317
2008	Wang and Zhang [19]	Continuum shell	0.1	3.5	0.350
2008	Jinan Lu and Chen [20]	FEM	0.34	1.274	0.4332
2010	Mahmood et al. [4]	Continuum mechanics	0.34	1.033–1.042	0.3512–0.3543
2012	Xiaoxing and Zhong Hu [21]	Equivalent-continuum	0.33	1.058	0.349

Table 1.

Elastic modulus and wall thickness of nanotubes.

From Table 1, the Young’s moduli of CNT with considering various thicknesses using various methods are not all the same. Actually, there is no agreement of Young’s modulus nowadays due to choosing various wall thickness, zigzag or armchair, D/L ratio, indirect measurements, dimensions, methods, etc.

Thus, a new so-called atomic finite element method (AFEM) as a bridge between molecular dynamics and continuum mechanics [22] can be used to investigate elastic properties of covalent bond structures. In 2011, the AFEM [22] is proposed, in which the interatomic bonds are modeled as nonlinear spring elements. Besides, for the typical C▬C frame-like structures (e.g., graphene and SWCNT), a structural mechanics approach proposed by Li and Chou [8] has good accuracy and stability to calculate their Young’s modulus, and the modeling part of covalent bond can be modeled as an equivalent beam element. As a structural mechanics approach, Li and Chou [8] have successfully modeled the C▬C covalent bond structures by equivalent beam elements and a so-called AFEM [22] has solved the interatomic bonds modeled by nonlinear spring elements, modeling of nanostructure realized by an equivalent beam and nonlinear spring elements becomes possible.

In this chapter, carbon nanotubes (CNT), hexagonal graphene, and portlandite (CH) are particularly investigated. CNT, has its unique physical and chemical characteristics, has higher mechanical stability properties and impact resistance, thus can be used as a reinforcing material or composite material matrix. CH is a typical constituent of cement paste and its Young’s moduli are needed in the MuMoCC platform. Above all, Young’s modulus formula of zigzag CNT can be obtained and then the AFEM method is verified. Modeling of TLGSs structure and CH crystal using the equivalent beam-spring element is also investigated. This chapter aims to develop a modeling method of nanostructure realized by an equivalent beam and nonlinear spring elements, of which the nonlinear spring element’s definition of nonbonding can be realized according to a kind of Lennard-Jones potential in various corresponding equilibrium conditions.

2. Modeling methods of continuum mechanics and AFEM modeling

2.1 Classical modeling methods by continuum mechanics method

2.1.1 Classical molecular mechanics method (MM)

Generally, the force field is expressed in the form of steric potential energy. It depends solely on the positions of the nuclei constituting the molecule. The total steric potential energy, omitting the electrostatic interaction, is a sum of energies due to valence or bonded interactions and nonbonded interactions [23]:

$$U=\sum U_r+\sum U_{\theta }+\sum U_{\varphi }+\sum U_{\omega }+\sum U_{\mathit{vdw}}\hspace{1em}\mathrm{(1)}$$

where U_r is for a bond stretch interaction, U_θ for a bond angle bending, U_ϕ for a dihedral angle torsion, U_ω for an improper (out-of-plane) torsion, U_vdw for a nonbonded van der Waals interaction. Interatomic interactions in molecular mechanics by Li and Chou are shown in Figure 1.

Figure 1.

Interatomic interactions in MM by Li and Chou.

Li and Chou have proposed the molecular mechanics method (MM) [8], under the small deformation premise, to simplify the calculation, the simplest harmonic forms to merge the dihedral angle torsion and the improper torsion into a single equivalent term are adopted, which are as follows:

$$U_r=\frac{1}{2}{k}_r{\left(r-r_0\right)}^2=\frac{1}{2}{k}_r{\left(\mathrm{\Delta }r\right)}^2\hspace{1em}\mathrm{(2)}$$

$$U_{\theta }=\frac{1}{2}{k}_{\theta }{\left(\theta -{\theta }_0\right)}^2=\frac{1}{2}{k}_{\theta }{\left(\mathrm{\Delta }\theta \right)}^2\hspace{1em}\mathrm{(3)}$$

$$U_{\tau }=U_{\varphi }+U_{\omega }=\frac{1}{2}{k}_{\tau }{\left(\mathrm{\Delta }\varphi \right)}^2\hspace{1em}\mathrm{(4)}$$

where k_r, k_θ, and k_τ are the bond stretching force constant, bond angle bending force constant, and torsional resistance, respectively, and the symbols Δr, Δθ, and Δϕ represent the bond stretching increment, the bond angle change, and the angle change of bond twisting, respectively.

Assuming that the sections of carbon-carbon bonds are uniformly circular, thus $I_x=I_y=I$. So parameters of EA, EI, and GJ can be determined. As the molecules’ chemical energy is equivalent to the strain energy of beam, thus, there are the following approximate relationships:

$$\frac{\mathit{EA}}{L}=k_r,\frac{\mathit{EI}}{L}=k_{\theta },\frac{\mathit{GJ}}{L}=k_{\tau }\hspace{1em}\mathrm{(5)}$$

where L is the beam length of the C▬C bond, EA represents tensile stiffness, EI represents bending stiffness, GJ represents torsional stiffness of equivalent beam, k_r is the stretching constant, k_θ is the bending constant, and k_τ is the torsional constant. Obviously, these formulas have established the equivalent relation between the macrobeam structure and microscopic molecular structure. As long as the force constants k_r, k_θ, and k_τ are known, the sectional stiffness parameters EA, EI, and GJ can be readily obtained. Then, a mathematical modeling—AFEM modeling—of ionic crystals is proposed by means of interatomic potentials.

2.1.2 The chemical bond element method

The chemical bond element method has been proposed to link the relation between finite element method and continuum mechanics considering atomic interaction, where the local geometry of chemical bonds of most materials is characterized by a bond length r, bend angle θ, and dihedral angle φ, shown in Figure 2(a).

Figure 2.

Geometry of chemical bonds, nodal forces, and nodal displacements. (a) Geometry of chemical bonds. (b) Nodal forces and nodal displacements.

Chemical bond element nodal forces and displacements are shown in Figure 2(b). As is shown in Figure 2(b), there are three nodal forces for each node: F_xi, F_yi, M_i for node i and F_xj, F_yj, M_j for node j. There are three nodal displacements for each node: u_i, v_i, _θi for node i and u_j, v_j, _θj for node j.

The nodal force increment and the displacement increment in each loading step in the x-direction are:

$$\mathrm{\Delta }{F}_{\mathit{xi}}=k_r\left(\mathrm{\Delta }{u}_i-\mathrm{\Delta }{u}_j\right)\hspace{1em}\mathrm{(6)}$$

$$\mathrm{\Delta }{F}_{\mathit{xj}}=k_r\left(-\mathrm{\Delta }{u}_i+\mathrm{\Delta }{u}_j\right)\hspace{1em}\mathrm{(7)}$$

$$\mathrm{\Delta }{M}_i=-2k_{\mathit{\theta i}}\left(\mathrm{\Delta }{\theta }_L-\mathrm{\Delta }{\theta }_i\right)\hspace{1em}\mathrm{(8)}$$

$$\mathrm{\Delta }{M}_j=2k_{\mathit{\theta j}}\left(\mathrm{\Delta }{\theta }_j-\mathrm{\Delta }{\theta }_L\right)\hspace{1em}\mathrm{(9)}$$

where $\mathrm{\Delta }{\theta }_L=\frac{\mathrm{\Delta }{v}_j-\mathrm{\Delta }{v}_i}{L}$, k_r is the stretching stiffness, k_yi and k_yj are the bending stiffness of joints i and j, and L is the length of the element.

For the C▬C bonds, parameters by references [8, 24] are as follows: $k_r=2.78\mathrm{aJ}/{\mathring{\mathrm{A}}}^2$, $k_{\theta }=0.498\mathrm{aJ}/\mathrm{rad}$, $\epsilon =0.585\times {10}^{-3}\mathrm{aJ}$, $r_0=1.47\times {10}^{-10}\mathrm{m}$, ${\theta }_0=1.88\mathrm{rad}$, $\sigma =3.53\times {10}^{-10}\mathrm{m}$, $k_r=938\text{kcal}\times {\mathrm{mol}}^{-1}\times {\mathring{\mathrm{A}}}^{-2}$, $k_{\theta }=126\text{kcal}\cdot {\mathrm{mol}}^{-1}\cdot {\mathrm{rad}}^{-2}$, $k_{\tau }=40\text{kcal}\cdot {\mathrm{mol}}^{-1}\cdot {\mathrm{rad}}^{-2}$, $a_{C-C}=L=0.142\mathrm{nm}$. For the round section of C▬C bonds, when t = 0.34 nm, thus A = 0.090792 (nm²), E = 1019.259 (GPa), and I = 1.21959 × 10⁻⁴ (nm⁴).

2.2 AFEM realization at nanoscale considering atomic interaction

Continuum methods (such as FEM) use averaged description for modeling the material properties. As a result, they are not capable of modeling phenomena at the atomic scale. Localized nonlinear deformation, defects, nanoscale materials, and structures (e.g., CNTs) as well as various nanoscale phenomena need to be modeled with the atomic scale accuracy. CNT types and characteristic parameters are classified in Table 2.

Table 2.

Different types, chiral angle, and chiral vector of CNT.

Where chiral vector ${\overrightarrow{C}}_h=n{\overrightarrow{a}}_1+m{\overrightarrow{a}}_2$, $a_1=a_2=\sqrt{3}{a}_{c-c}$, so $C_h=\sqrt{3}{a}_{c-c}\sqrt{n^2+\mathit{nm}+m^2}$; chiral angle $cos\theta =\frac{2n+m}{2\sqrt{n^2+\mathit{nm}+m^2}}$ or $\theta ={\mathit{tg}}^{-1}\left[\frac{\sqrt{3}m}{m+2n}\right]$; diameter $d=\frac{C_h}{\pi }=\frac{\sqrt{3}}{\pi }{a}_{\mathit{cc}}\sqrt{n^2+\mathit{nm}+m^2}$.

2.2.1 AFEM modeling based on molecular mechanics

Molecular mechanics method mainly includes two basic assumptions. The first assumption is the Born-Oppenheimer approximation, of which the electronic function is centered on the nuclear positions, with the electrons around them in an optimal distribution [25]. The second assumption is that each nucleus is seen as a classical particle. Electrons are not taken into account and atoms are simulated as spheres, with an assigned radius and charge, interacting by means of a collection of interatomic potentials [26, 27].

The AFEM is proposed based on the development of nanotechnology across multiple length scales [22]. The AFEM can be used to calculate the van der Waals force in the plane. For a system composed of N atoms, the total energy of the system can be expressed as a sum of all the chemical energy stored:

$$E_{\mathit{\text{total}}}\left(x\right)=U_{\mathit{\text{total}}}\left(x\right)-\sum _{i=1}^N{\overline{F}}_i\cdot x_i\hspace{1em}\mathrm{(10)}$$

where $x_i$ is the coordinate of atom i vector and ${\overline{F}}_i$ is the external force exerted on the atom i vector. The minimum energy state is:

$$\frac{\partial E_{\mathit{\text{total}}}\left(x\right)}{\partial x}=0\hspace{1em}\mathrm{(11)}$$

For $\partial F_{\mathit{\text{total}}}\left(x\right)$, the initial state $x^{\left(0\right)}$ can be expanded by the developed Taylor, we obtain:

$$E_{\mathit{\text{total}}}\left(x\right)\approx E_{\mathit{\text{total}}}\left(x^{\left(0\right)}\right)+{\left.\frac{\partial E_{\mathit{\text{total}}}\left(x\right)}{\partial x}\right|}_{x=x^{\left(0\right)}}\cdot \left(x-x^{\left(0\right)}\right)+\frac{1}{2}{\left(x-x^{\left(0\right)}\right)}^T{\left.\frac{{\partial }^2{E}_{\mathit{\text{total}}}\left(x\right)}{\partial x\partial x}\right|}_{x=x^{\left(0\right)}}\left(x-x^{\left(0\right)}\right)\hspace{1em}\mathrm{(12)}$$

Each atom not only forms a direct bond with its three adjacent atoms, but also has an interaction with its six subadjacent atoms. Newton-Raphson iteration is as follows:

$$\psi \left(x^n\right)=-{\left(\frac{{\mathit{dE}}_{\mathit{\text{total}}}}{\mathit{dx}}\right)}_n+\overline{F}\hspace{1em}\mathrm{(13)}$$

$$\psi \left(x^{n+1}\right)=\psi \left(x^n\right)+{\left(\frac{\mathit{d\psi }}{\mathit{dx}}\right)}_n\mathrm{\Delta }{x}^n=0\hspace{1em}\mathrm{(14)}$$

According to Eq. (14), we can obtain:

$$\mathrm{\Delta }{x}^n=-{\left(\frac{\mathit{d\psi }}{\mathit{dx}}\right)}^{-1}\cdot \psi \left(x^n\right)\hspace{1em}\mathrm{(15)}$$

$$x^{n+1}=x^n+\mathrm{\Delta }{x}^n\hspace{1em}\mathrm{(16)}$$

where $\overline{F}$ is the external force vector, $\psi \left(x^n\right)$ is imbalance force vector after n iterations, $x^n$ is the vector of generalized coordinates of atoms after n iterations. Each iteration needs to compute the total potential energy of the system and coordinates of the first and the second derivatives.

For the expression of K and F, the corresponding expressions are as:

$$K=-\frac{\mathit{d\psi }}{\mathit{dx}}=\frac{d^2{E}_{\mathit{\text{total}}}}{{\mathit{dx}}^2}\hspace{1em}\mathrm{(17)}$$

$$F=-\frac{{\mathit{dE}}_{\mathit{\text{total}}}}{\mathit{dx}}+\overline{F}.\hspace{1em}\mathrm{(18)}$$

The above equation can be written as:

$$K\cdot u=P\hspace{1em}\mathrm{(19)}$$

For the atoms before deformation, the expression can be described as follows:

$$K={\left.\frac{{\partial }^2{E}_{\mathit{\text{total}}}}{\partial x\partial x}\right|}_{x=x^{\left(0\right)}}={\left.\frac{{\partial }^2{U}_{\mathit{\text{total}}}}{\partial x\partial x}\right|}_{x=x^{\left(0\right)}}\hspace{1em}\mathrm{(20)}$$

$$P={\left.\frac{\partial E_{\mathit{\text{total}}}}{\partial x}\right|}_{x=x^{\left(0\right)}}=\overline{F}-{\left.\frac{\partial U_{\mathit{\text{total}}}}{\partial x}\right|}_{x=x^{\left(0\right)}}\hspace{1em}\mathrm{(21)}$$

K and P can be directly obtained by the interatomic interaction potential. The appropriate boundary conditions and finite elements are chosen to calculate the minimum energy state during iterative calculation, and the coordinates of each atom in the system can be obtained.

Figure 3 shows the interaction of van der Waals force of atoms A and atoms L among atoms.

Figure 3.

Interaction of van der Waals force among atoms.

As is shown in Figure 3, stiffness matrix for unit atom A is thus expressed as:

$$K_A=\left[\begin{array}{cc}{\left(\frac{{\partial }^2{E}_{\mathit{\text{total}}}}{\partial a_A\partial a_A}\right)}_{2\times 2}& {\left(\frac{1}{2}\frac{{\partial }^2{E}_{\mathit{\text{total}}}}{\partial a_A\partial a_i}\right)}_{2\times 18}\\ {\left(\frac{1}{2}\frac{{\partial }^2{E}_{\mathit{\text{total}}}}{\partial a_i\partial a_A}\right)}_{18\times 2}& 0_{18\times 18}\end{array}\right]\hspace{1em}\mathrm{(22)}$$

The right endpoint of the vector of atom A can be expressed as:

$$F_A=\left[\begin{array}{l}{\left({\overline{F}}_A-\frac{\partial E_{\mathit{\text{total}}}}{\partial a_A}\right)}_{2\times 1}\\ 0_{18\times 1}\end{array}\right]\hspace{1em}\mathrm{(23)}$$

where $a_A$ represents the coordinates of atomic vector, and in the two-dimensional case, $a_A={\left[x_A{y}_A\right]}^T$; the range of i is B-J, and the sum of the order number of $K_A$ and $F_A$ is related to the interaction between atoms A and B.

For $F_A$, the component $\frac{\partial E_{\mathit{\text{total}}}}{\partial a_A}={\left[\frac{\partial E_{\mathit{\text{total}}}}{\partial x_A}\frac{\partial E_{\mathit{\text{total}}}}{\partial y_A}\right]}^T$; For $K_A$, the component $\frac{{\partial }^2{E}_{\mathit{\text{total}}}}{\partial a_A\partial a_A}$.

For component F_A:

$$\frac{\partial U\left(r_{\mathit{AL}}\right)}{\partial x_A}=\frac{\partial U\left(r_{\mathit{AL}}\right)}{\partial r_{\mathit{AL}}}\frac{\partial r_{\mathit{AL}}}{\partial x_A}\hspace{1em}\mathrm{(24)}$$

$$\frac{\partial U\left(r_{\mathit{AL}}\right)}{\partial y_A}=\frac{\partial U\left(r_{\mathit{AL}}\right)}{\partial r_{\mathit{AL}}}\frac{\partial r_{\mathit{AL}}}{\partial y_A}\hspace{1em}\mathrm{(25)}$$

For component K_A:

$$\frac{\partial U^2\left(r_{\mathit{AL}}\right)}{\partial x_A\partial x_A}=\frac{{\partial }^{2}U\left(r_{\mathit{AL}}\right)}{\partial r_{\mathit{AL}}^2}{\left(\frac{\partial r_{\mathit{AL}}}{\partial x_A}\right)}^2+\frac{\partial U\left(r_{\mathit{AL}}\right)}{\partial r_{\mathit{AL}}}\frac{{\partial }^2{r}_{\mathit{AL}}}{\partial x_A^2}\hspace{1em}\mathrm{(26)}$$

$$\frac{\partial U^2\left(r_{\mathit{AL}}\right)}{\partial y_A\partial y_A}=\frac{{\partial }^{2}U\left(r_{\mathit{AL}}\right)}{\partial r_{\mathit{AL}}^2}{\left(\frac{\partial r_{\mathit{AL}}}{\partial y_A}\right)}^2+\frac{\partial U\left(r_{\mathit{AL}}\right)}{\partial r_{\mathit{AL}}}\frac{{\partial }^2{r}_{\mathit{AL}}}{\partial y_A^2}\hspace{1em}\mathrm{(27)}$$

$$\frac{\partial U^2\left(r_{\mathit{AL}}\right)}{\partial x_A\partial y_A}=\frac{\partial U^2\left(r_{\mathit{AL}}\right)}{\partial y_A\partial x_A}=\frac{{\partial }^{2}U\left(r_{\mathit{AL}}\right)}{\partial r_{\mathit{AL}}^2}\frac{\partial r_{\mathit{AL}}}{\partial x_A}\frac{\partial r_{\mathit{AL}}}{\partial y_A}+\frac{\partial U\left(r_{\mathit{AL}}\right)}{\partial r_{\mathit{AL}}}\frac{{\partial }^2{r}_{\mathit{AL}}}{\partial x_A\partial y_A}\hspace{1em}\mathrm{(28)}$$

As is known to us, $U^{\mathit{VdW}}\left(r\right)=4\epsilon \left[{\left(\frac{\sigma }{r_{\mathit{ij}}}\right)}^{12}-{\left(\frac{\sigma }{r_{\mathit{ij}}}\right)}^6\right]$, so:

$$\frac{\partial U\left(r_{\mathit{AL}}\right)}{\partial r_{\mathit{AL}}}=4\epsilon \left[-12\times \frac{{\sigma }^{12}}{r_{\mathit{AL}}^{13}}+6\times \frac{{\sigma }^6}{r_{\mathit{AL}}^7}\right]\hspace{1em}\mathrm{(29)}$$

$$\frac{{\partial }^{2}U\left(r_{\mathit{AL}}\right)}{\partial r_{\mathit{AL}}^2}=4\epsilon \left[156\times \frac{{\sigma }^{12}}{r_{\mathit{AL}}^{14}}-42\times \frac{{\sigma }^6}{r_{\mathit{AL}}^8}\right]\hspace{1em}\mathrm{(30)}$$

$$\frac{\partial r_{\mathit{AL}}}{\partial x_A}=\frac{x_A-x_L}{r_{\mathit{AL}}}\hspace{1em}\mathrm{(31)}$$

$$\frac{{\partial }^2{r}_{\mathit{AL}}}{\partial x_A^2}=\frac{1}{r_{\mathit{AL}}}-\frac{x_A-x_L}{r_{\mathit{AL}}^2}\frac{\partial r_{\mathit{AL}}}{\partial x_A}\hspace{1em}\mathrm{(32)}$$

$$\frac{\partial r_{\mathit{AL}}}{\partial y_A}=\frac{y_A-y_L}{r_{\mathit{AL}}}\hspace{1em}\mathrm{(33)}$$

$$\frac{{\partial }^2{r}_{\mathit{AL}}}{\partial y_A^2}=\frac{1}{r_{\mathit{AL}}}-\frac{y_A-y_L}{r_{\mathit{AL}}^2}\frac{\partial r_{\mathit{AL}}}{\partial y_A}\hspace{1em}\mathrm{(34)}$$

$$\frac{{\partial }^2{r}_{\mathit{AL}}}{\partial x_A\partial y_A}=\frac{{\partial }^2{r}_{\mathit{AL}}}{\partial y_A\partial x_A}=-\frac{\left(x_A-x_L\right)\left(y_A-y_L\right)}{r_{\mathit{AL}}^2}\hspace{1em}\mathrm{(35)}$$

In summary, there are two types of interactions, including bonded and nonbonded: (1) for bonded type, it usually corresponds to strong covalent bonds, the number of neighbors is preset and is dictated by the lattice. Take the CNT for example, the carbon covalent bond often has carbon atoms 3–4 neighbors, and preset list of neighbors for each atom may be used throughout the calculations, used for crystalline solids with well-defined bonds. Typical structure example is graphene; (2) for nonbonded type, it commonly corresponds to relatively weak forces, and the number of neighbors may be changing and is found based on the cutoff radius. It can be used to represent van der Waals and Coulomb forces of solid/crystal unit.

Background of atomic finite element method (AFEM) contains the C▬C covalent bonds [8] by beam element and the interatomic bonds modeled [22] by nonlinear spring elements of AFEM, shown in Figure 4.

Figure 4.

AFEM modeling hypotheses of covalent bond and ionic bond. (a) Covalent bond. (b) Ionic bond or metallic bond.

From Figure 4, the ionic bond or metallic bond is taken into consideration, van der Waals interaction and electrostatic interaction (Coulomb force) can be separately calculated. Besides, modeling of covalent bond is based on the theory of molecular mechanics method (MM). AFEM can be readily linked with the conventional continuum FEM since they are in the same theoretical framework. Based on the AFEM modeling, the elastic stiffness tensors can be described and defined [28].

2.2.2 Description of AFEM modeling hypotheses

The expression of the total steric potential energy U is given out [23]. For ionic crystal, the Coulomb force should also be considered. So the ideal expression of U, considering the electrostatic interaction, is a sum of energies due to valence or bonded interactions and nonbonded interactions:

$$U=\sum U_r+\sum U_{\theta }+\sum U_{\varphi }+\sum U_{\omega }+\sum U_{\mathit{vdw}}+\sum U_{\mathit{coul}}\hspace{1em}\mathrm{(36)}$$

where U_r is for a bond stretch interaction, U_θ for a bond angle bending, U_ϕ for dihedral angle torsion, U_ω for an improper (out-of-plane) torsion, U_vdw for a nonbonded van der Waals interaction, and U_coul for Coulombic (electrostatic) interactions. Equivalent beam elements [8] can be used to represent covalent bond to calculate U_r, U_θ, U_ϕ, and U_ω of Eq. (36), where the molecules’ chemical energy is equivalent to the beam strain energy. The coupled field of van der Waals force and Coulomb force between atoms or charged ions may be obtained by using Lennard-Jones potential or other potentials.

Nonlinear spring elements can represent U_vdw and U_coul of ionic bond. Around the equilibrium position, van der Waals force is developed by the first-order Taylor:

$$F^{\mathit{VdW}}\left(r_i\right)=F^{\mathit{VdW}}\left(r_0\right)+\frac{\mathit{dF}\left(r_0\right)}{{\mathit{dr}}_i}\left(r_i-r_0\right)=25\frac{\epsilon }{\sigma }\left[2{\left(\frac{\sigma }{r_0}\right)}^{13}-{\left(\frac{\sigma }{r_0}\right)}^7\right]-24\frac{\epsilon }{\sigma }\left[26{\left(\frac{\sigma }{r_0}\right)}^{14}-7{\left(\frac{\sigma }{r_0}\right)}^8\right]\left(r_i-r_0\right)\hspace{1em}\mathrm{(37)}$$

where r₀ (unit: Å) is the distance between two atoms around the equilibrium position without external force, r_i is the distance of the Lennard-Jones potential interaction, ε is the bond energy of van der Waals interaction between two balanced atoms, and σ is the value of r at which the potential is zero [29].

For arbitrary atom type of A-B, the van der Waals force can be obtained as follows:

$$F_{\mathit{VdW}}^{\mathit{\text{atomA}}-\mathit{\text{atomB}}}\left(r\right)=24\frac{{\epsilon }^{A-B}}{{\sigma }^{A-B}}\left[2{\left(\frac{{\sigma }^{A-B}}{r^{A-B}-r_0^{A-B}}\right)}^{13}-{\left(\frac{{\sigma }^{A-B}}{r^{A-B}-r_0^{A-B}}\right)}^7\right]-F_{\mathit{VdW}}^{\mathit{\text{atomA}}-\mathit{\text{atomB}}}\left(r_0^{A-B}\right)\hspace{1em}\mathrm{(38)}$$

For U_vdw, the nonlinear displacement-force characteristic of spring element as the first derivative of potentials must be defined by means of the program potentials [22], then by the coordination transformation.

Coulomb potential of ionic crystal according to Coulomb’s law is given by:

$$U_{\mathit{ij}}^{\mathit{\text{Coulom}}b}=q_i{q}_j\left(r_i-r_j\right)/\left(4{\mathit{\pi \epsilon }}_0\cdot {\left|r_i-r_j\right|}^2\right)\hspace{1em}\mathrm{(39)}$$

where q_i and q_j are partial charges, ε₀ is the dielectric permittivity of vacuum (8.85419 × 10⁻¹² F/m).

Nonlinear spring element is used to represent interatomic force of Lennard-Jones potential. For U_vdw and U_coul represented by nonlinear spring element, since transforming from $\left[x,y,1\right]$ to $\left[{\mathrm{x}}^{\prime },{\mathrm{y}}^{\prime },1\right]$ by the matrix [T], the expression is as:

$$\left[{\mathrm{x}}^{\prime },{\mathrm{y}}^{\prime },1\right]=\left[x,y,1\right]\cdot T=\left[x,y,1\right]\cdot \left[\begin{array}{c}1,0,0\\ 0,1,0\\ r_0,F_0,0\end{array}\right]=\left[x+r_0,y+F_0,1\right]\hspace{1em}\mathrm{(40)}$$

where r₀ and F₀ are the offset values in the x- and y-directions, respectively.

2.2.3 Interactions between atoms by Lennard-Jones potential

Generally, the expression of van der Waals force can be obtained by the Lennard-Jones potential [30]:

$$U_{\mathit{ij}}^{\mathit{VdW}}\left(r\right)=U_{\mathit{repuls\hspace{0.17em}if}}+U_{\mathit{attract\hspace{0.17em}if}}=\sum _i\sum _{j}4{\epsilon }_{\mathit{ij}}\left[{\left(\frac{{\sigma }_{\mathit{ij}}}{r_{\mathit{ij}}}\right)}^{12}-{\left(\frac{{\sigma }_{\mathit{ij}}}{r_{\mathit{ij}}}\right)}^6\right]\hspace{1em}\mathrm{(41)}$$

where r_ij is the distance between two atoms i and j; ε_ij and σ_ij are parameters of Lennard-Jones. ε_ij represents the minimum potential energy between atoms i and j, which corresponds to the most stable interaction (the depth of the potential well under the distance of σ_ij). σ_ij is the equilibrium distance where atomic energy between two atoms is zero.

Schematic illustration of Lennard-Jones potential and van der Waals force is shown in Figure 5.

Figure 5.

Schematic illustration of Lennard-Jones potential and van der Waals force. (a) Schematic of Lennard-Jones potential. (b) van der Waals force of C▬C bond.

As shown in Figure 5(a), repulsive force is more intense but the range is very short that prevents exchange of electrons from occupying the same region of space, since the forces exerted at larger distances are very small and can be neglected [23, 31]. Besides, from Figure 5(b), the equilibrium position of Lennard-Jones potential is at the place of r/σ = 1. The potential must be truncated at some point named the cutoff radius, R_cut, usually set to R_cut = 2.5σ. When the atomic distance is greater than 2.5σ, the Lennard-Jones potential is approximate to zero. As is shown in Figure 5(b), if the average distance r between the atoms increases, the attractive force of van der Waals force between pairs of atoms will resist the force caused by the increase of distance r. And if the average distance r between the atoms increases, the repulsive force between pairs of atoms will be able to resist the force caused by the decrease of distance r.

Interactions of van der Waals consist of two terms: a repulsive term and an attractive term. The function of Lennard-Jones is known as a function of the potential of van der Waals type. It is used to describe the intermolecular force between argon atoms i and j separated by a distance r_ij. Any two molecules attract each other at long separation distance and repel each other when they come closer [32].

Moreover, coordinate transformation is used, and these curves are translated so that the system at the beginning of the calculation is at the equilibrium (no more attractive/repulsive force). When the interatomic spacing is greater than its unstressed value, the attractive forces between atoms must be greater than the repulsive forces (the attractive forces balance both the repulsive forces and the external forces). Repulsive forces are taken as positive and attractive as negative. This conveniently makes their potential energies positive and negative, respectively, in which the zero of potential energy is taken at large separation. Thus, AFEM can be used in multiscale computation with applications to carbon nanotubes and CH [28, 33].

2.2.4 Modeling of ionic crystal based on AFEM method

AFEM modeling of zigzag CNT and zigzag graphene thus can be described by the algorithm, in Figure 6.

Figure 6.

Modeling of hexagonal structure based on AFEM and solving process.

As is shown in Figure 6, covalent bonds (e.g., for C▬C in the studied structures) are characterized by an equivalent beam, and nonbonded potential is characterized by using nonlinear spring elements to represent the coupled field of van der Waals force. Then, axial moduli of the structure, as a function of elastic constants, are derived and verified by cubic linear elements volume unit on ABAQUS. Thus, elastic moduli of TLGSs and CH hexagonal structures can be determined.

In the AFEM modeling, O▬H bond and C▬C bond as covalent bonds are characterized to be equivalent beam based on MM method, and nonlinear spring element is used to characterize the coupled field of van der Waals force and Coulomb force for ionic or metallic bond by using Lennard-Jones potential or Buckingham potential, which is calculated by MATLAB and then imported into input program. Moreover, axial moduli of triple-layer graphene sheet (TLGSs) and CH hexagonal structure based on AFEM as a function of elastic constants are derivated and verified by FE simulation on ABAQUS. Elasticity coefficient matrix of TLGSs and CH crystals are solved reversely, thus elastic properties of TLGSs and CH polycrystals are furtherly investigated by the well-known Reuss-Voigt-Hill estimation.

2.3 Elastic modulus determination of hexagonal crystal structures

ABAQUS as a powerful nonlinear software is used to the continuation of quantum mechanics and an approximate equivalent method—AFEM. Nonlinear springs are introduced to describe the relationship between atoms of different initial distances. When the van der Waals force is considered, the coordinate transformation needs to be done before carrying out the INP file. Input file usage:

*SPRING, NONLINEAR, DEPENDENCIES = n

first data line

force, relative displacement, field variable 1, etc.

…

The nonlinear spring behavior is not supported in ABAQUS/CAE when you define springs as engineering features. Instead, we can define connectors that have spring-like elastic behavior. To obtain the homogenized moduli, the notion of homogenization has been introduced using the representative volume element (RVE) of the material [34]. When a solid model is under external load, microstresses and strains are induced; according to the statistical homogeneity assumption, an appropriate RVE can be defined and isolated. On the RVE boundary, there exist definitive surface displacements and surface tractions. Within the RVE, there exist definitive stress field σ_ij and strain field ε_ij. Through homogenization [35], mechanical behavior is described by a definitive constitutive law. So, the mechanical properties of anisotropic solids can be investigated using FEM simulation by RVE unit on ABAQUS.

2.3.1 Axial moduli determination of hexagonal structures

Stress-strain relation in an orthotropic hexagonal crystal can be defined by the independent elastic stiffness parameters of ABAQUS manual [36]:

$$\begin{array}{c}\left[\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{12}\\ {\sigma }_{13}\\ {\sigma }_{23}\end{array}\right]=\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& 0& 0& 0\\ & c_{11}& c_{13}& 0& 0& 0\\ & & c_{33}& 0& 0& 0\\ & \text{sym}& & c_{44}& 0& 0\\ & & & & c_{44}& 0\\ & & & & & 1/2\times \left(c_{11}-c_{12}\right)\end{array}\right]\left[\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ {\epsilon }_{33}\\ {\gamma }_{12}\\ {\gamma }_{13}\\ {\gamma }_{23}\end{array}\right]\\ =\left[\begin{array}{cccccc}{D}_{1111}& D_{1122}& D_{1133}& 0& 0& 0\\ & D_{2222}& D_{2233}& 0& 0& 0\\ & & D_{3333}& 0& 0& 0\\ & \text{sym}& & D_{1212}& 0& 0\\ & & & & D_{1313}& 0\\ & & & & & D_{2323}\end{array}\right]\left[\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ {\epsilon }_{33}\\ {\gamma }_{12}\\ {\gamma }_{13}\\ {\gamma }_{23}\end{array}\right](42)\end{array}$$

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. The formula developed by finite element technique to solve the elastic constants by volume unit has been described [28], as in a manual of theory details [37].

If a unidirectional force is applied to obtain the strain value, the elastic constants of the matrix can be obtained by simultaneous equations. Figure 7 shows schematic deformation component.

Figure 7.

Schematic deformation component and loading modes of linear/shear strain. (a) Volume unit and stress direction under unidirectional stress. (b) Loading modes of linear/shear strain.

A finite element method for solving the elastic coefficient matrix C_ij is proposed according to the schematic volume unit and stress direction after strain-loading (Figure 7(a)) and the different components of the deformation through different loading modes (Figure 7(b)). Supposing $A={\epsilon }_{11},B={\epsilon }_{22},C={\epsilon }_{33},D={\epsilon }_{31},E={\epsilon }_{11}+{\epsilon }_{21}$ (here ${\epsilon }_{\mathit{ij}}$ independent strain of each unit in all directions), as elastic constants C_ij = D^el, Eq. (42) can be changed as:

$$\left(\begin{array}{ccccc}A& B& C& 0& 0\\ B& A& C& 0& 0\\ 0& 0& A+B& C& 0\\ 0& 0& ED& 0& \\ 0& 0& 0& 0& {\gamma }_{12}\end{array}\right)\left(\begin{array}{c}{c}_{11}\\ c_{12}\\ c_{13}\\ c_{33}\\ c_{44}\end{array}\right)=\left(\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ 0\\ {\sigma }_{12}\end{array}\right)\to \left\{\begin{array}{l}{c}_{11}=D_{1111}=D_{2222}=\frac{A{\sigma }_{11}-B{\sigma }_{22}}{A^2-B^2}-\frac{\mathit{CD}{\sigma }_{33}}{D{\left(A+B\right)}^2-C\left(A+B\right)E}\\ c_{12}=D_{1122}=\frac{A{\sigma }_{22}-B{\sigma }_{11}}{A^2-B^2}-\frac{\mathit{CD}{\sigma }_{33}}{D{\left(A+B\right)}^2-C\left(A+B\right)E}\\ c_{13}=D_{1133}=D_{2233}=\frac{D{\sigma }_{33}}{D\left(A+B\right)-\mathit{CE}};\\ c_{33}=D_{3333}=-\frac{E{\sigma }_{33}}{D\left(A+B\right)-\mathit{CE}};\\ c_{44}=D_{1212}=D_{1313}=\frac{{\sigma }_{12}}{{\gamma }_{12}};c_{66}=\frac{1}{2}\left(c_{11}-c_{12}\right)=D_{2323}\end{array}\right.\hspace{1em}\mathrm{(43)}$$

Stress-strain relation of hexagonal crystal is presented in Eq. (42). Assuming a normal stress σ₃ in z-axis direction is introduced, it will not result in shear strains, thus: γ₁₂ = γ₂₃ = γ₃₁ = 0. We assume σ₃₃ expressed by ε₃₃ as: σ₃₃ = E₃ε₃₃. Substituting this equation into the ternary equations of stress-strain relation Eq. (42), if ε₁₃, ε₂₃, ε₃₃ are not simultaneously zero, the determinant factor must equal to zero according to the theory of linear algebra, so homogeneous equations can be changed:

$$\left\{\begin{array}{l}0=c_{11}{\epsilon }_{13}+c_{12}{\epsilon }_{23}+c_{13}{\epsilon }_{33}\\ 0=c_{12}{\epsilon }_{12}+c_{11}{\epsilon }_{23}+c_{13}{\epsilon }_{33}\\ 0=c_{13}{\epsilon }_{13}+c_{13}{\epsilon }_{23}+\left(c_{33}-E_3\right){\epsilon }_{33}\end{array}\to \left|\begin{array}{ccc}{c}_{11}& c_{12}& c_{13}\\ c_{12}& c_{11}& c_{13}\\ c_{13}& c_{13}& c_{33}-E_3\end{array}\right|=0\right.\hspace{1em}\mathrm{(44)}$$

Eq. (44) can be simplified as:

$$\begin{array}{l}\left|\begin{array}{ccc}{c}_{11}& c_{12}& c_{13}\\ c_{12}& c_{11}& c_{13}\\ c_{13}& c_{13}& c_{33}-E_3\end{array}\right|\underset{\_}{\underset{\_}{R_2-R_1}}\left|\begin{array}{ccc}{c}_{11}& c_{12}& c_{13}\\ c_{12}-c_{11}& c_{11}-c_{12}& 0\\ c_{13}& c_{13}& c_{33-E_3}\end{array}\right|\underset{\_}{\underset{\_}{C_2+C_1}}\left|\begin{array}{ccc}{c}_{11}& c_{11}+c_{12}& c_{13}\\ c_{12}-c_{11}& 0& 0\\ c_{13}& 2c_{13}& c_{33}-E_3\end{array}\right|\hfill \\ =\left|\begin{array}{cc}{c}_{11}+c_{12}& c_{13}\\ 2c_{13}& c_{33-E_3}\end{array}\right|=0.\hfill \end{array}$$

Then, we can get the relationship between z-axial modulus E₃ and elastic constants C_ij, the formula of z-axial modulus E₃ is thus obtained:

$$E_3=c_{33}-\frac{2c_{13}^2}{c_{11}+c_{12}}\hspace{1em}\mathrm{(45)}$$

Similarly, we can obtain relationship between y-axial modulus E₂ and elastic constants as:

$$E_2=\left(c_{11}-c_{12}\right)\left(1+\frac{c_{12}{c}_{33}-c_{13}^2}{c_{11}{c}_{33}-c_{13}^2}\right)\hspace{1em}\mathrm{(46)}$$

Finally, we can see that formulas of E₁ and E₂ axial moduli in x- and y-directions are the same, so:

$$E_1=E_2=\left(c_{11}-c_{12}\right)\left[1+\left(c_{12}{c}_{33}-c_{13}^2\right)/\left(c_{11}{c}_{33}-c_{13}^2\right)\right]\hspace{1em}\mathrm{(47)}$$

Based on elastic constants calculated, homogenized elastic properties of polycrystals can be calculated. Firstly, shear and bulk moduli are obtained by means of Voigt and Reuss bounds, which are based on the c_ij calculated by Wu et al. [38]. Then, Young’s modulus can be finally determined from shear/bulk modulus formulas.

2.3.2 Elastic coefficient solution by finite element method

As ${\epsilon }_1={\epsilon }_a,{\epsilon }_3={\epsilon }_c$, so ${\epsilon }_b$ can be easily obtained by the relative displacement in b-direction.

$$\overrightarrow{b}\left(\begin{array}{c}-cos{60}^{°}\\ sin{60}^{°}\\ 0\end{array}\right)={\displaystyle \overrightarrow{b}}\left(\begin{array}{c}-1/2\\ \sqrt{3}/2\\ 0\end{array}\right),{\epsilon }_b=\left(\begin{array}{c}-1/2\\ \sqrt{3}/2\\ 0\end{array}\right)\cdot \left(\begin{array}{ccc}{c}_{11}& c_{12}& c_{13}\\ c_{12}& c_{11}& c_{13}\\ c_{13}& c_{13}& c_{33}\end{array}\right)\cdot \left(\begin{array}{c}-1/2\\ \sqrt{3}/2\\ 0\end{array}\right)\hspace{1em}\mathrm{(48)}$$

For hexagonal structure, the coordinate transformation from 90 to 120° (in Figure 8), the relation is as:

Figure 8.

Solving equivalent strain conversion from 90 to 120°.

Under uniaxial condition (${\epsilon }_{12}={\epsilon }_{13}={\epsilon }_{23}=0$), so

$${\epsilon }_b=\frac{1}{4}{\epsilon }_{11}+\frac{3}{4}{\epsilon }_{22},{\epsilon }_{22}=\frac{4}{3}{\epsilon }_b-\frac{1}{3}{\epsilon }_{11}\hspace{1em}\mathrm{(49)}$$

Linear elasticity in an orthotropic material can also be defined by giving the independent elastic stiffness parameters, as functions of predefined fields. ABAQUS/standard spatially varying orthotropic elastic behavior can be defined for homogeneous solid continuum elements by using a distribution. The distribution must include default values for the elastic moduli. If a distribution is used, no dependencies on temperature and/or field variables for the elastic constants can be defined and determined on ABAQUS.

3. Application to structures of zigzag CNT and zigzag SLGSs

For application to structures of CNT and graphene by AFEM, a calculation formula of Young’s modulus of zigzag SWCNT and zigzag SLGSs is given out and verified. As the total potential energy is equal to the thin cylindrical shell strain energy, the elastic modulus of (n, 0) zigzag nanotube is developed by reference [13] without considering the thickness as:

$$E_{\mathit{CNT}}^{\mathit{\text{zigzag}}}=\frac{4\sqrt{3}\mathit{KC}}{9C+4{\mathit{Kl}}_0^2\left({\lambda }_{z1}^2+2{\lambda }_{z2}^2\right)}\hspace{1em}\mathrm{(50)}$$

where K, C, respectively, correspond to constants of bond angle changes and bond stretching; λ_z1 and λ_z2 are structure parameters of zigzag nanotube, ${\lambda }_{z1}=-3\sqrt{4-3{cos}^2\left(\frac{\pi }{2n}\right)}\cdot cos\left(\frac{\pi }{2n}\right)/\left[8\sqrt{3}-2\sqrt{3}{cos}^2\left(\frac{\pi }{2n}\right)\right]$, and ${\lambda }_{z2}=12-9{cos}^2\left(\frac{\pi }{2n}\right)/\left[16\sqrt{3}-4\sqrt{3}{cos}^2\left(\frac{\pi }{2n}\right)\right]$. Young’s modulus formulas of the graphene sheet can be derived when n becomes ∞ in zigzag nanotube parameters λ_z1 and λ_z2 expression, which is as follows:

$$E_{\mathit{\text{SLGSs}}}^{\mathit{\text{zigzag}}}=\frac{8\sqrt{3}\mathit{KC}}{18C+{\mathit{Kl}}_0^2}\hspace{1em}\mathrm{(51)}$$

where K, C, respectively, correspond to constants of bond angle changes and bond stretching. So, the Young’s modules of C▬C covalent bonds of zigzag SWCNT and SLGSs can be verified by reference [13].

Axial direction restraint and transverse direction with freedom (to ensure axial cross section during tensile can be freely deformed with shrinkage), the concentrated tensile force is added on each node of the other end. The sum of axial centralized force vectors at each node of carbon nanotubes is added up to F, and the finite element calculation is continuously carried out on ABAQUS. Young’s modulus is obtained as:

$$E=\frac{\sigma }{\epsilon }=\frac{F/A}{\mathrm{\Delta }l/l}=\frac{F\cdot l}{\mathit{\pi D}\cdot t\cdot \mathrm{\Delta }l}\hspace{1em}\mathrm{(52)}$$

where A is the cross-sectional area, A = πDt; D is the carbon nanotube diameter, and t is the wall thickness. The boundary conditions are free-fixed.

3.1 Nanoscale modeling of zigzag SWCNT and SLGSs

3.1.1 Definition of C▬C bond by coordination transformation

The transformed van der Waals curves are calculated, and the one-spring and two-spring of VDW forces are tested on ABAQUS software, as in Figure 9.

Figure 9.

One-spring and two-spring tests of transformed van der Waals curves. (a) True VDW of one-spring test. (b) True VDW of two-spring test.

From Figure 9(a) and (b), the definition of van der Waals forces and the coordinate conversion of C▬C bond are tested on ABAQUS, spring step-election △x is 0.001 nm in order to meet the accuracy requirements. For C▬C bond, parameters are as: r_max = 0.425 nm, F_max = 7.0174 nN. When the distance r = 0.3833208 nm, the van der Waals force is zero, whereas the equilibrium distance r₀ is 0.3415 nm, and the F₀ is −3.9247 nN.

3.1.2 Modeling of zigzag CNT and zigzag SLGSs

Modeling of a finite number graphene and (10, 0) zigzag CNT are shown in Figure 10.

Figure 10.

AFEM modeling of zigzag-CNT and single-layer graphene sheet (SLGSs). (a) Structures of zigzag graphene and SWNCT. (b) Models of zigzag graphene and SWNCT.

The boundary conditions of SLGSs and zigzag CNT with the fixed bound at x-direction and the tensional force applied at y-direction are shown in Figure 10(a). Moreover, the size 3.339 nm × 4.446 nm of SLGSs and (10, 0) zigzag CNT modeled with ABAQUS are shown in Figure 10(b) to apply the tensional force in y-direction. Parameters of C▬C covalent bonds used in simulation are as [13, 39]: C▬C bond length is 0.142 nm, thickness is 0.3414 nm, bond angle θ₀ is 120°, balanced distance constant σ_c-c is 0.3414 nm, and bond energy constant ε_c-c is 0.00239 eV. After we obtain the reaction force of y-direction, we substitute all parameters into the formula of E_SLGSs and E_CNT in Eq. (52), the results are as: E_SLGSs = 0.867 TPa, E_CNT = 0.847 TPa.

3.2 Young’s moduli of zigzag SWCNT and SLGSs

Here, we choose K = 742 nN/nm and C = 1.42 nN.nm; When n becomes ∞ (here λ_z1 = −0.28867; λ_z2 = 0.14433), Young’s modulus of graphene by Eq. (51) is calculated. Figure 11 shows the Young’s modulus variation of zigzag SWCNT with various diameters and the thickness and Young’s modulus relation.

Figure 11.

Young’s modulus variation and the thickness and Young’s modulus relation of zigzag SWCNT. (a) Young’s modulus variation with various diameters. (b) The thickness and Young’s modulus relation.

Figure 11(a) provides the result of the present modeling of zigzag SWCNT and SLGSs with a thickness of 0.34 nm, and the comparison of present results (E_SLGSs = 1.06 TPa) and literature results (E_CNT = 1.05 TPa [40]) for various thicknesses are also added. It can be seen in Figure 11(b) that the result of E.t (product of the Young’s modulus and the thickness) tends to the analytical value found by molecular mechanics approach [13]. Young’s modulus variation of zigzag nanotubes with the diameter is shown in Figure 11, the result is in good agreement with zigzag nanotubes [8, 9, 10, 13, 14, 20, 21] and no-distinguish CNT [5, 11, 12, 15, 16, 17, 18]. When the thickness of carbon nanotubes chooses to be the experience value 0.34 nm, the elasticity modulus E of carbon nanotubes calculated tends to 1.05 TPa, which is close to the consistent results of Yakobson et al. [5] and the result of the force-constant models by Popov et al. [10] and an analytical model based on a molecular mechanics approach by Chang and Gao [13]. The main difference with references [8, 13, 20] stems from the different values of the chemical bond parameters {k_r, k_θ, k_τ}. It can be seen in Figure 11(b) that the numerical result of E.t of zigzag CNT with various thicknesses is plotted, we can see that E.t is mainly distributed in the range of 0.33–0.37 kJ/m², and our result of E.t is also consistent with literature findings. It therefore shows that the atomic finite element analysis to study a covalent bond is feasible.

4. Application to the nanostructure of the triple-layer graphene

In general, the configuration of the equilibrium solid is decided by the minimum energy state. The basic idea of the continuum finite element is to divide into a certain number of units, of which the cell is represented by the nodes of the finite number of discrete. The location of these nodes is determined by the energy minimization state, which is similar to the basic idea of molecular statics to calculate the position of the atoms. The point is to minimize the energy of molecular static system for determining the position of each atom, that is why the concept of cohesion is introduced to determine the displacement of each point. For multilayered graphene, two interatomic potentials are used [41, 42]: one is for in-layer bonded interactions—multibody (Tersoff-Brenner), the other for interlayer nonbonded interactions (Lennard-Jones). A multilayer graphene modeling considering van der Waals force can be described by introducing nonlinear spring element between layers, which can be modeled by the AFEM method.

4.1 C▬C bond definition and nanoscale modeling of graphene

AFEM modeling of triple-layer graphene sheet (TLGSs) is shown in Figure 12.

Figure 12.

AFEM modeling of triple-layer graphene sheet (TLGSs). (a) Graphene crystal. (b) TLGSs model with beam-spring element. (c) van der Waals force of C▬C bond.

Graphene crystal is shown in Figure 12(a). A round-section definition of equivalent graphene-beam is: k_r = 669 nN/(nm²), k_θ = 0.876 nN.nm, and k_τ = 0.278 nN.nm [8]. The size 3.41 × 3.46 × 0.680 nm of TLGSs is modeled with ABAQUS in Figure 12(b) to determine its elastic constants. van der Waals force of C▬C bond is shown in Figure 12(c). Parameters of Lennard-Jones potential is as [30]: $U_{\mathit{vdw}}=-15.2/r^6+24.1\times {10}^3/r^{12},r_0=\sqrt[6]{2}\sigma =0.383\mathit{nm}$. Parameters of graphene used in simulation are as [24, 43]: C▬C bond length is 0.142 nm, thickness is 0.3414 nm, bond angle θ₀ is 120°, balanced distance constant σ_c-c is 0.3414 nm, and bond energy constant ε_c-c is 0.00239 eV.

C▬C bonds are represented by means of equivalent beams, for triple-layer graphene sheets (TLGSs) model, interlayer is mainly maintained by van der Waals force, thus it can be defined by Lennard-Jones potential and connected by nonlinear-spring element. The cohesive energy between layers is also considered [44]. Here the C-C bond distance σ = 0.3415 nm, spacing distance is as: h₁ = 0.992σ_C-C = 0.3373 nm, h₂-h₁ = 0.997σ_C-C = 0.3390 nm.

4.2 Elastic constants and homogenized moduli of graphene

Similarly, elastic constants of TLGSs are verified by volume unit. Then, we extend to the equivalent beam-spring model by AFEM. Diagram of strain loading in triple-layer graphene sheets is as shown in Figure 13.

Figure 13.

Strain loading of triple-layer graphene sheets.

Strain loading of triple-layer graphene sheets is shown in Figure 13. Then, the strain-stress data of TLGSs with beam element and beam-spring element separately by FEM and AFEM are given in Table 3.

Strain-stress by volume unit anisotropy				Strain-stress by beam-spring element
σij (nN/nm^2)	ε11	ε22	ε33	σij (nN/nm^2)	ε11	ε22	ε33
σ11 = 902.32	0.879739	−0.145212	−0.301956	σ11 = 957.769	0.8797	−0.15202	−0.112044
σ22 = 889.883	−0.143092	0.8676238	−0.297795	σ22 = 944.594	−0.14993	0.86762	−0.112044
σ33 = 31.8732	−0.0106694	−0.0106694	0.882005	σ33 = 15.641	−5.278e-3	−5.403e-3	0.441003
σ12 = 3.0591	γ12 = 0.6294			σ12 = 2.256	γ12 = 0.8797

Table 3.

Strain-stress data of grapheme by volume unit using FEM and beam-spring element using AFEM modeling.

From Table 3, C_ij is calculated and axial modulus by Eqs. (45)–(47) is as: σ_x = σ_y = 1025.21ε_y, σ_z = 36.14ε_z. As the graphene layer plane exists in the X-Y plane and the thickness direction of Z is relatively small, axial modulus of SLGSs can be averaged as: E = (E_x + E_y)/2 = 1025.81 GPa, with comparison of 996.49 GPa by Berryman [45]. Homogenized moduli of TLGSs based on AFEM can be calculated, and elastic moduli of polycrystalline graphene are given in Table 4.

Parameters	C11	C12	C13	C33	C44	C66	Bv	Br	Gv	Gr	E	μ
FEM	1057.21	183.49	15.00	36.49	4.86	436.84	286.43	35.75	218.48	11.00	507.27	0.212
AFEM	1123.02	201.33	5.50	35.72	2.56	460.82	300.71	34.39	231.15	6.07	525.01	0.214
Equation-AFEM	1123.10	201.39	5.12	38.24	2.61	460.86	300.86	36.65	231.40	6.20	526.39	0.215
Equation-FEM	1057.51	183.83	13.29	45.24	4.86	436.83	286.78	43.64	219.30	11.19	511.53	0.217
Graphene [45]	1060.00	180.00	15.00	36.50	4.5	440.00	286.28	35.76	219.57	10.26	507.86	0.212

Table 4.

Elastic constants of grapheme and homogenized moduli (bulk, shear, and Young?s moduli) of polycrystalline grapheme by FEM and AFEM (unit: GPa).

From references of Berryman [46, 47], elastic moduli for polycrystalline graphenes can be calculated based on the elastic constant matrix C_ij of graphene. Elastic constants and elastic modulus by AFEM are in agreement with results of Berryman [48]. Stress-strain and elastic moduli of TLGSs by AFEM are shown in Figure 14.

Figure 14.

Comparison of TLGSs stress-strain and elastic moduli by AFEM modeling. (a) Stress-strain and axial modulus by AFEM. (b) Comparison of Young’s and shear moduli.

From Figure 14(a), elastic constants are in quite good agreement; axial modulus in Z direction is much smaller than that of X, Y direction. From Figure 14(b), elastic moduli by AFEM are nevertheless slightly lower than values by Berryman. The nonspherical nature in Young’s modulus of TLGSs is anisotropic and the x-axis is more anisotropic than the z-axis because of the values of C₃₃ larger than C₁₁. Elastic moduli of TLGSs by Wu et al. [38] are as: B = 167.55 GPa, G = 118.21 GPa. Thus, Young’s modulus of TLGSs is averaged to be 0.525 TPa, which is lower than SLGSs; this trend is consistent with the trend between the bilayer graphene sheets (BLGSs) of 0.8 TPa [49] and SLGSs of 0.996 TPa [45], showing van der Waals forces between layers has weakened the Young’s modulus value of multigraphene.

5. Application to the hexagonal portlandite (CH) structure

CH is a typical nanostructure of the microscale-hydrated Portland cement paste [50] and its Young’s modulus is needed in a 3D multiscale mortar model [51], so we aim to propose a modeling method based on AFEM to investigate anisotropic stiffness of CH crystalline and elastic properties of polycrystals [52].

5.1 Nanoscale modeling and coupling force field of CH

5.1.1 Nanoscale modeling of CH ionic crystal

CH is a hexagonal structure; its crystal lattice [53] is: a = b = 3.5930 Å, c = 4.9090 Å, α = β = 90°, γ = 120°. AFEM modeling of portlandite (CH) is shown in Figure 15.

Figure 15.

AFEM modeling of portlandite (CH). (a) Portlandite crystal. (b) CH model with beam-spring element.

From Figure 15(a), the Ca²⁺ is octahedrally coordinated by oxygen. There is no hydrogen bond across the layer. The direction of hydroxyl groups formed by oxygen and hydrogen is vertical to the (001) plane direction. The atomic distance between O²⁻ and H⁺ is 2.757 Å. O▬H covalent bond can be solved by a linear bond energy ratio relative to C▬C bond energy. For O▬H covalent bond, L = 0.985 nm [54]. Portlandite is modeled with equivalent beam-spring element in Figure 15(b). The potential cutoff radius is 0.8 nm.

Parameters of weight, charges, and Lennard-Jones potential are given in Table 5.

Weight and charges [55]			Lennard-Jones potential
Symbol	Weight [u]	Charges [e]	Bond	Dij (eV)	Rij (nm)
Oh	16.00	−1.40	O▬O	0.15540	0.35532
Ho	1.00	0.40	O▬Ca	3.83368 × 10^−5	0.48980
Ca	39.90	2.00	Ca▬Ca	2.18105 × 10^−7	0.62428
Ohw	16.00	−0.95	O▬H	6.73854 × 10^−3	0.35532
How	1.00	0.42	H▬H	0.02000	0.20000

Table 5.

CH parameters of atomic charges and atomic force field employed (1e = 1.6021892 × 10⁻¹⁹ C).

5.1.2 Coupling force field definition of CH ionic crystal

For CH modeling, the initial conditions are as: (1) no hydrogen bond or other strong bonds across the layer; (2) without considering the effect of three-body potential and dihedral; (3) O▬H is approximately defined by estimation of C▬C bond according to the linear bond energy ratio. Moreover, CH parameters by AFEM modeling are as: O▬H bond length L_O▬H = 0.986 nm, potential cutoff radius is 0.8 nm, E_O▬H = 4.611 × 10³ nN/nm², d_O▬H = 0.175 nm, G_O▬H = 1.921 × 10³ nN/nm². Spring element to describe the interatomic forces under the equilibrium state is shown in Figure 16.

Figure 16.

Ionic bond’s definition and various atomic types of CH by AFEM modeling. (a) VDW force interaction. (b) Coulomb force of interaction. (c) Coupling force field.

From Figure 16(a), the curve of VDW force has been translated so that the system at the beginning of the calculation is at the equilibrium (no more attractive/repulsive force). From Figure 16(b), Coulomb force curve is given, then the coupling force field of Ca▬O and O▬H bond is decayed rapidly with increasing distance, and Ca▬Ca and O▬O bonds are weakened during the modeling of CH, as shown in Figure 16(c).

5.2 Elastic constants and homogenized moduli of CH

Similarly, elastic constants of CH crystal are verified by volume unit. Then, we extend to the equivalent beam-spring model by AFEM. Strain loading of CH crystal is shown in Figure 17.

Figure 17.

Strain loading of CH crystal. (a) Strain loading in x-direction. (b) Strain loading in y-direction. (c) Strain loading in z-direction.

Figure 17 shows the stretching schematic diagram of CH crystal in three directions. The CH crystal modeling based on AFEM is extensively concerned in this section. CH crystal modeling using the equivalent beam-spring element is mainly investigated, and elastic modulus of CH hexagonal polycrystals can be finally determined based on the calculated elastic constants matrix.

More particularly, the results of the AFEM are compared to previous results of experiment [56]. Comparison of CH stress-strain and elastic modulus by AFEM modeling is shown in Figure 18.

Figure 18.

Stress-strain and comparison of elastic moduli by various methods. (a) Stress-strain and axial modulus by AFEM. (b) Comparison of Young’s and shear moduli. (c) Comparison of Young’s modulus of single CH crystal. (d) Comparison of shear modulus of single CH crystal.

From Figure 18(a), it can be seen that, even if all the results are in quite good agreement, the numerical results using AFEM are slightly lower than those obtained by Brillouin spectroscopy measurement. Anisotropy of elastic modulus in each direction with various angles of z-axis is shown in Figure 18(b). The comparison of elastic modulus curves based on C_ij shows anisotropy in each direction with various angles of z-axis in Figure 18(c). As can be seen from Figure 18(d) and axial modulus formulas, even if all the results are in quite good agreement, it can be noted that DFT results are upper than those of Brillouin spectroscopy measurement, whereas AFEM ones are lower.

Elastic constant C_ij is then calculated and axial modulus can be obtained by Eqs. (45) and (47). Axial moduli based on elastic constants by AFEM are as: σ_x = σ_y = 86.4918ε_x, σ_z = 31.7002ε_z. From Figure 18, shear and Young’s moduli by AFEM are slightly lower than Brillouin spectroscopy measurement. Moreover, elastic moduli by Wu et al. [38] are as: B = 28.249 GPa, G = 17.314 GPa. Thus, Young’s modulus of CH polycrystals is averaged to be 43.130 GPa by Reuss-Voigt-Hill estimation, which is in quite good agreement with experimental averaged value (39.88 GPa) and with the literature values (45.94 GPa by Laugesen, 52.4 GPa by Speziale, 44.69 GPa by Kerisit, and 46.58 GPa by Holuj). AFEM has a certain meaning to provide a method to calculate elastic modulus of other hexagonal structures.

Strain-stress data and elastic moduli (bulk, shear, and Young’s moduli) by AFEM are given in Tables 6 and 7. Strain-stress data of CH crystal by AFEM are given in Table 6. Elastic modulus of polycrystalline CH compared with previous results of the literature is given in Table 7.

Strain-stress by volume unit anisotropy				Strain-stress by beam-spring element
σij (nN/nm^2)	ε11	ε22	ε33	σij (nN/nm^2)	ε11	ε22	ε33
σ11 = 74.3832	0.834493	−0.249651	−0.118612	σ11 = 41.0121	0.453152	−0.064053	−0.052899
σ22 = 59.7722	−0.200612	0.670576	−0.095314	σ22 = 32.3671	−0.057148	0.337453	−0.052226
σ33 = 36.1151	−0.057590	−0.057590	1.018542	σ33 = 11.3929	−0.041024	−0.041024	0.3483412
σ12 = 4.0644	γ12 = 0.515788			σ12 = 3.5601	γ12 = 0.501512

Table 6.

Strain-stress data of CH by volume unit using FEM and beam-spring element using AFEM modeling.

Method	C11	C12	C13	C33	C44	C66	Bv	Br	Gv	Gr	μ
FEM	99.38	30.78	7.36	36.29	7.88	34.31	36.23	26.63	22.65	13.92	0.2
AFEM	94.02	25.51	6.02	32.31	7.099	34.25	32.83	23.67	21.88	12.75	0.246
Equation-AFEM	95.806	25.294	4.626	35.671	7.098	35.256	32.93	24.59	22.74	12.99	0.243
Equation-FEM	100.22	31.611	6.389	37.430	7.878	34.305	36.29	26.79	22.91	14.00	0.255
Brillouin [56]	102.00	32.00	8.40	33.60	12.00	35.00	37.24	26.02	24.39	18.18	0.225

Table 7.

Elastic constants of CH and homogenized moduli (bulk, shear, and Young?s moduli) of polycrystalline CH by FEM and AFEM (unit: GPa).

6. Conclusions

Equivalent beam elements are applied in zigzag nanotubes of C▬C covalent bond, as well as the AFEM modeling of a tensile test on CH and TLGSs. These models at nanoscale are established and then elastic constants are determined based on several methods/algorithm existed and developed. Conclusions are as:

1. A theoretical formula between axial modulus E_x,_y,_z and c_ij of hexagonal crystal is derived by crystal elastic theory. Besides, a reliable equation to calculate these C_ij is also given out and verified by both volume unit (FEM) and equivalent beam-nonlinear spring elements (AFEM), which are as follows:

   1. The expressions of the x-axial modulus E₁, y-axial modulus E₂, and z-axial modulus E₃ of hexagonal crystal as a function of elastic constants C_ij are as:

      $$E_3=c_{33}-\frac{2c_{13}^2}{c_{11}+c_{12}};E_1=E_2=\left(c_{11}-c_{12}\right)\left[1+\left(c_{12}{c}_{33}-c_{13}^2\right)/\left(c_{11}{c}_{33}-c_{13}^2\right)\right]$$

   2. The equation of elastic constants Cij calculated by the stress components and strain components using the volume unit (FEM method) are derived and detailed as in Eq. (43), where $A={\epsilon }_{11},B={\epsilon }_{22},C={\epsilon }_{33},D={\epsilon }_{31},E={\epsilon }_{11}+{\epsilon }_{21}$ (${\epsilon }_{\mathit{ij}}$ are independent strains in all directions).

2. The product E.t of Young’s modulus and wall thickness for CNT is mainly distributed in the constant range of 0.33–0.37 kJ/m². Young’s modulus of zigzag SWCNT is 1.05 Tpa and E.t is close to be 0.3604 kJ/m². Besides, the calculated E of TLGSs model is close to the measurement of Berryman.

3. Young’s modulus of single-layer graphene calculated by AFEM is 1025.81 GPa, compared with 996.49 GPa by Berryman, and the relative errors of c_ij are in an acceptable range with experimental values.

4. CH crystal has shown anisotropy. Based on the stress-strain data of AFEM modeling, elastic constant c_ij is calculated, Young’s modulus is estimated to be 43.130 GPa using Reuss-Voigt-Hill estimation, which is close to the Brillouin spectroscopy result of 52.15 GPa.

As an attempt, the computational efficiency of AFEM method of CH crystal on ABAQUS software is greatly improved. Besides, the computation cycle (no more than 5 min in total) is greatly reduced. Thus, based on the content above, AFEM can be considered as a kind of alternative method to obtain mechanical properties of crystals at nanoscale.

However, it should be noted that molecular mechanics methods cannot provide information of the electronic structure as they omit the electronic nature of the atoms. That is why the first-principles simulation (the density functional theory, etc.) is also introduced to calculate the Young’s moduli of other typical structures considered in this thesis, under the situation of lacking a complete potential function parameters or more complex chemical bonds existed in the structures.

In short, we can see that AFEM modeling methodology at the atomic scale (or nanoscale) can also be an effective method for the specific models, since the macroscopic properties depend largely on the physicochemical properties of the interatomic bonds. Multiscale methods, including the AFEM method, therefore connect nanoscale modeling and continuous pattern of deformation behavior by identifying relevant parameters from small scales to larger scales, thus have their prospective applications. It is not far for us to deeply investigate mechanical properties with the rapid development of computer hardware level, the interatomic potential functions accurately described, new high-efficiency algorithms, etc. Perspectives of the computer tools, accurate interatomic potential functions, and new high-efficiency algorithms will complete this meaningful work and study the mechanical behavior of materials more deeply.

Acknowledgements

The author acknowledges the financial support provided by the China Scholarship Council and the startup foundation of Xi’an Shiyou University. Thanks Qiufeng Wang for her proofreading, and very sincere thankness to prof. Fabrice Bernard and prof. Siham Kamali-Bernard for the their enthusiastic and responsible guideness.

Classical Density Functional Theory (DFT): Elastic Constants of Typical Single Crystals

Abstract

Anisotropy becomes the focus of condensed matter physics, materials physics, and materials chemistry research, where molecular pressure characteristics are critical to better understand the mechanical properties of the materials at the atomic scale. So the influencing effect of pressure on structural, elastic properties of single cubic crystals (CaO and MgO), single hexagonal crystals (CH and calcite), and single monoclinic crystals (11 Å tobermorite and gypsum) are mainly investigated. This chapter is committed to model and calculate elastic constants of typical anisotropic crystal structures. Above all, the functional local density approximation (LDA) or generalized gradient approximation (GGA) is mainly used. Then, the elastic constants are calculated, which can be used as the measure criterion of the resistance of a crystal to an externally applied stress. Besides, the optimized structural lattice parameters at zero pressure are calculated to compare with the experimental parameters, which all have good agreement with the experimental and theoretical values. Therefore, it is fairly meaningful to study the elastic constants to understand the physical, chemical, and mechanical properties of these cubic, hexagonal, and monoclinic crystal structures. The results show that the applied pressure is beneficial to the elastic properties. These researches provide basic physical parameters for their homogenized moduli of these typical polycrystals.

Keywords: DFT, lattice parameters, elastic constants, GGA, LDA, pressure effect

1. Introduction

Anisotropy can effectively increase the freedom of material, thus can be used to provide controllability and design new materials, and anisotropy becomes the focus of condensed matter physics, materials physics, and materials chemistry research. Molecular pressure characteristics are critical to better understand the mechanical properties of the materials at the atomic scale.

First-principles simulation starting from quantum mechanics is a powerful tool for materials research, of which the density functional theory (DFT) is commonly used to study the crystal structure, lattice energy, the equation of state, the electronic bandgap, and vibration spectral properties [1]. DFT is confirmed to be one of the important nanoscale modeling methods. With the increasing interest in the advanced applications of cementitious materials at the nanoscale, elastic constants of single crystal are essential to the material performance because many other parameters related to mechanical properties can be derived from them. This resistant force of bond in crystal can be characterized by energy gap, and the number of bond per unit area can be determined by valence electron density.

DFT has a great importance of simplifying the solving process of the complex Schrödinger equation. The internal degrees of freedom of the nuclei extend to a scale of several of magnitude orders (smaller than the electron variables and concentrated at the bulk of the mass), which causes a significant slowdown compared to electrons. Based on the kinetic energy density functional of Thomas [2] and the exchange-correlation effects of Dirac [3], the density functional theory has been greatly developed by Kohn and Sham [4], who have established the fundamental approximation theorem on the functional status to describe real systems by electronic structure calculations. Assuming that the validity of the exchange-correlation functional and the accuracy of density and energy convergence are correct, we should also clear that the eigenvalues of KS equations have no physical meaning and the ionization energy is in the opposite state direction [5].

Many studies show that the band structure obtained by KS equations remains more reasonable prediction of experimental energy spectra in the main usage of functional local density approximation (LDA) or generalized gradient approximation (GGA). Moreover, the use of the potential and the disturbance of eigenfunctions allow us to determine the excitation energies during KS solving process. One proposed approach is to introduce the eigenstates to calculate multibody (many-body calculation) on the basis of Monte Carlo calculations [6] and perturbation theory [7]. The approach of Kohn and Sham (KS) [4] proposed in 1965 is still the most commonly used approach in DFT calculations, which aims to determine the specific properties of a multiparticle system using independent particle and owns satisfactory success.

So, we mainly investigate some typical crystals (such as: cubic CaO and MgO, hexagonal portlandite and calcite, monoclinic 11 Å tobermorite and gypsum), focusing on the elastic constants within a certain pressure range, so as to lay a foundation for the elastic modulus calculations.

2. Molecular modeling using DFT by materials studio

2.1 Introduction of density functional theory (DFT)

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 ${\sigma }_i$ is positioned $\overrightarrow{R}=\left\{{\overrightarrow{R}}_i;i=1,\dots N_n\right\}$ and $\overrightarrow{r}=\left\{\left({\overrightarrow{r}}_i,{\sigma }_i\right);i=1,\dots N_e\right\}$, respectively.

$$\mathit{H\psi }\left(\overrightarrow{R},\overrightarrow{r}\right)=\mathit{E\psi }\left(\overrightarrow{R},\overrightarrow{r}\right)\hspace{1em}\mathrm{(1)}$$

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

A multiparticle electronic structure satisfies the Schrödinger equation. Generally, Hamiltonian contains kinetic energy, Coulomb potential energy between electrons and the external potential. 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)\hspace{1em}\mathrm{(2)}$$

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 material elastic constants and Debye temperature with the accurate predictions. Therefore, DFT-QHA method is quite suitable to study physical properties of crystal structure under high pressure.

The Hamiltonian and Kohn-Sham equations as an approximation to simplification of Schrodinger equation to solve the multiparticle systems are [8]:

$$\widehat{H}=-\frac{\mathrm{\hslash }}{2m_e}\sum _i{\nabla }_i^2+\sum _i{V}_{\mathit{ex}}\left(r_i\right)+\frac{1}{2}\sum _{i\ne j}\frac{e^2}{\left|r_i-r_j\right|}\hspace{1em}\mathrm{(3)}$$

$$\left\{\begin{array}{l}\left[-\frac{\mathrm{\hslash }}{2m}{\nabla }^2+V_{\mathit{eff}}\left(\overrightarrow{r}\right)\right]{\phi }_i\left(\overrightarrow{r}\right)={\epsilon }_i{\phi }_i\left(\overrightarrow{r}\right)\\ V_{\mathit{eff}}\left(r\right)=V\left(\overrightarrow{r}\right)+V_{\mathit{ee}}\left(\overrightarrow{r}\right)+V_{\mathit{ex}}\left(\overrightarrow{r}\right)\\ V_{\mathit{xe}}\left(\overrightarrow{r}\right)=\frac{\delta E_{\mathit{xe}}\left[\rho \right]}{\delta \left(\overrightarrow{r}\right)}\\ \rho \left(\overrightarrow{r}\right)=\sum _i{\left|{\varphi }_i\left(\overrightarrow{r}\right)\right|}^2\end{array}\right.\hspace{1em}\mathrm{(4)}$$

where V_ex(r) is an external potential; r_i and r_j are the nucleus position vector; while m_ε, m, respectively, stand for the quality of nuclei and electrons; ρ(r) is the electron density. Eqs. (3) and (4) include not only the nucleus and electron kinetic energy term, but also interaction term of nucleus-nucleus and electron-electron [8].

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 is 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}}}{c}_{\mathit{\text{ijkl}}\dots }^T{\eta }_{\mathit{ij}}{\eta }_{\mathit{kl}}\dots \hspace{1em}\mathrm{(5)}$$

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:

$$\mathrm{\Delta }E=\frac{V}{2}\sum _{i=1}^6\sum _{j=1}^6{c}_{\mathit{ij}}{e}_i{e}_j\hspace{1em}\mathrm{(6)}$$

As the independent elastic constants of different crystals are different from each other, through the different matrix elements of the crystal elastic modulus, we can get different strain applied methods.

2.2 Introduction of the DFT and CASTEP module

At present, there are a lot of software for molecular modeling, such as: Ceruis2, InsightII, VASP, Gaussian, Material Studio, etc., of which Ceruis2, InsightII run under Unix mainframe systems and thus are inconvenient to use. Material Studio can run a set of quantum mechanics, molecular mechanics, mesoscopic model, and statistical correlation simulation. Moreover, Materials Studio has an environment of server/client mode, which has brought the advanced material simulation and modeling techniques.

Materials Studio CASTEP has been applied to a wide range of research problems such as surface chemistry, and chemisorption, heterogeneous catalysis, defects in semiconductors, grain boundaries, stacking faults, nanotechnology, molecular crystals, polymorphic studies, diffusion mechanisms, and molecular dynamics of liquids. Materials Studio CASTEP uses a total energy plane-wave pseudopotential method. In the mathematical model of the material, Materials Studio CASTEP replaces core electrons with effective potentials acting only on the valence electrons in the system. Electronic wavefunction are expanded through a plane-wave basis set, and exchange and correlation effects can be included within either the local density approximation (LDA) or generalized gradient approximation (GGA). Combining the use of pseudopotentials and plane wave basis sets enables extremely efficient geometry optimizations of molecules, solids, surfaces, and interfaces. Based on total energy pseudopotential methods, Materials Studio CASTEP can predict properties such as lattice constants, molecular geometry, elastic constants, band structures, density-of-states, charge densities and wave functions, and optical properties. The pseudopotential plane-wave technology underlying CASTEP is well validated, and efficient parallel versions of the code are also available for large systems involving hundreds of atoms.

CASTEP is capable of computing many electronic and optical properties using density functional perturbation theory (DFPT), known as the linear response method. The primary reason that CASTEP has become so powerful is that the numerical methods used to solve the underlying quantum mechanical calculations are both computationally efficient and extremely accurate. This approach makes possible a wider variety of properties than are possible using the so-called finite difference approaches, which require repeated computations on series systems.

For a detailed description of the quality settings, we can set the calculation parameters as in Table 1.

Value	Coarse	Medium	Fine	Ultrafine
SCF tolerance (eV/atom)	1.0 × 10^−5	2.0 × 10^−6	1.0 × 10^−6	5.0 × 10^−7
k-Points separation (1/Å)	0.07	0.05	0.04	0.04
Energy cutoff	Coarse	Medium	340 eV	380 eV
Energy tolerance (eV/atom)	5.0 × 10^−5	2.0 × 10^−5	1.0 × 10^−5	5.0 × 10^−6
Max. force tolerance (eV/Å)	0.1	0.05	0.03	0.01
Max. stress tolerance (GPa)	0.2	0.1	0.05	0.02
Max. displacement tolerance (Å)	5.0 × 10^−3	2.0 × 10^−3	1.0 × 10^−3	5.0 × 10^−4

Table 1.

Parameter accuracy of CASTEP module.

Convergence criteria of the accuracy of the calculation are controlled in Table 1 by the quality settings. Specify the number of steps for each strain and maximum strain amplitude in the appropriate text boxes. The minimizer tab of the CASTEP geometry optimization dialog is used to obtain the optimized structure. Specify the quality using either the dropdown list or by providing individual convergence tolerances for the energy, maximum force, and maximum displacement in the corresponding text boxes.

Elastic constants will be generated for the theoretical lattice constants after cell optimization. Theoretical calculation of Kohn-Sham equations is shown in Figure 1(a), and the corresponding DFT calculation process is shown in Figure 1(b).

Figure 1.

Solution of Kohn-Sham equations. (a) Theoretical calculation of Kohn-Sham equations. (b) Flowchart of solving Kohn-Sham equations by DFT.

From Figure 1, the calculation of elastic constants is proceeded by full geometry optimization, including cell optimization. Elastic constants are evaluated by calculating the stress tensor for a number of distorted structures. Internal coordinates are optimized in each run, while keeping the lattice parameters fixed. The accuracy of the elastic constants depends to a great extent on the accuracy of the self-consistent field (SCF) part and also on the level of convergence of geometry optimizations for each distorted structure.

3. Theoretical approximate solution of DFT calculation

The crystalline ion movement of the electron is as: $\psi \left(\overrightarrow{R},\overrightarrow{r}\right)=\chi \left(\overrightarrow{R}\right)\phi \left(\overrightarrow{R},\overrightarrow{r}\right)$ and assume that the electron mobility ($\phi$) does not depend on the speed nuclei but on their positions.

3.1 Solving equation of adiabatic approximation

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

$$\widehat{H}={\widehat{T}}_e+{\widehat{U}}_{\mathit{ne}}+{\widehat{U}}_{\mathit{ee}}+{\widehat{U}}_{\mathit{nn}}\hspace{1em}\mathrm{(7)}$$

$${\widehat{H}}_{\overrightarrow{R}}{\phi }_{\overrightarrow{R}}^0\left(\overrightarrow{r}\right)=E^{\mathit{BO}}\left(\overrightarrow{R}\right){\phi }_{\overrightarrow{R}}^0\left(\overrightarrow{r}\right)\left(\overrightarrow{r}\right)\hspace{1em}\mathrm{(8)}$$

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:

$$E^{\mathit{BO}}\left(\overrightarrow{R}\right)⟨{\phi }_{\overrightarrow{R}}^0\left|\widehat{H}\right|{\phi }_{\overrightarrow{R}}^0⟩=min⟨{\phi }_{\overrightarrow{R}}^0\left|\widehat{H}\right|{\phi }_{\overrightarrow{R}}^0⟩\hspace{1em}\mathrm{(9)}$$

This energy has a surface in the space coordinates that is said ionic Born Oppenheimer surface.

The ions move according to the effective potential energy, including Coulomb repulsion and the anchoring effect of the electron, as:

$${\widehat{H}}^{\mathit{BO}}={\widehat{T}}_n+E^{\mathit{BO}}\left(\overrightarrow{R}\right)\hspace{1em}\mathrm{(10)}$$

$${\widehat{H}}^{\mathit{BO}}\chi \left(\overrightarrow{R}\right)=\mathit{E\chi }\left(\overrightarrow{R}\right)\hspace{1em}\mathrm{(11)}$$

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 formulation applies to any system of mutually interacting particles in an external potential $V_{\mathit{ex}}\left(\overrightarrow{r}\right)$, where the Hamiltonian is written as:

$$\widehat{H}=-\frac{\mathrm{\hslash }}{2m_e}\sum _i{\nabla }_i^2+\sum _i{V}_{\mathit{ex}}\left(r_i\right)+\frac{1}{2}\sum _{i\ne j}\frac{e^2}{\left|r_i-r_j\right|}\hspace{1em}\mathrm{(12)}$$

DFT and its founding principle is summarized in two theorems, first introduced by Hohenberg and Kohn, which say there is bijection between the set of potential $V_{\mathit{ex}}\left({\overrightarrow{r}}_i\right)$ and the density minimizing the Eq. (9), based on the following points:

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}}_{\mathit{ex}}\left(r\right)+E_{\mathit{nm}}\left(\overrightarrow{R}\right)\cdot F_{\mathit{HK}}\left[n\right]+\int d^3{\mathit{rV}}_{\mathit{ex}}\left(r\right)+E_{\mathit{nn}}\left(\overrightarrow{R}\right)\hspace{1em}\mathrm{(13)}$$

As a result, the density n₀(r) minimizing the energy associated with the Hamiltonian in Eqs.(13) 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^{\mathit{BO}}\left[\overrightarrow{R}\right]=minE\left(\overrightarrow{R},n\left(\overrightarrow{r}\right)\right)\hspace{1em}\mathrm{(14)}$$

Because the ground state is concerned, it is possible to replace (in the Hilbert space) 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 $\int n\left(\overrightarrow{r}\right)d^{3}r=N_e$. In this stage, the DFT can reformulate the problem rather than solve an uncertain functional $F_{\mathit{HK}}\left(n\right)$.

3.2 The approximation approach of Kohn-Sham

This approach assumes that the density in the ground state of the system is equal to that in some systems composed of noninteracting particles. This involves independent particle equations for the noninteracting (numerically manageable) system, gathering all the terms complicated and difficult to assess, in a functional exchange-correlation $E_{\mathit{XC}}\left(n\right)$.

$$E_{\mathit{KS}}=F\left[n\right]+\int d^3{\mathit{rV}}_{\mathit{ex}}\left(r\right)=T_S\left[n\right]+E_H\left[n\right]+E_{\mathit{XC}}\left[n\right]+\int d^3{\mathit{rV}}_{\mathit{ex}}\left(r\right)\hspace{1em}\mathrm{(15)}$$

T is the kinetic energy of a system of particles (electrons) independent (noninteracting) embedded in an effective potential which is no other than the real system,

$$T_S\left[n\right]=⟨{\psi }_{\mathit{NI}}\left|{\widehat{T}}_e\right|{\psi }_{\mathit{NI}}⟩=\sum _{i=1}^{N_e}⟨{\phi }_i\left|-\frac{1}{2}{\nabla }^2\right|{\phi }_i⟩\hspace{1em}\mathrm{(16)}$$

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

$$E_{\mathit{\text{Hartree}}}\left[n\right]=\frac{1}{2}\int d^3{\mathit{rd}}^3{r}^{\prime }\frac{n\left(r\right)n\left(r^{\prime }\right)}{\left|r-r^{\prime }\right|}\hspace{1em}\mathrm{(17)}$$

$$n\left(r\right)=\sum _{i=1}^{N_e}{\left|{\phi }_i\left(r\right)\right|}^2\hspace{1em}\mathrm{(18)}$$

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 that T_S is orbital function, all other terms depend on the density. Therefore, it is possible to vary the functions of the wave and deriving 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\hspace{1em}\mathrm{(19)}$$

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

$$\left({\widehat{H}}_{\mathit{KS}}-{\epsilon }_i\right){\phi }_i\left(r\right)=0\hspace{1em}\mathrm{(20)}$$

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

$${\widehat{H}}_{\mathit{KS}}\left(r\right)=-\frac{1}{2}{\nabla }^2+V_{\mathit{KS}}\left(r\right)\hspace{1em}\mathrm{(21)}$$

$$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)}\hspace{1em}\mathrm{(22)}$$

Eqs. (20)–(22) are known equations of Kohn Sham, the density n(r), and the resulting total energy $E_{\mathit{KS}}$. These equations are independent of any approximation on the functional $E_{\mathit{XC}}\left(n\right)$, resolution provides the exact values of the density and the energy of the ground state of the interacting system, provided that $E_{\mathit{XC}}\left(n\right)$ is known. The latter can be described in terms of Hohenberg Kohn function (12)

$$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)\hspace{1em}\mathrm{(23)}$$

$$E_{\mathit{xc}}\left[n\right]=⟨\widehat{T}⟩-T_S\left[n\right]+⟨{\widehat{V}}_{\mathrm{int}}⟩-E_{\mathit{\text{Hartree}}}\left[n\right]\hspace{1em}\mathrm{(24)}$$

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 (LDA and GGA) and the parameter set of Bloch theorem and Brillouin zone, the content is detailed in our previous chapter [10].

4. Elastic coefficients of typical crystals solved by DFT method

4.1 Elastic strain energy formulas of cubic crystals

The matrix element of cubic crystal system is $\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{12}& 0& 0& 0\\ c_{12}& c_{11}& c_{12}& 0& 0& 0\\ c_{12}& c_{12}& c_{11}& 0& 0& 0\\ 0& 0& 0& c_{44}& 0& 0\\ 0& 0& 0& 0& c_{44}& 0\\ 0& 0& 0& 0& 0& c_{44}\end{array}\right]$, where the subscript of matrix C corresponds to i and j. In the calculation of the elastic modulus, the strain component needs to be applied to the same $\delta$.

Strain tensor exhibits a total of six independent components. When a strain is applied, the corresponding strain component e_i (i = 1~6) can be expressed as:

$$e=\left(e_1,e_2,e_3,e_4,e_5,e_6\right)\hspace{1em}\mathrm{(25)}$$

Submitting the cubic crystal matrix into Eqs. (6), we can obtain the following formula:

$$\begin{array}{c}\text{Δ}E=\frac{V}{2}(c_{11}{e}_1{e}_1+c_{11}{e}_2{e}_2+c_{11}{e}_3{e}_3+c_{12}{e}_1{e}_2+c_{12}{e}_1{e}_3+c_{12}{e}_2{e}_1+c_{12}{e}_2{e}_3+c_{12}{e}_3{e}_1\\ +c_{12}{e}_3{e}_2+c_{44}{e}_4{e}_4+c_{44}{e}_5{e}_5+c_{44}{e}_6{e}_6)\end{array}(26)$$

Firstly, considering the simplest constant C₄₄, i and j are equal to 4, and by the cubic matrix elements, C₄₄ = C₅₅ = C₆₆. By setting $e=\left(0,0,0,\delta ,\delta ,\delta \right)$ of Eq. (25) and submitting into Eq. (5), we get:

$$\mathrm{\Delta }E=\frac{V_0}{2}\left(c_{44}{e}_4{e}_4+c_{44}{e}_5{e}_5+c_{44}{e}_6{e}_6\right)\hspace{1em}\mathrm{(27)}$$

Therefore,

$$\frac{\mathrm{\Delta }E}{V_0}=\frac{3}{2}{c}_{44}{\delta }^2\hspace{1em}\mathrm{(28)}$$

Similarly, setting $e=\left(\delta ,\delta ,0,0,0,0\right)$, we obtain:

$$\mathrm{\Delta }E=\frac{V_0}{2}\left(c_{11}{e}_1{e}_1+c_{11}{e}_2{e}_2+c_{12}{e}_1{e}_2+c_{12}{e}_2{e}_1\right)\hspace{1em}\mathrm{(29)}$$

Therefore,

$$\frac{\mathrm{\Delta }E}{V_0}=\left(c_{11}+c_{12}\right){\delta }^2\hspace{1em}\mathrm{(30)}$$

On this basis, if a strain $e=\left(\delta ,\delta ,\delta ,0,0,0\right)$ is imposed, we have:

$$\mathrm{\Delta }E=\frac{V_0}{2}\left(c_{11}{e}_1{e}_1+c_{11}{e}_2{e}_2+c_{11}{e}_3{e}_3+c_{12}{e}_1{e}_2+c_{12}{e}_1{e}_3+c_{12}{e}_2{e}_1+c_{12}{e}_2{e}_3+c_{12}{e}_3{e}_1+c_{12}{e}_3{e}_2\right)\hspace{1em}\mathrm{(31)}$$

Therefore,

$$\frac{\mathrm{\Delta }E}{V_0}=\frac{3}{2}\left(c_{11}+2c_{12}\right){\delta }^2\hspace{1em}\mathrm{(32)}$$

After fitting the three equations, we can obtain three independent elastic constants c₁₁, c₁₂, and c₄₄.

A transverse strain causes a change in shape without a change in volume. Thus, by applying these five specific strains mentioned above and selecting a series of $\delta$ strain amplitude, we can get the data from the ΔE-δ curve, respectively, and then obtain by fitting quadratic coefficient according to the corresponding applied strain methods mentioned in Table 2. Finally, we can get the independent elastic constants of cubic crystals by these simultaneous equations.

Deformation tensor	LCEC
$e=\left(\delta ,\delta ,0,0,0,0\right)$	2c11 + 2c12
$e=\left(\delta ,-\delta ,0,0,0,0\right)$	2c11−2c12
$e=\left(0,0,0,0,0,\delta \right)$	3c11 + 6c12
$e=\left(\delta ,\delta ,\delta ,0,0,0\right)$	3c11 + 6c12
$e=\left({\delta }^2/\left(1-{\delta }^2\right),0,0,2\delta ,0,0\right)$	4c44
$e=\left(0,0,0,\delta ,\delta ,\delta \right)$	3c44

Table 2.

Deformation tensors to calculate independent elastic constants of cubic crystal.

For the cubic phase, the criteria for mechanical stability are given by [11]:

$$C_{11}>0,C_{44}>0,C_{11}>\left|C_{12}\right|,\left(C_{11}+2C_{12}\right)>0\hspace{1em}\mathrm{(33)}$$

4.2 Elastic strain energy formulas of hexagonal crystal structures

Matrix element of hexagonal crystal structure is as: $\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& 0& 0& 0\\ c_{12}& c_{11}& c_{13}& 0& 0& 0\\ c_{13}& c_{13}& c_{33}& 0& 0& 0\\ 0& 0& 0& c_{44}& 0& 0\\ 0& 0& 0& 0& c_{44}& 0\\ 0& 0& 0& 0& 0& \frac{1}{2}\left(c_{11}-c_{12}\right)\end{array}\right]$, system has a total of five independent components, of which independent elastic constants are c₁₁, c₁₂, c₁₃, c₃₃, and c₄₄.

Similarly, we can get the strain tensor expansion by strain energy, thus leading to the relationship between strain energy and second-order elastic coefficients of the hexagonal crystal:

$$\begin{array}{c}{\rho }_{0}F\left({\eta }_J\text{,}T\right)={\rho }_{0}F\left({\eta }_J\text{,}T\right)+\frac{1}{2}\{c_{11}^T{\eta }_1{\eta }_1+c_{11}^T{\eta }_2{\eta }_2+c_{33}^T{\eta }_3{\eta }_3+c_{44}^T{\eta }_4{\eta }_4+c_{44}^T{\eta }_5{\eta }_5\\ +\frac{c_{11}^T-c_{12}^T}{2}{\eta }_6{\eta }_6+2c_{12}^T{\eta }_1{\eta }_2+2c_{13}^T{\eta }_1{\eta }_3+2c_{13}^T{\eta }_2{\eta }_3\}\begin{array}{l}+\dots +\frac{1}{n!}{\displaystyle \sum }_{JP\dots }{c}_{JP\dots }^T{\eta }_J{\eta }_{P\cdots }\hfill \end{array}\end{array}(34)$$

For hexagonal crystal, the basis vector of original cell can be taken as:

$$\overrightarrow{R}=\left(\begin{array}{c}{a}_1\\ a_2\\ a_3\end{array}\right)=\left(\begin{array}{ccc}\frac{\sqrt{3}}{2}a& \frac{1}{2}a& 0\\ -\frac{\sqrt{3}}{2}a& \frac{1}{2}a& 0\\ 0& 0& c\end{array}\right)\hspace{1em}\mathrm{(35)}$$

where a and c are the lattice constants of crystals.

We can apply strain $e=\left(\delta ,\delta ,0,0,0,0\right)$ to calculate the value of c₁₁ + c₁₂:

$$\frac{\mathrm{\Delta }E}{V_0}=\left(c_{11}+c_{12}\right){\delta }^2\hspace{1em}\mathrm{(36)}$$

Similarly, we can apply strain $e=\left(0,0,0,0,0,\delta \right)$ to calculate the value of c₁₁−c₁₂:

$$\frac{\mathrm{\Delta }E}{V_0}=\frac{1}{4}\left(c_{11}-c_{12}\right){\delta }^2\hspace{1em}\mathrm{(37)}$$

Also for c₃₃, we can apply strain $e=\left(0,0,\delta ,0,0,0\right)$:

$$\frac{\mathrm{\Delta }E}{V_0}=\frac{1}{2}{c}_{33}{\delta }^2\hspace{1em}\mathrm{(38)}$$

Then, for c₄₄, we can apply strain $e=\left(0,0,0,\delta ,\delta ,0\right)$:

$$\frac{\mathrm{\Delta }E}{V_0}=c_{44}{\delta }^2\hspace{1em}\mathrm{(39)}$$

Finally, we can apply strain $e=\left(\delta ,\delta ,\delta ,0,0,0\right)$ to calculate the values of c₁₁, c₁₂, c₁₃, and c_33:

$$\frac{\mathrm{\Delta }E}{V_0}=\left(c_{11}+c_{12}+2c_{13}+\frac{1}{2}{c}_{33}\right){\delta }^2\hspace{1em}\mathrm{(40)}$$

For hexagonal structure, mechanical stability condition is as [10]:

$$c_{44}>0,c_{11}>\left|c_{12}\right|,\left(c_{11}+2c_{12}\right)c_{33}>2c_{13}^2\hspace{1em}\mathrm{(41)}$$

Similarly, by applying these five specific strains mentioned above and selecting a series of $\delta$ strain amplitude, the independent elastic constants of hexagonal crystals by these simultaneous equations of ΔE-δ data in Table 3 can be obtained.

Distortion matrix	Energy (Taylor’s expansion)	$\frac{1}{V_0}\frac{{\partial }^{2}E(V,\mathit{\delta })}{\partial {\delta }^2}{|}_{\delta }$
$\left(\begin{array}{ccc}1+\delta & 0& 0\\ 0& 1+\delta & 0\\ 0& 0& 1\end{array}\right)$	$\begin{array}{l}E\left(V,\delta \right)=E\left(V_0,0\right)\\ +V_0\left[\left({\tau }_1+{\tau }_2\right)\delta +\left(c_{11}+c_{12}\right){\delta }^2\right]\end{array}$	$2\left(c_{11}+c_{12}\right)$
$\left(\begin{array}{ccc}1+\delta & 0& 0\\ 0& 1-\delta & 0\\ 0& 0& 1\end{array}\right)$	$\begin{array}{l}E\left(V,\delta \right)=E\left(V_0,0\right)\\ +V_0\left[\left({\tau }_1-{\tau }_2\right)\delta +\left(c_{11}-c_{12}\right){\delta }^2\right]\end{array}$	$2\left(c_{11}-c_{12}\right)$
$\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1+\delta \end{array}\right)$	$\begin{array}{l}E\left(V,\delta \right)=E\left(V_0,0\right)\\ +V_0\left[{\tau }_3\delta +1/2c_{33}{\delta }^2\right]\end{array}$	$c_{33}$
$\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& \delta \\ 0& \delta & 1\end{array}\right)$	$\begin{array}{l}E\left(V,\delta \right)=E\left(V_0,0\right)\\ +V_0\left[{\tau }_5\delta +2c_{44}{\delta }^2\right]\end{array}$	$4c_{44}$
$\left(\begin{array}{ccc}1+\delta & 0& 0\\ 0& 1& 0\\ 0& 0& 1+\delta \end{array}\right)$	$\begin{array}{l}E\left(V,\delta \right)=E\left(V_0,0\right)\\ +V_0\left[\left({\tau }_1+{\tau }_3\right)\delta +1/2\left(c_{11}+2c_{13}+c_{33}\right){\delta }^2\right]\end{array}$	$c_{11}+2c_{13}+c_{33}$
$\left(\begin{array}{ccc}1+\delta & 0& 0\\ 0& 1& 0\\ 0& 0& 1-\delta \end{array}\right)$	$\begin{array}{l}E\left(V,\delta \right)=E\left(V_0,0\right)\\ +V_0\left[\left({\tau }_1+{\tau }_3\right)\delta +1/2\left(c_{11}-2c_{13}+c_{33}\right){\delta }^2\right]\end{array}$	$c_{11}-2c_{13}+c_{33}$
$\left(\begin{array}{ccc}1+\delta & 0& 0\\ 0& 1+\delta & 0\\ 0& 0& 1+\delta \end{array}\right)$	$\begin{array}{l}E\left(V,\delta \right)=E\left(V_0,0\right)\\ +V_0\left[\left({\tau }_1+{\tau }_2+{\tau }_3\right)\delta + 1/2\left(2c_{11}+2c_{12}+4c_{13}+c_{33}\right){\delta }^2\right]\end{array}$	$2c_{11}+2c_{12}+c_{33}+4c_{13}$

Table 3.

Deformation tensor to determine independent elastic constants of hexagonal crystal [12, 13].

The corresponding applied strain methods are mentioned in Table 2. Elastic strain energy and strain relations for hexagonal structure can be given in Table 3.

Last but not the least, for elastic strain energy formulas of monoclinic crystals, the content is detailed in our previous chapter work [9].

5. Elastic coefficients of typical single crystals at nanoscale

5.1 Elastic coefficients C_ij of cubic crystals: CaO and MgO

All DFT calculations were performed by the total energy of the linear combination of atomic orbitals with a Hartree-Fock self-consistent field [14, 15] method implementation. Monkhorst Pack using the Brillouin zone integration method [16] is used to set K point to verify the convergence of different sampling meshes.

5.1.1 Modeling of cubic crystals: CaO and MgO

5.1.1.1 Lime (CaO) cubic structure

CaO is an alkaline earth oxide of native lime. In fact, there are many studies of CaO phase transition characteristics using quantum mechanics. The high-pressure structure phases and the elastic properties of CaO are important for mortar modeling in civil engineering. CaO has the B₁–B₂ transition at 61 GPa [17], and phase transformation from B₁ to B₂ takes over a range of pressure 60–70 GPa [18] from shockwave measurements. Kalpana et al. [19] and Baltache et al. [20] also investigated this transition for CaO. Mehl et al. [21] predicted elastic constants and their pressure dependences of CaO using potential-induced breathing (PIB) models. The ab initio pseudopotential calculations with plane-wave basis (PWPP) and ab initio full-potential linear muffin-tin-orbitals (FP-LMTO) are used to study high-pressure elastic properties of CaO by Karki and Crain [22] and Tsuchiya and Kawamura [23], respectively. Tsuchiya et al. [23] have used the full-potential linearized muffin-tin-orbital method to study elastic constants of MgO and CaO.

Lime (CaO) is a cubic structure, the initial lime crystal lattice by Fiquet et al. [24] is as: a = b = c = 4.8330 Å, α = β = γ = 90°, with space group Fm3m.

As is shown in Figure 2, Ca and O are arranged alternately to occupy cube vertex, edge, and face-center. From the lattice unit terms, the ion at the apex, on the edge, on the surface, and inside of cell, each ion belonging to the cell separately have the occupancy rate of 1/8, 1/4, 1/2, and 1. Atomic coordinates and displacement parameters of lime (CaO) are given in Table 4.

Figure 2.

Crystal structure of CaO. a) CaO cubic crystal b) crystal in x-direction c) crystal in y-direction d) crystal in z-direction.

Atom	x	y	z	Occupancy rate	Uiso or Ueq
Ca	0.000	0.000	0.000	1.00	1.00
O	0.500	0.500	0.500	1.00	1.00

Table 4.

Atomic coordinates and displacement parameters of lime (CaO) [24].

The elastic properties of the B₁-phase CaO are calculated by using the plane-wave pseudopotential density functional theory (DFT) method. The Cambridge Serial Total Energy Package (CASTEP) program is used to obtain the calculated energy-volume points. The lattice constants, the bulk modulus, and the pressure derivatives of bulk modulus of the B₁-CaO at zero pressure are given. In addition, the elastic constants, shear modulus, and bulk modulus of the B₁-phase CaO are investigated over a wide pressure range from 0 to 1.0 GPa.

5.1.1.2 Periclase (MgO) cubic structure

In 1995, Duffy et al. [25] have found, for pressures up to 227 GPa, that B1-type MgO remained stable, while at 300 K, the stable pressure can reach up to 199 GPa. By Anvil hydrostatic test, Fei [26] has measured and got high-pressure state equation as well as elastic modulus of MgO. By ultrasonic interferometry method, Chen et al. [27] have found that the elastic anisotropy of MgO structure decreases with an increase of pressure under normal temperature conditions, while the elastic anisotropy increases with temperature at high pressure. Jaekson and Niesler [28] have obtained isothermal bulk modulus of MgO and its first derivative of pressure. Sinogeikin and Bass [29] have studied the pressure state equation of MgO by Anvil Brillouin scattering experiments. Zha et al. [30] have experimentally measured elastic properties with the pressure up to 55 GPa by scattering spectroscopy X Brillouin zone. Merkel et al. [31] have studied the pressure up to 47 GPa and measured elastic modulus, shear modulus, and elastic anisotropy. Wolf and Bukowinski [32] have investigated the stabilization of the ionic charge densities in MgO and CaO crystals by an electron gas model. Matsui et al. [33] have calculated elastic constants and thermodynamic properties of MgO structure over a wide temperature and pressure range. Karki et al. [34], using the local density approximation of density functional perturbation theory, have studied the MgO thermal expansion coefficient and elastic properties. Karki et al. [35] have used a plane-wave qualitative study of the elastic properties of the potential B1-type MgO in the range of 0–25 GPa.

Periclase (MgO) is a cubic structure, the initial lime crystal lattice by Kalpana [19] is as: a = b = c = 4.2130 Å, α = β = γ = 90°, with space group Fm3m.

As is shown in Figure 3, Mg and O are arranged alternately to occupy cube vertex, edge, and face-center. Mg²⁺ ions will fill into the O⁻² ion accumulation to form MgO crystal.

Figure 3.

Crystal structure of MgO. (a) MgO cubic crystal. (b) Crystal in x-direction. (c) Crystal in y-direction. (d) Crystal in z-direction.

Atomic coordinates and displacement parameters of MgO are given in Table 5.

Atom	x	y	z	Occupancy rate	Uiso or Ueq
Mg	0.000	0.000	0.000	1.00	1.00
O	0.500	0.500	0.500	1.00	1.00

Table 5.

Atomic coordinates and displacement parameters of MgO [19].

Electronic basis is optimized in calculation using a linear combination of atomic orbitals by self-consistent field (SCF) method of Hartree Fock. The K points of convergence have been met when taking a 16 × 16 × 16 mesh in our calculations. For complex space, we use parameter of Gilat and Raubenheimer [36]: when the total energy reaches less than 10⁻⁶ eV/atom, the self-consistent convergence is more effective.

In addition, the elastic constants, Young’s modulus, and shear modulus are investigated over a wide pressure range from 0 to 1.0 GPa.

5.1.2 Initial conditions and elastic constants of CaO and MgO

5.1.2.1 Initial conditions and elastic constants of CaO

All the total energy electronic structure calculations are implemented through the CASTEP code. In CaO crystal modeling, we adopt the exchange-correlation interaction by using the GGA and LDA of the exchange-correlation function. Ultrasoft pseudopotentials with cutoff energy of the plane waves of 500 eV are used to describe the electron-ion interaction. For the Brillouin zone sampling, the 17 × 17 × 17 of Monkhorst-Pack mesh [15] is used, where the self-consistent convergence of the total energy is at 10⁻⁶ eV per atom.

For CaO crystal at zero temperature and pressure equilibrium, lattice constant a and bulk modulus B₀ are as shown in Table 6.

	a	C11	C12	C44	B	G
Present (LDA)	4.7124	289.497	71.432	76.652	144.12	88.295
Present (GGA)	4.8339	200.879	56.902	75.600	104.89	74.134
Ref. [21]	4.72	274	54	53	128	—
Ref. [22]	4.820	206	50	66	102	71
Ref. [23]	4.838	239	51.6	77.4	117	83.6
Ref. [24]	4.840	223	53	84	109	85
Exp. [37]		221.89	57.81	80.32		80.59
Exp. [38]		220.53	57.67	80.03		80.59

Table 6.

Lattice a, elastic constants, bulk and shear moduli of CaO at zero pressure.

The structural and the elastic properties of CaO are investigated by ab initio plane-wave pseudopotential density functional theory method. Lattice a, elastic constants, bulk and shear moduli of CaO under different pressures are shown in Table 7.

P	a	V(Å^3)	c11	c12	c44
0.0 (LDA)	4.7124	104.647	289.497	71.432	76.652
0.0 (GGA)	4.8339	112.951	200.879	56.902	75.600
0.1 (GGA)	4.8324	112.847	201.548	56.742	75.621
0.2 (GGA)	4.8310	112.748	202.538	56.986	75.639
0.3 (GGA)	4.8299	112.665	203.062	56.893	75.659
0.4 (GGA)	4.8279	112.534	204.506	57.361	75.683
0.5 (GGA)	4.8264	112.428	205.667	57.736	75.717
0.6 (GGA)	4.8251	112.336	206.105	57.440	75.717
0.7 (GGA)	4.8238	112.249	207.001	57.664	75.753
0.8 (GGA)	4.8219	112.110	208.202	57.762	75.784
0.9 (GGA)	4.8204	112.008	209.415	58.192	75.890
1.0 (GGA)	4.8193	111.934	209.830	58.059	75.823

Table 7.

Lattice a, elastic constants, bulk and shear moduli of CaO under different pressures.

As can be seen from the Table 8, the obtained a, C_ij, B, and G of the B₁-phase CaO agree well with the experimental data [37, 38] and other calculated results [21, 22, 23, 24]. CaO cubic crystal under pressure 0–1.0 GPa by GGA is shown in Figure 4.

	a	C11	C12	C44	B	G	A
Present (LDA)	4.2130	298.508	92.228	141.197	160.988	124.506	—
Present (GGA)	4.2130	343.468	106.358	140.823	185.395	131.449	—
Ref. [19]	4.2130	—	—	—	167.600	—	—
Ref. [22]	4.30	270	73	127	139	115	0.21
Ref. [24]	4.2590	318	87	144	164.000	132	0.18
Ref. [29]	—	298.96	96.42	157.13	163.930	131.80	—
Ref. [33]	4.2210	300	93.6	147	162.4	127.6	0.29
Ref. [35]	4.2510	291	91	139	158.00	121.800	0.27
Ref. [39]	4.16–4.24	286–352	91–108	158–188	157.0–198.0	—	—
Ref. [40]	4.212	—	—	—	—	—	—

Table 8.

Lattice a, elastic constants C_ij, bulk and shear moduli of MgO at zero pressure.

Figure 4.

CaO cubic crystal under pressure 0–1.0 GPa by DFT. (a) Relative change of V. (b) Elastic constants.

From Figure 4, we can see the pressure dependences of the elastic constants c_ij of CaO at zero temperature and different pressures. It is found that c₁₁ varies largely under the effect of pressure, compared with the variations in c₁₂ and c₄₄. The change of elastic constant c₁₁ represents a produced change with a longitudinal strain in length elasticity. The elastic constants c₁₂ and c₄₄ are related to the elasticity in shape, which is a shear constant. Besides, c₁₂ and c₄₄ are less sensitive of pressure as compared with c₁₁. As pressure increases, c₁₁ and c₁₂ increase monotonically, but c₄₄ slightly decreases monotonically.

5.1.2.2 Initial conditions and elastic constants of MgO

Similarly, the exchange-correlation interaction by using the GGA and LDA of the exchange-correlation functions is adopted. Ultrasoft pseudopotentials with cutoff energy of the plane waves are set to 700 eV to describe the electron-ion interaction. For Brillouin zone, the 15 × 15 × 15 of Monkhorst-Pack mesh [15] is used, where the self-consistent convergence of the total energy is at 10⁻⁶ eV per atom.

As can be seen in Table 8, the obtained a, C_ij, B, and G of the MgO agree well with the experimental data [29] and other calculated results [19, 22, 24, 33, 35, 39, 40]. It is noted that c₁₂ and c₄₄ for MgO are positive indicating that the MgO is stable at zero pressure (Table 9).

P	a	V (Å^3)	c11	c12	c44
0.0 (GGA)	4.2130	74.778	343.468	106.358	140.823
0.0 (LDA)	4.2130	74.778	298.508	92.228	141.197
0.1 (GGA)	4.2038	74.289	309.892	95.878	142.181
0.2 (GGA)	4.2030	74.246	308.989	94.200	142.178
0.3 (GGA)	4.2022	74.202	310.994	95.635	142.300
0.4 (GGA)	4.2012	74.153	312.885	96.635	142.371
0.5 (GGA)	4.2005	74.115	313.050	96.187	142.484
0.6 (GGA)	4.1996	74.069	313.565	96.092	142.519
0.7 (GGA)	4.1988	74.025	314.766	96.493	142.597
0.8 (GGA)	4.1979	73.980	314.612	95.599	142.677
0.9 (GGA)	4.1971	73.937	316.353	96.652	142.857
1.0 (GGA)	4.1963	73.893	317.067	96.738	142.940

Table 9.

Lattice a, elastic constants, bulk and shear moduli of MgO under different pressures.

In the case of cubic crystals, there are only three independent elastic constants of c₁₁, c₁₂, and c₄₄. MgO cubic crystal under pressure 0–1.0 GPa by GGA is shown in Figure 5.

Figure 5.

MgO cubic crystal under pressure 0–1.0 GPa by DFT. (a) Relative change of V. (b) Elastic constants.

From Figure 5, we can see that c₁₁ varies largely under the effect of pressure as compared with the variations in c₁₂ and c₄₄.

5.2 Elastic coefficients C_ij of hexagonal crystals: CH and calcite

Portlandite is one of the main typical constituents of hydrated cementitious systems. CH represents 17–25% of the volume fraction of portland cement paste, and its Young’s modulus is needed in the modeling of cement systems on the macro- and microscales [41, 42]. Portlandite influences the physical and mechanical properties as well as the durability of cement-based materials. Calcite has very large industrial applications in the field of civil engineering construction, and it can also be used as a flux in glass and in the field of metallurgy.

5.2.1 Modeling of hexagonal crystals: CH and calcite

5.2.1.1 Portlandite hexagonal structure

Laugesen [43] has calculated the elastic properties of portlandite using the density functional calculations, where elastic constants are given by GGA method. Speziale et al. [44] have determined the elastic constants of portlandite by Brillouin spectroscopy. Besides, Busing et al. [45] have investigated the calcium hydroxide structure by the neutron diffraction, of which the atomic information is given and thus is used to establish the CH model. The morphology of the portlandite crystal is hexagonal when the structure of the crystal is trigonal [45] with space group P3m1. Initial model of CH is as shown in Figure 6.

Figure 6.

Final optimized parameters of CH hexagonal crystal. (a) CH hexagonal crystal. (b) Real CH cell. (c) Cell in x-direction. (d) Cell in y-direction. (e) Cell in z-direction.

As shown in Figure 6, Ca²⁺ is octahedrally coordinated by oxygen. CH is a hexagonal structure, the initial crystal lattice [43] is: a = b = 3.5930 Å, c = 4.9090 Å, α = β = 90°, γ = 120°. There is no hydrogen bond or other strong bonds across the layer. Ca²⁺ is octahedrally coordinated by oxygen and each O²⁻ is tetrahedrally coordinated by calcium and hydrogen. The direction of hydroxyl groups formed by oxygen and hydrogen is vertical to the (001) plane direction. The atomic distance between O²⁻ and H⁺ is 2.757 Å (Table 10).

Atomic species	X	Y	Z	Occupancy rate	Uiso or Ueq
Ca	0.0000	0.0000	0.0000	1.00	1.00
O	0.3333	0.6667	0.2346	1.00	1.00
H	0.3333	0.6667	0.4280	1.00	1.00

Table 10.

Atomic coordinates and displacement parameters of portlandite [45].

In order to determine whether the elastic coefficients under various pressures are stable, the pressure region of 0–5.0 GPa is added, with the same initial conditions of Laugesen [43].

5.2.1.2 Calcite hexagonal structure

Several researchers [46, 47] have been interested in elastic constants of calcite and its structure. Each carbon atom of calcite is the center of an equilateral triangle, and each vertex of the triangle is occupied by the oxygen atom. The calcite is crystallizing in the trigonal system “or rhombohedral” system. Zhang and Reeder [48] have investigated compressibilities of calcite-structure carbonates and confirmed the results found by Aydinol et al. [49] on the post- and pyroxene-aragonite phase type. Rohl et al. [50] have developed a new force field for calcite model. Effenberger et al. [51] discussed some aspects of the stereochemistry of calcite-type carbonates. Pilati et al. [52] have given the atomic displacement parameters in the calcite structure. Zaoui and Sekkal [53] have investigated the mechanisms behind the ikaite-to-calcite phase transformation by using molecular dynamics calculations. Le Page et al. [54] have investigated calcite using ab initio calculation to study their stiffness properties. Meanwhile, the elastic properties of calcite structure have been studied by Chen et al. [55] using Brillouin spectroscopy. Dandekar and Ruoff [56] have investigated the MgO structure on the temperature dependence of the elastic constants. Prencipe et al. [57] have studied its structure and properties using an ab initio quantum-mechanical calculation.

Initial hexagonal model of calcite by Prencipe [57] is as: a = b = 5.0492 Å, c = 17.343 Å, α = β = 90°, γ = 120°, with space group R-3C. The calcite model is as shown in Figure 7.

Figure 7.

Calcite hexagonal model. (a) Calcite hexagonal model. (b) In x-direction. (c) In y-direction. (d) In z-direction. (e) CO₃ group.

The initial lattice of calcite hexagonal crystal is as shown in Figure 7(a). From Figure 7(a) and (b), calcite crystal has the following features: (1) two Ca atoms are connected to the anionic groups (CO₃⁻²), form together an equilateral triangle, including the central carbon atom and each vertex totalizing three oxygen atoms. (2) CO₃ group interatomic distances are as: C▬O is 1.2815 Å, Ca▬O₁ is 3.4589 Å, O▬O₁ is 2.2197 Å, O▬O₂ is 4.7911 Å, O▬O₃ is 3.4107 Å, and O▬O₄ is 4.2459 Å (Table 11).

Atomic species	X	Y	Z	Occupancy rate	Uiso or Ueq
Ca	0.0000	0.0000	0.0000	1.00	1.00
C	0.0000	0.0000	0.2500	1.00	1.00
O	0.25642	0.0000	0.2500	1.00	1.00

Table 11.

Atomic coordinates and displacement parameters of calcite [57].

In order to determine whether the elastic coefficients under various pressures are stable, the pressure region of 0–5.0 GPa using DFT methods is added.

5.2.2 Initial conditions and elastic constants of CH and calcite

5.2.2.1 Initial conditions and elastic constants of CH

Initial conditions are as: GGA of the exchange-correlation function is used. Ultrasoft pseudopotentials with cutoff energy of the plane waves are 700 eV. Brillouin zone is 12 × 12 × 12. The self-consistent convergence of the total energy is at 10⁻⁶ eV per atom. The simulation of CH hexagonal crystal under pressure 0–5.0 GPa by GGA method is shown in Figure 8.

Figure 8.

CH hexagonal crystal under pressure 0–5.0 GPa by DFT. (a) Relative change of a, c, and V. (b) Elastic constants.

The changes of elastic constants (GPa) of CH under various pressures are separately calculated, and the results are shown in Table 12.

	Pressure	a0(Å)	c0(Å)	V0 (Å^3)	C11	C12	C13	C14	C33	C44	C66
Cal. [43]	0	3.5930	4.9090	54.8830	99.39	30.78	7.36	0	36.29	7.88	34.305
Exp. [44]	0	3.5927 (±0.0009)	4.908 (±0.005)	54.87 (±0.06)	102.0 (±2.0)	32.0 (±1.0)	8.4 (±0.4)	4.5 (±0.2)	33.6 (±0.7)	12.0 (±0.3)	35.0 (±1.1)
Cal. [58]	0	3.6580	4.7930	55.5000	135.1	45.6	6.6	0	26.0	4.1	44.7
Exp. [59]	0	3.5890	4.9110	54.8000	99.28 (±1.5)	36.21 (±2.0)	29.65 (±15.0)	≈0 (assumed)	32.6 (±2.0)	9.864 (±2.0)	31.55 (±1.5)
Present	0.0 (GGA)	3.5918	4.9063	54.8162	103.948	30.923	2.564	4.336	42.670	14.140	36.513
	0.0 (LDA)	3.6233	5.0086	56.9445	95.417	27.162	1.313	0.242	27.585	8.998	34.127
	0.5 (GGA)	3.6118	4.9383	55.7899	97.140	27.042	9.553	0.503	35.899	12.003	35.049
	1.0 (GGA)	3.5999	4.86716	54.6274	101.463	33.613	10.193	2.725	54.182	10.904	33.925
	1.5 (GGA)	3.5903	4.8207	53.8137	107.685	33.889	12.593	3.396	50.867	12.045	36.898
	2.0 (GGA)	3.5815	4.7855	53.1590	107.236	34.681	17.650	0.285	49.072	20.818	36.278
	2.5 (GGA)	3.5733	4.755	52.5737	112.249	37.190	17.139	2.661	58.666	19.959	37.529
	3.0 (GGA)	3.5651	4.7245	52.0042	114.388	39.664	22.111	2.372	64.582	16.340	37.362
	3.5 (GGA)	3.5572	4.6961	51.4622	117.505	41.661	26.485	2.846	65.960	16.476	37.922
	4.0 (GGA)	3.5492	4.6718	50.9649	119.696	44.869	26.498	0.648	73.285	17.543	37.414
	4.5 (GGA)	3.5415	4.6461	50.4651	120.400	44.261	25.089	2.851	73.884	22.543	38.070
	5.0 (GGA)	3.5342	4.6258	50.0389	128.865	49.445	36.702	4.781	83.203	20.346	39.710

Table 12.

Elastic constants (GPa) of CH under various pressures.

From Figure 8, for CH crystal at 0 GPa, elastic coefficients under initial conditions by DFT are as: c₁₁ = 103.948 GPa, c₁₂ = 30.923 GPa, c₁₃ = 2.564 GPa, c₃₃ = 42.670 GPa, and c₄₄ = 14.140 GPa. C₁₄ is turned out not to be zero, which is different to the result of Laugesen. We can see that all elastic coefficients (c₁₄ not included) under various pressures are stable except for a pressure of 1.0 and 4.5 GPa.

5.2.2.2 Initial conditions and elastic constants of calcite

Generalized gradient approximation (GGA) of the exchange-correlation function is used. Ultrasoft pseudopotentials with cutoff energy of the plane waves are 700 eV. Brillouin zone is 12 × 12 × 12. The self-consistent convergence of the total energy is 10⁻⁶ eV per atom.

Results of hexagonal calcite crystal under pressure 0–5.0 GPa by GGA is shown in Figure 9.

Figure 9.

Calcite hexagonal crystal under pressure 0–5.0 GPa by DFT. (a) Relative change of a, c, and V. (b) Elastic constants.

From Figure 9, for calcite crystal, elastic coefficients under initial conditions at 0 GPa by DFT are given as follows: c₁₁ = 143.938 GPa, c₁₂ = 55.579 GPa, c₁₃ = 50.970 GPa, c₃₃ = 81.175 GPa, and c₄₄ = 32.658 GPa. The elastic coefficients of calcite under pressure 0–5 GPa show slight growing trends except for c₁₁ and c₁₂ within a pressure region of 2.5–3.5 GPa.

The calculated elastic constants of calcite are given in Table 13.

P (GPa)	a0 (Å)	c0 (Å)	c0/a0	V0 (Å^3)	c11	c12	c13	c14	c33	c44	c66
Ref. [46]	—	—	—	—	146.26	59.70	50.76	−20.76	85.31	34.05	43.28
Ref. [48]	4.989	17.044	—	367.39	—	—	—	—	—	—	—
Ref. [49]	4.981	15.902	—	341.676	—	—	—	—	—	—	—
Ref. [51]	4.990	17.061	—	367.8	—	—	—	—	—	—	—
Ref. [54]	4.988	17.068	—	367.761	159	70	63	19.5	98	39	44.5
Ref. [55]	—	—	—	—	149.4	57.9	53.5	20	85.2	34.1	45.75
Ref. [56]	—	—	—	—	145.7	55.9	53.5	20.5	85.3	33.4	44.9
Ref. [8]	5.061	17.097	—	379.279	152.3	57.05	54.83	17.14	87.77	36.00	47.625
0.0 (LDA)	4.9718	16.5296	3.3247	353.526	153.9261	55.8064	50.9698	13.6541	82.7015	33.9451	59.9905
0.0 (GGA)	5.0548	17.2551	3.4136	381.818	143.9383	55.5791	50.9698	16.8678	81.1753	32.6583	44.1796
0.5 (GGA)	5.0478	17.1808	3.4036	379.127	148.479	57.7574	53.4967	17.6855	83.2565	33.6887	45.3608
1.0 (GGA)	5.0411	17.1089	3.3939	376.536	150.1629	60.2232	55.5284	17.3573	84.8005	34.5971	44.9699
1.5 (GGA)	5.0357	17.0338	3.3826	374.074	152.0584	60.9377	58.008	17.8412	87.0412	35.0652	45.5604
2.0 (GGA)	5.0293	16.9642	3.3731	371.609	159.5886	63.1081	60.7533	17.8817	88.253	34.6942	48.2403
2.5 (GGA)	5.0245	16.8939	3.3623	369.350	158.2239	75.896	61.9659	19.5687	90.0313	35.5319	41.1640
3.0 (GGA)	5.0194	16.8242	3.3518	367.090	164.1449	68.8673	65.1085	21.3005	91.4266	36.7669	47.6388
3.5 (GGA)	5.0144	16.7569	3.3418	364.885	167.8984	72.1265	67.4523	21.4184	92.9625	37.7036	47.8860
4.0 (GGA)	5.0101	16.6887	3.3310	362.789	170.1009	74.8428	69.618	21.4534	94.683	37.3224	47.6291
4.5 (GGA)	5.0049	16.6285	3.3224	360.724	172.8895	76.9844	72.4475	22.8814	96.4539	37.6871	47.9526
5.0 (GGA)	5.0009	16.5620	3.3118	358.704	177.1157	78.6186	74.8454	23.6571	97.6551	37.548	49.2486

Table 13.

The calculated elastic constants of calcite (unit: GPa).

From Table 13, elastic constants have been calculated by GGA method as follows: c₁₁ = 143.938 GPa, c₁₂ = 55.579 GPa, c₁₃ = 50.970 Pa, c₃₃ = 81.175 GPa, and c₄₄ = 32.658 GPa. By LDA method, elastic constants have been calculated as follows: c₁₁ = 153.926 GPa, c₁₂ = 55.806 GPa, c₁₃ = 50.970 Pa, c₃₃ = 82.702 GPa, and c₄₄ = 59.991 GPa.

Similarly, for elastic constants C_ij of monoclinic crystals (gypsum and 11 Å tobermorite), the content is detailed in our previous chapter work [9]. Furthermore, elastic moduli of cubic, hexagonal, and monoclinic structures at larger scale can be homogenized by the calculated elastic constants mentioned, as in Chapter 4.

6. Conclusions

Elastic constants of typical crystal structures under a certain pressure region are calculated by DFT method. Three types of the cubic, hexagonal, and monoclinic crystals are considered, which have a certain value for both application and reference. Results are as:

1. For cubic structure, results are as:

   1. Elastic constants of CaO single crystal are as: c₁₁ = 200.879 GPa, c₁₂ = 56.902 GPa, and c₄₄ = 75.600 GPa.

   2. Elastic constants of our MgO single crystal are as: c₁₁ = 298.508 GPa, c₁₂ = 92.228 GPa, and c₄₄ = 141.197 GPa.

2. For hexagonal structure, results are as:

   1. For CH crystal, elastic constants under initial conditions at 0 GPa based on GGA method are as: c₁₁ = 103.948 GPa, c₁₂ = 30.923 GPa, c₁₃ = 2.564 GPa, c₁₄ = 4.336 GPa, c₃₃ = 42.670 GPa, and c₄₄ = 14.140 GPa. Besides, elastic constants based on LDA method are as: c₁₁ = 95.417 GPa, c₁₂ = 27.162 GPa, c₁₃ = 1.313 GPa, c₁₄ = 0.242 GPa, c₃₃ = 27.585 GPa, and c₄₄ = 8.998 GPa. We can see that C₁₄ is not zero, which is different to the result of Laugesen. The elastic coefficients under various pressures are stable except for a pressure of 0.8 GPa.

   2. For calcite crystal, elastic coefficients under initial conditions at 0 GPa by DFT are given as follows: c₁₁ = 143.938 GPa, c₁₂ = 55.579 GPa, c₁₃ = 50.970 GPa, c₃₃ = 81.175 GPa, and c₄₄ = 32.658 GPa. The elastic coefficients of calcite under pressure 0–5 GPa show slight growing trends except for a pressure of 2.5–3.5 GPa.

Moreover, for monoclinic single-crystal structure (11 Å tobermorite and gypsum), elastic coefficients are detailed in our previous chapter work of “Elastic constants and homogenized moduli of monoclinic structures based on density functional theory.”

Acknowledgements

The author acknowledges the financial support provided by the China Scholarship Council and the startup foundation of Xi’an Shiyou University. Thanks Qiufeng Wang for her proofreading, and very sincere thankness to prof. Fabrice Bernard and prof. Siham Kamali-Bernard for the their enthusiastic and responsible guideness.

Homogenized Moduli of Typical Polycrystals by Elastic Constants Based on DFT

Abstract

X-ray method to test the material properties and to obtain elastic constants is commonly based on the Reuss and Kroner models. Y parameter has been turned out to be an effective method to estimate elastic properties of polycrystalline material. Since Y parameters of cubic polycrystalline material based on the certain uniform stress (Reuss model) have not been given, our work aims to complete this part of the theoretical analysis which can effectively compare elastic constants measured by the X-ray diffraction method. The structural and the elastic properties of cubic structures (CaO and MgO) and hexagonal structures (CH and Calcite CaCO₃) are investigated by the density functional theory method, and then, the credibility of Y parameters for determining elastic moduli of cubic structures is proved and elastic properties in typical crystallographic planes of [100], [110], and [111] are also calculated. Meanwhile, Young’s moduli of CaO, MgO, CH, and calcite structure are 175.76, 293.17, 58.08, 84.549, 54.30, and 79.51 GPa, which are all close to references. Elastic properties of cubic and hexagonal structures under various pressures are calculated and the surface constructions of elastic moduli are drawn, showing the anisotropy at various directions. Then, the bulk modulus, shear modulus, Young’s modulus, Poisson’s ratio, and anisotropy factor are determined.

Keywords: anisotropic elasticity, Young’s modulus, Y parameter, cubic crystal, hexagonal crystal, polycrystalline materials

1. Introduction

The development of computational methods to describe and predict the mechanical properties of materials is of obvious practical importance. The accuracy of the prediction of mechanical or, from a more general point of view, physical properties depends largely on the knowledge of the intrinsic properties of the various constitutive phases. Up to now, these properties have been solely obtained by experimental procedures. The tremendous increase of computational capabilities has largely favored the development of numerical modeling based on a realistic and multiscale description of these kinds of materials. Since several years, some of the authors of this chapter have been involved in the development of such numerical tools from the micro- to the macrolevel (see, for instance, [1] or [2]). The intent is to assist the material developer by providing a rational approach to material development and concurrently assist the structural designer by providing an integrated analysis tool that incorporates fundamental material behavior. It is now clear that any realistic attempts to accomplish such an objective should be based on appropriate nanostructure/performance relationships, and also on an efficient method to upscale the properties determined on the lowest level to the next higher one. This is the context of the work detailed in this chapter.

Generally, many industrial and engineering materials consisted of dissimilar constituents are essentially inhomogeneous, which are distinguishable at a small length scale. Each constituent composed by anisotropic crystals shows different material properties and/or material orientations at some smaller length scales. At a sufficiently fine scale, all materials can be considered to be heterogeneous. The crystal anisotropy would be very difficult to handle computationally when the crystals are overlapping, the structure of isotropic polycrystals as the real materials at large scale are usually isotropic. The length-scale problem is considered as a great importance factor in material mechanics. Each material has structural features from the microstructure when studied at sufficiently physical parameters, from the view of homogenization. Polycrystalline structure contains a variety of information (e.g., orientation) and the properties of a single crystal, such as anisotropy. Within the mechanics of typical crystal structures, the transition from the micro- to the mesoscale (homogenization) and vice versa (localization) can be estimated.

Homogenization is an idealized description of a statistical distribution inside the actual heterogeneous material. At the engineering level, there exists such a scale (intermediate between the microscopic scale—that of the constituents—and the scale of the structure); it is said that the material can be homogenized. As opposed to the scale of constituents, the macroscopic homogeneity is called to be of statistical homogeneity or microscopic. The concept of homogenization is gradually proposed, the Gologanu-Leblond-Devaux’s (GLD) analysis of an rigid-ideal plastic (von Mises) spheroidal volume gradually is extended to the case when the matrix is anisotropic (obeying Hill’s anisotropic yield criterion), and the representative volume element (RVE) is subjected to arbitrary deformation. For quasibrittle materials, Zhu et al. [3] have formulated the anisotropic model in the framework of Eshelby-based homogenization methods. However, for the homogenization of elastic deformation, especially for polycrystalline structures, it is still the traditional Reuss-Voigt-Hill method. The homogenized moduli can be determined using the homogenization theory where periodic media are implemented. The hypothesis of continuity of a material is introduced to overcome theoretical difficulty, which implies a notion of statistical average, and the actual constitution of the material is idealized by considering the material to be continuous. Once the continuity model is admitted, the concept of homogeneity is deduced from it [4]. A homogeneous medium is then characterized by properties that are identical at every point. This hypothesis is particularly suitable for polycrystalline model.

1.1 Limits of X-ray diffraction methods and Voigt-Reuss-Hill averaging

Macroscopic averaging assuming a random distribution of initial crystals leads to isotropic mechanical properties such as Young’s modulus. The upscaling to higher levels is then done by means of the Y-parameter theory which has been introduced by Zheng et al. [5, 6] a few years ago to overcome the limits of X-ray diffraction methods and conventional Voigt-Reuss-Hill averaging, as it is recalled in the following of the introduction, but which can also be seen as an efficient way to estimate homogenized properties.

X-ray method is commonly used to test the physical properties of materials. By satisfying the Bragg crystal diffraction conditions, the values, the shapes, and the position changes of the diffraction peak are tested to calculate the distance between the crystal surfaces. Then, the strain value at the normal direction of the crystal surface can be solved through the distance between the crystal surfaces mentioned above. Normally, elastic constants measured by X-ray method contain the function of diffraction surface, which can be averaged within the total space to determine the macroscopic elastic properties of the material [7]. By multiplying the relevant X-ray elastic constant and the strain value, the stress in the material can be thus obtained. Until now, for the commonly used X-ray stress analyzer, the principle of X-ray diffraction to test the material stress and to obtain elastic constants [8] is based on the uniform certain stress Reuss model. However, for the uniform certain strain Voigt model, there are still some difficulties in explaining its principle by X-ray diffraction method. Based on the content mentioned above, Y parameter is introduced on the basis of theoretical formulas’ derivation to explain both Reuss and Voigt models.

However, for the uniform certain strain Voigt model, there are still some difficulties in explaining its principle by X-ray diffraction method. X-ray method is associated to the material surface, of which X-ray absorption and X-ray magnetic circular dichroism (XMCD) are very powerful tools [9]. Diffraction-based stress analysis depends critically on the use of the correct diffraction elastic constants [10]. By comparing various calculation method to determine homogenized moduli of the polycrystalline material composed of a single crystal, for example, the certain stress of Reuss model [11], the certain strain of Voigt model [12] and Kröner-Voigt model [7] taking into account these two models and interaction between crystals, the Y parameter, in the theoretical calculation to forecast the elastic modulus of polycrystalline material, has a high consistency. Y parameter makes up the methodological disadvantage of X-ray diffraction method (only the stress determined by strain). It can be used to estimate the homogenized properties under various crystal planes, which is the function associated with the crystal plane index.

Y parameters of both cubic and hexagonal polycrystalline materials based on the Voigt model of the certain strain [3, 4] and Y parameters of only hexagonal polycrystalline materials based on the certain stress of the Reuss model [4] have been given. In this work, it is first proposed to complete this part of the theoretical analysis and to obtain the definition of the Y parameter based on the Reuss model [13]. Here, the expression form of the 6 × 6 matrix (36 parameters) [14] is used instead of the four-rank tensor (81 parameters) in order to simplify the actual computing work and the derivation [15]. Based on the content mentioned above, the concept of Y parameter, as new evaluation criteria, has been proposed on the basis of theoretical formulas’ derivation to explain both Reuss and Voigt models.

1.2 Macroscopic mechanical behavior of cubic and hexagonal polycrystals

Polycrystalline material composed of a single crystal, which generally exhibits anisotropic mechanics, with a probability distribution depicting upward almost the same as in three-dimensional space, exhibits an isotropic behavior at macroscopic scale. Using Voigt model, elastic stiffness tensor is obtained by Morris [16] and Sayers [17], using the volume averaging with setting orthogonal constitutive relation of cubic crystal. Based on Voigt and Reuss models, Li and Thompson [18] have deduced elastic stiffness with setting orthogonal constitutive relations of hexagonal crystal for determining the upper and lower bounds. Anderson et al. [19] have verified these bounds and obtained the orthogonal polycrystalline texture coefficient of hexagonal crystal using ultrasonic measurement. Then, the compressibility, shear modulus, bulk modulus, Young’s modulus, and Poisson’s ratio can be calculated based on the Voigt-Reuss-Hill approximation [20, 21].

Y parameter can be seen as a function of the crystal plane index. Its theoretical derivation is thus relative to the type of crystal. Only cubic and hexagonal structures have been concerned by this concept. Y parameter of cubic and hexagonal structures was first introduced and presented by Zheng et al. [3, 4] in 2009. Elasticity of single crystal and mechanical properties of polycrystalline material has been closely integrated. For that, a preliminary study by ab initio plane-wave pseudopotential density functional theory method is first performed for each of the considered crystalline structures mentioned in Chapter 3; then, this chapter focuses on some cubic, hexagonal, and monoclinic polycrystals to estimate their homogenized elastic moduli. In this part, the methodology of Y parameter is applied to typical cubic, hexagonal, and monoclinic structures, respectively, CaO and MgO, CH, and calcite CaCO₃. Elastic moduli of the monoclinic tobermorite and gypsum polycrystals structures are calculated by Reuss-Voigt-Hill estimation [13]. Based on the elastic constants of cubic and hexagonal crystals, the corresponding Y parameter, as the new evaluation criteria, has been proposed. Thus, the credibility of Y parameters for determining elastic moduli of both cubic and hexagonal structures is proved. Elastic moduli of the monoclinic tobermorite structures calculated by Reuss-Voigt-Hill estimation are found to be further investigated for an amorphous structure in Chapter 5.

2. Y parameter of cubic structure by Voigt and Reuss models

Voigt and Reuss models are separately based on the assumption of the uniform deformation for upper bound and the uniform stress for lower bound. After that, Hill proposed averaged Voigt and Reuss models with a certain value in engineering application. However, according to the elastic constant definition of Hill model, that is to say $S_{\mathit{\text{Hill}}}=\frac{1}{2}\left(S_{Re\mathit{uss}}+C_{\mathit{\text{Voigt}}}^{-1}\right)$ and $C_{\mathit{\text{Hill}}}=\frac{1}{2}\left(C_{\mathit{\text{Voigt}}}+S_{Re\mathit{uss}}^{-1}\right)$, it can be seen that the stiffness coefficient C and the compliance coefficient S do not satisfy the relation $C\cdot S=I$. In other words, the elastic constants of Hill model are not self-consistent. That is why the Y parameter has been introduced to redefine the bounds of cubic and hexagonal polycrystals.

For Y-parameter description, macroanalysis of a cubic and hexagonal polycrystalline system is shown in Figure 1. At L₃ direction, there exists an orientation relationship between the global coordinate system [S] and the experimental coordinate system [L]. If set the direction cosine of the direction L₃ within system [S] relative to l, m, n, one can get [3]:

$$\begin{array}{l}l=sin\psi cos\varphi \\ m=sin\psi sin\varphi \\ n=cos\psi \end{array}\hspace{1em}\mathrm{(1)}$$

Figure 1.

Macro analysis of a polycrystalline system. (a) Macro analysis of a polycrystalline system and (b) schematic shear crystal plane.

Therefore, whole strain ${\epsilon }_{\mathit{\varphi \psi }}$ at the direction L₃ is as [22, 23]:

$${\epsilon }_{\mathit{\varphi \psi }}=l^2{\epsilon }_{11}+m^2{\epsilon }_{22}+n^2{\epsilon }_{33}+\mathit{lm}{\gamma }_{12}+\mathit{mn}{\gamma }_{23}+\mathit{nl}{\gamma }_{31}\hspace{1em}\mathrm{(2)}$$

Therefore, whole stress ${\sigma }_{\mathit{\varphi \psi }}$ at the direction L₃ is as [22, 23]:

$${\sigma }_{\mathit{\varphi \psi }}=l^2{\sigma }_{11}+m^2{\sigma }_{22}+n^2{\sigma }_{33}+2\left(\mathit{lm}{\tau }_{12}+\mathit{mn}{\tau }_{23}+\mathit{nl}{\tau }_{31}\right)\hspace{1em}\mathrm{(3)}$$

For cubic and hexagonal crystals, when the polycrystalline material as a whole is isotropic, Lamé constant λ as a polycrystalline material parameter satisfies the relationship: $\lambda =\frac{2\mathit{\mu G}}{1-2\mu }$, $G=\frac{E}{2\left(1+\mu \right)}$; here, E, $\mu$, G are, respectively, Young’s modulus, Poisson’s ratio, and shear modulus of the polycrystalline materials. When the polycrystalline material as a whole is only under a normal stress, thus ${\tau }_{12}={\tau }_{23}={\tau }_{31}=0$; then, we can see that ${\gamma }_{12}={\gamma }_{23}={\gamma }_{31}=0$. If suppose $\varphi =0^{°}$, Eq. (3) can be simplified as [3]:

$${\sigma }_{\mathit{\varphi \psi }}=l^2{\sigma }_{11}+m^2{\sigma }_{22}+n^2{\sigma }_{33}=\lambda \left({\epsilon }_{11}+{\epsilon }_{22}+{\epsilon }_{33}\right)+2G\left[\left({\epsilon }_{11}-{\epsilon }_{33}\right){sin}^2\psi +{\epsilon }_{33}\right]\hspace{1em}\mathrm{(4)}$$

From Figure 1, for Reuss bounds of cubic crystal, when the isotropic material as a whole is only under a uniaxial stress (i.e., ${\sigma }_{11}\ne 0$, other stress components are zero) and under the conditions of $\varphi =0^{°}$, the ratio expression between the whole strain and the uniaxial stress from Eq. (4) will be eventually developed as [3]:

$$\frac{{\epsilon }_{\mathit{\varphi \psi }}}{{\sigma }_{11}}=\left(\frac{1+\mu }{E}\right){sin}^2\psi -\left(\frac{\mu }{E}\right)\hspace{1em}\mathrm{(5)}$$

2.1 Reuss bounds of cubic polycrystalline structure

For the Reuss bound of cubic polycrystalline structure, the content is detailed in our previous work [13]. The main crystallographic planes of a cubic crystal are shown in Figure 2.

Figure 2.

Some crystallographic planes of a cubic crystal. (a) The directions in a cubic lattice and (b) the {100}, {110}, and {111} planes.

As in Figure 2, if the uniaxial tensions are separately applied in the direction of [100], [110], and [111] for a cubic crystal, according to the Ref. [24], Young’s moduli by elastic constants C₁₁, C₁₂, and C₄₄ can be calculated as follows:

$${\left.E\right|}_{\left[100\right]}=\frac{c_{11}^2+c_{11}{c}_{12}-2c_{12}^2}{c_{11}+c_{12}}\hspace{1em}\mathrm{(6)}$$

$${\left.E\right|}_{\left[110\right]}=\frac{4\left(c_{11}^2+c_{11}{c}_{12}-2c_{12}^2\right)c_{44}}{c_{11}^2+c_{11}{c}_{12}+2c_{11}{c}_{44}-2c_{12}^2}\hspace{1em}\mathrm{(7)}$$

$${\left.E\right|}_{\left[111\right]}=\frac{3\left(c_{11}+2c_{12}\right)c_{44}}{c_{11}+2c_{12}+c_{44}}\hspace{1em}\mathrm{(8)}$$

For cubic crystal, Y elastic constants based on Reuss model which has contained crystal surface normal (n = (u, v, w)) are as:

$$G_Y\left(u,v,w\right)=\frac{1}{2\left(s_{11}-s_{12}\right)-3\left[2\left(s_{11}-s_{12}\right)-s_{44}\right]\mathrm{\Gamma }}\hspace{1em}\mathrm{(9)}$$

$$E_Y\left(u,v,w\right)=\frac{1}{s_{11}-\left[2\left(s_{11}-s_{12}\right)-s_{44}\right]\mathrm{\Gamma }}\hspace{1em}\mathrm{(10)}$$

$${\mu }_Y\left(u,v,w\right)=-\frac{2s_{12}+\left(2\left(s_{11}-s_{12}\right)-s_{44}\right)\mathrm{\Gamma }}{2\left\{s_{11}-\left[2\left(s_{11}-s_{12}\right)-s_{44}\right]\mathrm{\Gamma }\right\}}\hspace{1em}\mathrm{(11)}$$

where $\overline{\mathrm{\Gamma }}=u^2{v}^2+v^2{w}^2+w^2{u}^2$. From Eqs. (9)–(11), elastic moduli of cubic structures by Y parameters can be calculated, and the relationship of G-3Γ, E-3Γ, and μ-3Γ can be calculated and the limitation of elastic moduli can be determined by Reuss and Voigt models.

2.2 Voigt bounds of cubic polycrystalline structure

As is shown in Figure 2, supposing that the direction cosines are L₃ (u, v, w) relative to the cubic crystal coordinates [K], the relationship between the macroscopic stress σ^L3 in L₃ direction and the stress component (${\sigma }_x$, ${\sigma }_y$, ${\sigma }_z$, ${\tau }_{\mathit{yz}}$, ${\tau }_{\mathit{zx}}$, ${\tau }_{\mathit{xy}}$) of the coordinate system [K] is as [22, 23]:

$${\sigma }^{L_3}=u^2{\sigma }_x+v^2{\sigma }_y+w^2{\sigma }_z+2\left(\mathit{vw}{\tau }_{\mathit{yz}}+\mathit{wu}{\tau }_{\mathit{zx}}+\mathit{uv}{\tau }_{\mathit{xy}}\right)\hspace{1em}\mathrm{(12)}$$

With respect to the coordinate system [K], the relationship between stress and strain of single cubic crystal structure is as:

$$\left[\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\sigma }_z\\ {\tau }_{\mathit{yz}}\\ {\tau }_{\mathit{zx}}\\ {\tau }_{\mathit{xy}}\end{array}\right]=\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{12}& 0& 0& 0\\ c_{12}& c_{11}& c_{12}& 0& 0& 0\\ c_{12}& c_{12}& c_{11}& 0& 0& 0\\ 0& 0& 0& c_{44}& 0& 0\\ 0& 0& 0& 0& c_{44}& 0\\ 0& 0& 0& 0& 0& c_{44}\end{array}\right]\left[\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\epsilon }_z\\ {\gamma }_{\mathit{yz}}\\ {\gamma }_{\mathit{zx}}\\ {\gamma }_{\mathit{xy}}\end{array}\right]\hspace{1em}\mathrm{(13)}$$

Submitting Eq. (13) into Eq. (12), one can get:

$$\begin{array}{cc}{\sigma }^{L_3}& =c_{11}\left(u^2{\epsilon }_x+v^2{\epsilon }_y+w^2{\epsilon }_z\right)+c_{12}\left[\left(v^2+w^2\right){\epsilon }_x+\left(u^2+w^2\right){\epsilon }_y+\left(u^2+v^2\right){\epsilon }_z\right]\\ & +2c_{44}\left(\mathit{vw}{\gamma }_{\mathit{yz}}+\mathit{uw}{\gamma }_{\mathit{zx}}+\mathit{uv}{\gamma }_{\mathit{xy}}\right)\end{array}\hspace{1em}\mathrm{(14)}$$

With respect to the coordinate system [K], assuming the three direction cosines of coordinate axes [S] (respectively S₁, S₂, S₃) are S₁ (l₁, m₁, n₁), S₂ (l₂, m₂, n₂), S₃ (l₃, m₃, n₃), and considering relationship between strain components in coordinates system [S] and strain component of coordinates system [K], we can get relation ${\gamma }_{12}={\gamma }_{23}={\gamma }_{31}=0$ when the polycrystalline material is limited only by the principal stress (${\tau }_{12}={\tau }_{23}={\tau }_{31}=0$). So the following formula of strain holds [22]:

$$\begin{array}{l}{\epsilon }_x=l_1^2{\epsilon }_{11}+l_2^2{\epsilon }_{22}+l_3^2{\epsilon }_{33}\\ {\epsilon }_y=m_1^2{\epsilon }_{11}+m_2^2{\epsilon }_{22}+m_3^2{\epsilon }_{33}\\ {\epsilon }_z=n_1^2{\epsilon }_{11}+n_2^2{\epsilon }_{22}+n_3^2{\epsilon }_{33}\\ {\gamma }_{\mathit{yz}}=2\left(m_1{n}_1{\epsilon }_{11}+m_2{n}_2{\epsilon }_{22}+m_3{n}_3{\epsilon }_{33}\right)\\ {\gamma }_{\mathit{zx}}=2\left(n_1{l}_1{\epsilon }_{11}+n_2{l}_2{\epsilon }_{22}+n_3{l}_3{\epsilon }_{33}\right)\\ {\gamma }_{\mathit{xy}}=2\left(l_1{m}_1{\epsilon }_{11}+l_2{m}_2{\epsilon }_{22}+l_3{m}_3{\epsilon }_{33}\right)\end{array}\hspace{1em}\mathrm{(15)}$$

Substituting Eq. (15) into Eq. (14), one can get:

$$\begin{array}{c}\begin{array}{l}{\sigma }^{L_3}=c_{11}\left[u^2\left(l_1^2{\epsilon }_{11}+l_2^2{\epsilon }_{22}+l_3^2{\epsilon }_{33}\right)+v^2\left(m_1^2{\epsilon }_{11}+m_2^2{\epsilon }_{22}+m_3^2{\epsilon }_{33}\right)+w^2\left(n_1^2{\epsilon }_{11}+n_2^2{\epsilon }_{22}+n_3^2{\epsilon }_{33}\right)\right]\hfill \end{array}\\ \begin{array}{c}+c_{12}[\left(v^2+w^2\right)\left(l_1^2{\epsilon }_{11}+l_2^2{\epsilon }_{22}+l_3^2{\epsilon }_{33}\right)+\left(u^2+w^2\right)\left(m_1^2{\epsilon }_{11}+m_2^2{\epsilon }_{22}+m_3^2{\epsilon }_{33}\right)\end{array}\\ +\left(u^2+v^2\right)\left(n_1^2{\epsilon }_{11}+n_2^2{\epsilon }_{22}+n_3^2{\epsilon }_{33}\right)]\\ \begin{array}{l}+4c_{44}[vw\left(m_1{n}_1{\epsilon }_{11}+m_2{n}_2{\epsilon }_{22}+m_3{n}_3{\epsilon }_{33}\right)+uw\left(n_1{l}_1{\epsilon }_{11}+n_2{l}_2{\epsilon }_{22}+n_3{l}_3{\epsilon }_{33}\right)\hfill \end{array}\\ +uv\left(l_1{m}_1{\epsilon }_{11}+l_2{m}_2{\epsilon }_{22}+l_3{m}_3{\epsilon }_{33}\right)]\end{array}(16)$$

In practice, for cubic polycrystalline structure, the single crystal has a haphazard distribution, and the probability distribution of the crystals depicting up almost the same in the three-dimensional space. So, we can average all directions in three-dimensional space, and then obtain the mean stress $\overline{{\sigma }^{L_3}}$ so as to determine its macroscopic stress. Here, if $\varphi =0^{°}$, we can obtain the mean stress $\overline{{\sigma }_{Y^{L_3}}^{}}$ by averaging the orientation in the crystal plane:

$$\overline{{\sigma }_Y^{L_3}}=\raisebox{1ex}{${\int }_0^{2\pi }{\sigma }^{L_3}\mathit{d\theta }$}\!\left/ \!\raisebox{-1ex}{${\int }_0^{2\pi }\mathit{d\theta }$}\right.\hspace{1em}\mathrm{(17)}$$

Here, we use the certain strain of Voigt model, where the material is only under the uniaxial strain [12]. For polycrystalline, the stress can be obtained by averaging integral within a crystal surface according to the formula (17):

$$\begin{array}{l}\overline{{\sigma }_Y^{L_3}}=\left[c_{11}-c_{12}-3\left(c_{11}-c_{12}-2c_{44}\right)\mathrm{\Gamma }\right]\left({\epsilon }_{11}-{\epsilon }_{33}\right){sin}^2\psi +\left[c_{11}-c_{12}-3\left(c_{11}-c_{12}-2c_{44}\right)\mathrm{\Gamma }\right]{\epsilon }_{33}\\ +\left({\epsilon }_{11}+{\epsilon }_{22}+{\epsilon }_{33}\right)\cdot \left[c_{12}+\left(c_{11}-c_{12}-2c_{44}\right)\mathrm{\Gamma }\right]\end{array}\hspace{1em}\mathrm{(18)}$$

where $\mathrm{\Gamma }=u^2{v}^2+v^2{w}^2+u^2{w}^2,\left(u,v,w\right)$ is the cosine normals within planes of cubic crystal. With comparison of Eq. (18) and Eq. (16), we can obtain Y elastic constants containing crystal surface normal (n = (u, v, w)) as [3]:

$$G_Y\left(u,v,w\right)=\frac{c_{11}-c_{12}}{2}-\frac{3\left(c_{11}-c_{12}-2c_{44}\right)\mathrm{\Gamma }}{2}\hspace{1em}\mathrm{(19)}$$

$${\lambda }_Y\left(u,v,w\right)=c_{12}+\left(c_{11}-c_{12}-2c_{44}\right)\mathrm{\Gamma }\hspace{1em}\mathrm{(20)}$$

As $\lambda =\frac{2\mathit{\mu G}}{1-2\mu }$, $G=\frac{E}{2\left(1+\mu \right)}$, formulas can be developed as:

$${\left(\frac{1+\mu }{E}\right)}_Y=\frac{1}{c_{11}-c_{12}-\left(c_{11}-c_{12}-2c_{44}\right)\cdot 3\mathrm{\Gamma }}\hspace{1em}\mathrm{(21)}$$

$$-{\left(\frac{\mu }{E}\right)}_Y=-\frac{c_{12}+\left(c_{11}-c_{12}-2c_{44}\right)\mathrm{\Gamma }}{\left(c_{11}+2c_{12}\right)\left[c_{11}-c_{12}-3\left(c_{11}-c_{12}-2c_{44}\right)\mathrm{\Gamma }\right]}\hspace{1em}\mathrm{(22)}$$

$$E_Y\left(u,v,w\right)=\frac{\left(c_{11}+2c_{12}\right)\left[c_{11}-c_{12}-3\left(c_{11}-c_{12}-2c_{44}\right)\mathrm{\Gamma }\right]}{c_{11}+c_{12}-\left(c_{11}-c_{12}-2c_{44}\right)\mathrm{\Gamma }}\hspace{1em}\mathrm{(23)}$$

$${\mu }_Y\left(u,v,w\right)=\frac{c_{12}+\left(c_{11}-c_{12}-2c_{44}\right)\mathrm{\Gamma }}{c_{11}+c_{12}-\left(c_{11}-c_{12}-2c_{44}\right)\mathrm{\Gamma }}\hspace{1em}\mathrm{(24)}$$

Mechanical properties of polycrystalline can be obtained by Y elastic constant with averaging in the (u, v, w) entire 3-D space orientation, which are all elastic parameters having nothing to do with the crystal surface index. Assuming $\overline{\mathrm{\Gamma }}=1/5$ [25], the mechanical properties of polycrystalline, such as Young’s modulus E, Poisson’s ratio μ, and the shear modulus G, are as [3]:

$$G_M=\frac{c_{11}-c_{12}+3c_{44}}{5},{\lambda }_M=\frac{c_{11}+4c_{12}-2c_{44}}{5}\hspace{1em}\mathrm{(25)}$$

$${\left(\frac{1+\mu }{E}\right)}_M=\frac{5}{2\left(c_{11}-c_{12}+3c_{44}\right)}\hspace{1em}\mathrm{(26)}$$

$$-{\left(\frac{\mu }{E}\right)}_M=-\frac{c_{11}+4c_{12}-2c_{44}}{2\left(c_{11}+2c_{12}\right)\left(c_{11}-c_{12}+3c_{44}\right)}\hspace{1em}\mathrm{(27)}$$

$$E_M=\frac{\left(c_{11}+2c_{12}\right)\left(c_{11}-c_{12}+3c_{44}\right)}{2c_{11}+3c_{12}+c_{44}}\hspace{1em}\mathrm{(28)}$$

$${\mu }_M=\frac{c_{11}+4c_{12}-2c_{44}}{2\left(2c_{11}+3c_{12}+c_{44}\right)}\hspace{1em}\mathrm{(29)}$$

Eqs. (19)–(24) are the theoretical formula of mechanical elastic constants and Y parameters of cubic polycrystal structures based on a certain strain of Voigt model.

For cubic polycrystalline model, at first, we set a face normal as the basis axis on surface of a polycrystalline, and then average the normals in all 360° orientation within the polycrystalline surface. After obtaining Y elastic constants, we choose the average orientation of crystal surface normals on the entire three-dimensional space.

3. Y parameter of hexagonal structure by Voigt and Reuss models

Similarly for hexagonal crystal from Figure 1, for Voigt bounds, when the isotropic material as a whole is only under a uniaxial stress (i.e., ${\sigma }_{11}\ne 0$, other stress components are zero) and under the conditions of $\varphi =0^{°}$, the ratio expression between the whole strain and the uniaxial stress from Eq. (4) will be eventually developed as [4]:

$${\sigma }_{\mathit{\varphi \psi }}=\lambda \left({\epsilon }_{11}+{\epsilon }_{22}+{\epsilon }_{33}\right)+2G\left[\left({\epsilon }_{11}-{\epsilon }_{33}\right){sin}^2\psi +{\epsilon }_{33}\right]\hspace{1em}\mathrm{(30)}$$

3.1 Reuss bounds of hexagonal polycrystal structure

As in Figure 1, supposing that the direction cosines are L₃ (u, v, w) relative to the hexagonal crystal coordinates [K], the relationship between the macroscopic strain ${\epsilon }^{L_3}$ (equivalent to ${\epsilon }_{\mathit{\varphi \psi }}$) in L₃ direction and the strain component (${\epsilon }_x$, ${\epsilon }_y$, ${\epsilon }_z$, ${\gamma }_{\mathit{xy}}$, ${\gamma }_{\mathit{yz}}$, ${\gamma }_{\mathit{zx}}$) of the coordinate system [K] is as [22, 23]:

$${\epsilon }^{L_3}=u^2{\epsilon }_x+v^2{\epsilon }_y+w^2{\epsilon }_z+\mathit{uv}{\gamma }_{\mathit{xy}}+\mathit{vw}{\gamma }_{\mathit{yz}}+\mathit{wu}{\gamma }_{\mathit{zx}}\hspace{1em}\mathrm{(31)}$$

With respect to the coordinate system [K], the relationship between strain and stress of a single hexagonal crystal structure is as [23]:

$$\left[\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\epsilon }_z\\ {\gamma }_{\mathit{yz}}\\ {\gamma }_{\mathit{zx}}\\ {\gamma }_{\mathit{xy}}\end{array}\right]=\left[\begin{array}{cccccc}{s}_{11}& s_{12}& s_{13}& 0& 0& 0\\ s_{12}& s_{11}& s_{13}& 0& 0& 0\\ s_{13}& s_{13}& s_{33}& 0& 0& 0\\ 0& 0& 0& s_{44}& 0& 0\\ 0& 0& 0& 0& s_{44}& 0\\ 0& 0& 0& 0& 0& 2\left(s_{11}-s_{12}\right)\end{array}\right]\left[\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\sigma }_z\\ {\tau }_{\mathit{yz}}\\ {\tau }_{\mathit{zx}}\\ {\tau }_{\mathit{xy}}\end{array}\right]\hspace{1em}\mathrm{(32)}$$

Substituting Eq. (32) into Eq. (31), one can get:

$$\begin{array}{l}{\epsilon }^{L_3}=s_{11}\left(u^2{\sigma }_x+v^2{\sigma }_y+2\mathit{uv}{\tau }_{\mathit{xy}}\right)+s_{12}\left(u^2{\sigma }_y+v^2{\sigma }_x-2\mathit{uv}{\tau }_{\mathit{xy}}\right)\\ +s_{13}\left[\left(u^2+v^2\right){\sigma }_z+w^2\left(\sigma {}_x+\sigma {}_y\right)\right]+s_{33}{w}^2{\sigma }_z+s_{44}w\left(v{\tau }_{\mathit{yz}}+u{\tau }_{\mathit{zx}}\right)\end{array}\hspace{1em}\mathrm{(33)}$$

With respect to the coordinate system [K], assuming the three direction cosines of coordinate axes [S] (respectively S₁, S₂, S₃) are S₁ (l₁, m₁, n₁), S₂ (l₂, m₂, n₂), S₃ (l₃, m₃, n₃), the relationship between the stress components in the coordinate system [K] and the stress component of the coordinate system [S] is as follows [26]:

$${\sigma }_x=\sum l_i{l}_j{\sigma }_{\mathit{ij}}\hspace{1em}\mathrm{(34)}$$

$${\sigma }_y=\sum m_i{m}_j{\sigma }_{\mathit{ij}}\hspace{1em}\mathrm{(35)}$$

$${\sigma }_z=\sum n_i{n}_j{\sigma }_{\mathit{ij}}\hspace{1em}\mathrm{(36)}$$

$${\tau }_{\mathit{xy}}=\sum l_i{m}_j{\sigma }_{\mathit{ij}}\hspace{1em}\mathrm{(37)}$$

$${\tau }_{\mathit{yz}}=\sum m_i{n}_j{\sigma }_{\mathit{ij}}\hspace{1em}\mathrm{(38)}$$

$${\tau }_{\mathit{zx}}=\sum n_i{l}_j{\sigma }_{\mathit{ij}}\hspace{1em}\mathrm{(39)}$$

When the material is only under the uniaxial stress ${\sigma }_{11}\ne 0$ and other stress components are zero, under these conditions, substituting Eqs. (34)–(39) into Eq. (54), one can get:

$$\begin{array}{l}{\epsilon }^{L_3}=\left\{s_{11}\left(u^2{l}_1^2+v^2{m}_1^2+2{\mathit{uvl}}_1{m}_1\right)\right.+s_{12}\left(u^2{m}_1^2+v^2{l}_1^2-2{\mathit{uvl}}_1{m}_1\right)+\\ \left.s_{13}\left[\left(u^2+v^2\right)n_1^2+w^2\left(l_1^2+m_1^2\right)\right]+s_{33}{w}^2{n}_1^2+s_{44}w\left({\mathit{vm}}_1{n}_1+{\mathit{ul}}_1{n}_1\right)\right\}{\sigma }_{11}\end{array}\hspace{1em}\mathrm{(40)}$$

In practice, for hexagonal polycrystalline structure, the single crystal has a haphazard distribution, and the probability distribution of the crystals depicting up almost the same in the three-dimensional space. So, we can average all directions in three-dimensional space, and then obtain the mean strain $\overline{{\epsilon }^{L_3}}$ so as to determine its macroscopic strain. Here, if $\varphi =0^{°}$, we can obtain the mean strain $\overline{{\epsilon }_{Y^{L_3}}^{}}$ by averaging the orientation in the crystal plane:

$$\overline{{\epsilon }_Y^{L_3}}=\raisebox{1ex}{${\int }_0^{2\pi }{\epsilon }^{L_3}\mathit{d\theta }$}\!\left/ \!\raisebox{-1ex}{${\int }_0^{2\pi }\mathit{d\theta }$}\right.\hspace{1em}\mathrm{(41)}$$

Here, (u, v, w) are normal direction cosines of hexagonal crystal surface, so we take Reuss model using a certain stress. For polycrystalline, the stress can be obtained by averaging integral within a crystal surface according to the formula (41):

$$\begin{array}{l}\frac{\overline{{\epsilon }_Y^{L_3}}}{{\sigma }_{11}}=\left[\frac{1}{2}\left(2s_{11}-s_{12}-s_{13}\right)-\frac{w^2}{2}\left(5s_{11}-s_{12}-5s_{13}+s_{33}-3s_{44}\right)+\frac{3w^4}{2}\left(s_{11}-2s_{13}+s_{33}-s_{44}\right)\right]\\ \times {sin}^2\psi +\frac{1}{2}\left(s_{12}+s_{13}\right)+\frac{w^2}{2}\left(s_{11}-s_{12}-s_{13}+s_{33}-s_{44}\right)-\frac{w^4}{2}\left(s_{11}-2s_{13}+s_{33}-s_{44}\right)\end{array}\hspace{1em}\mathrm{(42)}$$

With comparison of Eq. (42) and Eq. (33), we can obtain Y elastic constants containing crystal surface normal (n = (u, v, w)) as:

$${\left(\frac{1+\mu }{E}\right)}_Y=\frac{1}{2}\left(2s_{11}-s_{12}-s_{13}\right)-\frac{w^2}{2}\left(5s_{11}-s_{12}-5s_{13}+s_{33}-3s_{44}\right)+\frac{3w^4}{2}\left(s_{11}-2s_{13}+s_{33}-s_{44}\right)\hspace{1em}\mathrm{(43)}$$

$$-{\left(\frac{\mu }{E}\right)}_Y=\frac{1}{2}\left(s_{12}+s_{13}\right)+\frac{w^2}{2}\left(s_{11}-s_{12}-s_{13}+s_{33}-s_{44}\right)-\frac{w^4}{2}\left(s_{11}-2s_{13}+s_{33}-s_{44}\right)\hspace{1em}\mathrm{(44)}$$

As, according to Eq. (43) and Eq. (44), we can get

$$E_Y\left(u,v,w\right)=\frac{1}{s_{11}-\left[2\left(s_{11}-s_{13}\right)-s_{44}\right]w^2+\left(s_{11}-2s_{13}+s_{33}-s_{44}\right)w^4}\hspace{1em}\mathrm{(45)}$$

$${\mu }_Y\left(u,v,w\right)=-\frac{s_{12}+s_{13}+\left(s_{11}-s_{12}-s_{13}+s_{33}-s_{44}\right)w^2-\left(s_{11}-2s_{13}+s_{33}-s_{44}\right)w^4}{2\left\{s_{11}-\left[2\left(s_{11}-s_{13}\right)-s_{44}\right]w^2+\left(s_{11}-2s_{13}+s_{33}-s_{44}\right)w^4\right\}}\hspace{1em}\mathrm{(46)}$$

$$G_Y\left(u,v,w\right)=\frac{1}{\left(2s_{11}-s_{12}-s_{13}\right)-\left(5s_{11}-s_{12}-5s_{13}+s_{33}-3s_{44}\right)w^2+3\left(s_{11}-2s_{13}+s_{33}-s_{44}\right)w^4}\hspace{1em}\mathrm{(47)}$$

Mechanical properties of polycrystalline can be obtained by Y elastic constant with averaging in the (u, v, w) entire 3-D space orientation, which are all elastic parameters having nothing to do with the crystal surface index. Assuming $\overline{w^2}=\frac{1}{3}$ and $\overline{w^4}=\frac{1}{5}$ [25], the mechanical properties of polycrystalline, such as Young’s modulus E, Poisson’s ratio μ, and the shear modulus G, are as [4]:

$${\left(\frac{1+\mu }{E}\right)}_M=\frac{7s_{11}-5s_{12}-4s_{13}+2s_{33}+3s_{44}}{15}\hspace{1em}\mathrm{(48)}$$

$$-{\left(\frac{\mu }{E}\right)}_M=\frac{s_{11}+5s_{12}+8s_{13}+s_{33}-s_{44}}{15}\hspace{1em}\mathrm{(49)}$$

$$E=\frac{15}{2\left(4s_{11}+2s_{13}+s_{44}\right)+3s_{33}}\hspace{1em}\mathrm{(50)}$$

$$\mu =\frac{s_{11}+5s_{12}+8s_{13}+s_{33}-s_{44}}{2\left(4s_{11}+2s_{13}+s_{44}\right)+3s_{33}}\hspace{1em}\mathrm{(51)}$$

$$G=\frac{15}{2\left(7s_{11}-5s_{12}-4s_{13}+2s_{33}+3s_{44}\right)}\hspace{1em}\mathrm{(52)}$$

Based on the Reuss bounds and elastic constants of a single crystal, the theoretical formulas (48)–(52) of mechanical elastic moduli of hexagonal polycrystalline are finally determined.

3.2 Voigt bounds of hexagonal polycrystal structure

With respect to the coordinate system [K], the relationship between strain and stress of a single hexagonal crystal is as in Eq. (54). By Eq. (54) and Ref. [3], one can get:

$$\begin{array}{l}{\sigma }^{L_3}=c_{11}\left(u^2{\epsilon }_x+v^2{\epsilon }_y+\mathit{uv}{\gamma }_{\mathit{xy}}\right)+c_{12}\left(u^2{\epsilon }_y+v^2{\epsilon }_x-\mathit{uv}{\gamma }_{\mathit{xy}}\right)\\ +c_{13}\left[\left(u^2+v^2\right){\epsilon }_z+w^2\left(\epsilon {}_x+\epsilon {}_y\right)\right]+c_{33}{w}^2{\epsilon }_z+2c_{44}w\left(v{\gamma }_{\mathit{yz}}+u{\gamma }_{\mathit{zx}}\right)\end{array}\hspace{1em}\mathrm{(53)}$$

Similar to Ref. [3], according to Voigt model [12], we can average oriented polycrystalline material within a crystal surface and obtain the mean stress:

$$\begin{array}{l}\overline{{\sigma }_Y^{L_3}}=\frac{1}{2}\left[2c_{11}-c_{12}-c_{13}-w^2\left(5c_{11}-c_{12}-5c_{13}+c_{33}-12c_{44}\right)+3w^4\left(c_{11}-2c_{13}+c_{33}-4c_{44}\right)\right]\\ \cdot \left({\epsilon }_{11}-{\epsilon }_{33}\right){sin}^2\psi \\ +\frac{1}{2}\left[2c_{11}-c_{12}-c_{13}-w^2\left(5c_{11}-c_{12}-5c_{13}+c_{33}-12c_{44}\right)+3w^4\left(c_{11}-2c_{13}+c_{33}-4c_{44}\right)\right]\\ \cdot {\epsilon }_{33}+\frac{1}{2}\left({\epsilon }_{11}+{\epsilon }_{22}+{\epsilon }_{33}\right)\\ \cdot \left[c_{12}+c_{13}+w^2\left(c_{11}-c_{12}-c_{13}+c_{33}-4c_{44}\right)-w^4\left(c_{11}-2c_{13}+c_{33}-4c_{44}\right)\right]\end{array}\hspace{1em}\mathrm{(54)}$$

According to Eq. (53) and Ref. [3], we can obtain Y elastic constants containing crystal surface normal (n = (u, v, w)) as:

$$G_Y\left(u,v,w\right)=\frac{1}{4}\left[2c_{11}-c_{12}-c_{13}-w^2\left(5c_{11}-c_{12}-5c_{13}+c_{33}-12c_{44}\right)+3w^4\left(c_{11}-2c_{13}+c_{33}-4c_{44}\right)\right]\hspace{1em}\mathrm{(55)}$$

$${\lambda }_Y\left(u,v,w\right)=\frac{1}{4}\left[c_{12}+c_{13}+w^2\left(c_{11}-c_{12}-c_{13}+c_{33}-4c_{44}\right)-w^4\left(c_{11}-2c_{13}+c_{33}-4c_{44}\right)\right]\hspace{1em}\mathrm{(56)}$$

As $G=\frac{E}{2\left(1+\mu \right)}$, $\lambda =\frac{2\mathit{\mu G}}{1-2\mu }$, formulas can be developed as follows [4]:

$${\left(\frac{1+\mu }{E}\right)}_Y=\frac{2}{2c_{11}-c_{12}-c_{13}-w^2\left(5c_{11}-c_{12}-5c_{13}+c_{33}-12c_{44}\right)+3w^4\left(c_{11}-2c_{13}+c_{33}-4c_{44}\right)}\hspace{1em}\mathrm{(57)}$$

$$-{\left(\frac{\mu }{E}\right)}_Y=-\frac{c_{12}+c_{13}+w^2\left(c_{11}-c_{12}-c_{13}+c_{33}-4c_{44}\right)-w^4\left(c_{11}-2c_{13}+c_{33}-4c_{44}\right)}{\left[c_{11}+c_{12}+c_{13}-w^2\left(c_{11}+c_{12}-c_{13}-c_{33}\right)\right]\cdot \left[2c_{11}-c_{12}-c_{13}-w^2\left(5c_{11}-c_{12}-5c_{13}+c_{33}-12c_{44}\right)+3w^4\left(c_{11}-2c_{13}+c_{33}-4c_{44}\right)\right]}\hspace{1em}\mathrm{(58)}$$

$$\begin{array}{l}{E}_Y\left(u,v,w\right)=\left[c_{11}+c_{12}+c_{13}-w^2\left(c_{11}+c_{12}-c_{13}-c_{33}\right)\right]\cdot \\ \frac{\left[2c_{11}-c_{12}-c_{13}-w^2\left(5c_{11}-c_{12}-5c_{13}+c_{33}-12c_{44}\right)+3w^4\left(c_{11}-2c_{13}+c_{33}-4c_{44}\right)\right]}{2c_{11}+c_{12}+c_{13}-w^2\left(3c_{11}+c_{12}-3c_{13}-c_{33}-4c_{44}\right)+w^4\left(c_{11}-2c_{13}+c_{33}-4c_{44}\right)}\end{array}\hspace{1em}\mathrm{(59)}$$

$${\mu }_Y\left(u,v,w\right)=\frac{c_{12}+c_{13}+w^2\left(c_{11}-c_{12}-c_{13}+c_{33}-4c_{44}\right)-w^4\left(c_{11}-2c_{13}+c_{33}-4c_{44}\right)}{2c_{11}+c_{12}+c_{13}-w^2\left(3c_{11}+c_{12}-3c_{13}-c_{33}-4c_{44}\right)+w^4\left(c_{11}-2c_{13}+c_{33}-4c_{44}\right)}\hspace{1em}\mathrm{(60)}$$

Mechanical properties of polycrystalline can be obtained by Y elastic constant with averaging in the (u, v, w) entire 3-D space orientation, which are all elastic parameters having nothing to do with the crystal surface index. Assuming $\overline{w^2}=\frac{1}{3}$ and $\overline{w^4}=\frac{1}{5}$ [25], the mechanical properties of polycrystalline, such as Young’s modulus E, Poisson’s ratio μ, and shear modulus G, are as [4]:

$$G_M=\frac{7c_{11}-5c_{12}-4c_{13}+2c_{33}+12c_{44}}{30}\hspace{1em}\mathrm{(61)}$$

$${\lambda }_M=\frac{c_{11}+5c_{12}+8c_{13}+c_{33}-4c_{44}}{15}\hspace{1em}\mathrm{(62)}$$

$${\left(\frac{1+\mu }{E}\right)}_M=\frac{15}{7c_{11}-5c_{12}-4c_{13}+2c_{33}+12c_{44}}\hspace{1em}\mathrm{(63)}$$

$$-{\left(\frac{\mu }{E}\right)}_M=\frac{3\left(c_{11}+5c_{12}+8c_{13}+c_{33}-4c_{44}\right)}{\left[2\left(c_{11}+c_{12}+2c_{13}\right)+c_{33}\right]\left(7c_{11}-5c_{12}-4c_{13}+2c_{33}+12c_{44}\right)}\hspace{1em}\mathrm{(64)}$$

$$E=\frac{\left[2\left(c_{11}+c_{12}+2c_{13}\right)+c_{33}\right]\left(7c_{11}-5c_{12}-4c_{13}+2c_{33}+12c_{44}\right)}{3\left[9c_{11}+5c_{12}+4\left(3c_{13}+c_{33}+c_{44}\right)\right]}\hspace{1em}\mathrm{(65)}$$

$$\mu =\frac{c_{11}+5c_{12}+8c_{13}+c_{33}-4c_{44}}{9c_{11}+5c_{12}+4\left(3c_{13}+c_{33}+c_{44}\right)}\hspace{1em}\mathrm{(66)}$$

Based on the Voigt bounds, elastic constants of a single crystal and the theoretical formulas (61)–(66), mechanical elastic moduli of polycrystalline hexagonal structure are finally determined.

For the relative formulas of cubic and hexagonal material, the elastic moduli (Young’s modulus E_φ and shear modulus G_φ) can be determined along any orientation, from the elastic constants or compliance constants [13].

For monoclinic structures, the relationship between elastic compliance and Young’s modulus E_Ф in either direction of crystal c-axis with the angle Ф can be described as [1]:

$$\begin{array}{l}\frac{1}{E_{\varphi }}=S_{11}{l}_1^4+2S_{12}{l}_1^2{l}_2^2+2S_{13}{l}_1^2{l}_3^2+2S_{15}{l}_1^3{l}_3+S_{22}{l}_2^4+2S_{23}{l}_2^2{l}_3^2+2S_{25}{l}_1{l}_2^2{l}_3\\ +S_{33}{l}_3^4+2S_{35}{l}_1{l}_3^3+S_{44}{l}_2^2{l}_3^2+2S_{46}{l}_1{l}_2^2{l}_3+S_{55}{l}_1^2{l}_3^2+S_{66}{l}_1^2{l}_2^2\end{array}\hspace{1em}\mathrm{(67)}$$

where $\left(\begin{array}{l}{l}_1\\ l_2\\ l_3\end{array}\right)=\left(\begin{array}{l}sin\theta cos\varphi \\ sin\theta sin\varphi \\ cos\theta \end{array}\right)$, S_ij is the single crystal compliance tensor in collapsed matrix notation and l₁, l₂, and l₃ are the direction cosines between the direction normal to the hkl plane and the lattice vectors (a, b, and c). So, for a given hkl plane, the elastic modulus E_hkl in a normal direction to hkl plane in polycrystalline can be determined from Eq. (67).

Meanwhile, for cubic and hexagonal material, the interchangeable relationship between elastic stiffness constants C_ij and the elastic compliance constants S_ij is detailed in Ref. [13].

4. Elastic moduli of cubic polycrystals structures at nanoscale

The homogenized moduli of typical cubic polycrystals can be calculated by Reuss-Voigt-Hill estimation [13].

4.1 Elastic moduli of typical cubic CaO structure

Elastic constants of CaO single crystal by GGA are as: c₁₁ = 200.879 GPa, c₁₂ = 56.902 GPa, and c₄₄ = 75.600 GPa. By LDA, the C_ij are as: c₁₁ = 289.497 GPa, c₁₂ = 71.432 GPa, and c₄₄ = 76.652 GPa. The elastic moduli of CaO polycrystalline based on Y parameters were calculated, shown in Figure 3.

Figure 3.

Poisson ratio, Young’s, and shear moduli of CaO by Y parameter.

Mechanical moduli of CaO polycrystalline by Reuss-Voigt-Hill method are shown in Table 1.

P	Bv	Br	Gv	Gr	B	G	E	μ
Ref. [27]	127.3333	127.3333	75.8000	66.8578	127.3333	71.3289	180.3170	0.2640
Ref. [28]	102	102	70.8000	70.3279	102	70.5639	172.0231	0.2189
Ref. [29]	114.0667	114.0667	83.9200	83.1886	114.0667	83.5543	201.4703	0.2056
Ref. [30]	109.6667	109.6667	84.4000	84.3972	109.6667	84.3986	201.5038	0.1938
0 (GGA)	104.894	104.894	74.155	74.113	104.894	74.134	179.9978	0.2140
0 (LDA)	144.120	144.120	89.604	86.985	144.120	88.295	219.9640	0.2456
0.1 (GGA)	105.011	105.011	74.334	74.300	105.011	74.317	176.617	0.2197
0.2 (GGA)	105.503	105.503	74.494	74.467	105.503	74.481	177.512	0.2196
0.3 (GGA)	105.616	105.616	74.629	74.608	105.616	74.619	178.159	0.2189
0.4 (GGA)	106.409	106.409	74.839	74.825	106.409	74.832	179.377	0.2190
0.5 (GGA)	107.046	107.046	75.016	75.006	107.046	75.011	180.357	0.2192
0.6 (GGA)	106.995	106.995	75.163	75.157	106.995	75.160	181.067	0.2180
0.7 (GGA)	107.443	107.443	75.319	75.315	107.443	75.317	181.873	0.2179
0.8 (GGA)	107.908	107.908	75.558	75.557	107.908	75.558	183.113	0.2172
0.9 (GGA)	108.600	108.600	75.730	75.730	108.600	75.730	184.108	0.2175
1.0 (GGA)	108.650	108.650	75.848	75.848	108.650	75.848	184.664	0.2167

Table 1.

Mechanical moduli of CaO polycrystalline by different methods.

From Figure 3(a)–(c), for CaO polycrystals, Y parameters of Reuss and Voigt models are very close to each other based on the elastic constants by GGA method. Besides, the Poisson ratio, Young’s, and shear moduli by GGA method have a very good agreement with other Refs. [27, 28, 29, 30]. By Y parameter, shear modulus is between 71.99 and 75.61 GPa, while Young’s modulus is between 175.78 and 182.86 GPa under overall crystal plane orientation.

From Table 1, the elastic moduli B and G of CaO can be obtained based on the dependences of the elastic constants C_ij. Mechanical moduli of CaO polycrystalline by different methods are shown in Table 1. Elastic moduli at 0 GPa are verified and averaged as: B = 104.894 GPa, G = 74.134 GPa, and μ = 0.2209 in Table 1. Young’s modulus is about 175.758 GPa by Reuss-Voigt-Hill estimation.

According to Eqs. (26)–(28), Young’s moduli of CaO at 0 GPa in directions of [100], [110], and [111] are determined as: E_[100] = 175.758 GPa, E_[110] = 181.037 GPa, and E_[111] = 182.868 GPa.

4.2 Elastic moduli of typical cubic MgO structure

Y parameters calculated by the polycrystalline bounds of MgO curves are shown in Figure 4.

Figure 4.

Poisson ratio, Young’s, and shear moduli of MgO by Y parameter.

As can be seen from Figure 4(a)–(c), MgO polycrystals, Y parameters of Reuss and Voigt models are very close to each other based on the elastic constants by GGA method. Besides, the Poisson ratio, Young’s, and shear moduli by GGA method have a very good agreement with other Refs. [28, 30, 31, 32]. Moreover, shear modulus of MgO by Reuss and Voigt bounds of Y parameter is between 118.78 and 140.96 GPa under overall crystal plane orientation. Moreover, Young’s modulus of MgO is between 293.62 and 337.21 GPa.

Mechanical moduli of MgO polycrystalline by Reuss-Voigt-Hill method are shown in Table 2.

P	Bv	Br	Gv	Gr	B	G	E	μ
Ref. [28]	138.6667	138.6667	115.6000	113.8262	138.6667	114.7131	269.7539	0.1758
Ref. [30]	164.0000	164.0000	132.6000	131.0638	164.000	131.8319	311.9172	0.1830
Ref. [31]	163.9333	163.9333	134.7860	128.7278	163.9333	131.7569	311.7504	0.1831
Ref. [32]	162.4000	162.4000	129.4800	125.6660	162.4000	127.5730	303.3001	0.1887
0 (GGA)	185.395	185.395	131.9160	130.9822	185.395	131.449	293.172	0.2364
0 (LDA)	160.988	160.988	125.974	123.037	160.988	124.506	254.969	0.2360
0.1 (GGA)	167.217	167.217	128.111	125.659	167.217	126.885	264.582	0.2363
0.2 (GGA)	165.796	165.796	128.265	125.871	165.796	127.068	264.972	0.2336
0.3 (GGA)	167.421	167.421	128.452	126.085	167.421	127.268	266.009	0.2352
0.4 (GGA)	168.718	168.718	128.673	126.362	168.718	127.517	267.279	0.2360
0.5 (GGA)	168.474	168.474	128.863	126.583	168.474	127.723	267.835	0.2350
0.6 (GGA)	168.583	168.583	129.006	126.766	168.583	127.886	268.484	0.2346
0.7 (GGA)	169.250	169.250	129.213	127.020	169.250	128.116	269.486	0.2346
0.8 (GGA)	168.603	168.603	129.409	127.258	168.603	128.334	270.054	0.2330
0.9 (GGA)	169.886	169.886	129.654	127.530	169.886	128.592	271.116	0.2340
1.0 (GGA)	170.181	170.181	129.830	127.738	170.181	128.784	271.836	0.2338

Table 2.

Mechanical moduli of MgO polycrystalline by different methods.

From Table 2, the MgO elastic moduli of B and G on the dependences of the elastic constants c_ij, the elastic anisotropic parameter A on pressure is discussed. Mechanical moduli of MgO polycrystalline by different methods are shown in Table 2. Elastic moduli of MgO at 0 GPa are verified and averaged as: B = 160.988 GPa, G = 124.506 GPa, and μ = 0.2360 in Table 2. Young’s modulus of MgO is about 254.969 GPa by Reuss-Voigt-Hill estimation.

Elastic constants of MgO single crystal by GGA are as: c₁₁ = 343.468 GPa, c₁₂ = 106.358 GPa, and c₄₄ = 140.823 GPa. By LDA, the C_ij are as: c₁₁ = 298.508 GPa, c₁₂ = 92.228 GPa, and c₄₄ = 141.197 GPa. By Eqs. (26)–(28), Young’s moduli of MgO at 0 GPa by GGA method in directions of [100], [110], and [111] are determined as: E_[100] = 293.173 GPa, E_[110] = 324.938 GPa, and E_[111] = 337.114 GPa.

5. Elastic moduli of typical hexagonal polycrystals at nanoscale

The homogenized moduli of typical cubic polycrystals can be calculated by Reuss-Voigt-Hill estimation [13].

5.1 Elastic moduli of hexagonal CH structures

Elastic constants of CH single crystal at 0 GPa have been calculated by LDA calculation method and are given in Chapter 3. Figure 5 also shows the anisotropy of the CH Young’s modulus.

Figure 5.

Poisson ratio, Young’s, and shear moduli of CH by Y parameter.

From Figure 5(a)–(c), Y parameters of Reuss and Voigt models are very close to each other based on the elastic constants by LDA method. Besides, the Poisson ratio, Young’s modulus, and shear modulus of CH by LDA method have a very good agreement with other Refs. [33, 34, 35, 36]. Moreover, shear modulus of CH by Reuss and Voigt bounds of Y parameter is between 6.35 and 40.53 GPa under overall crystal plane orientation. Moreover, it has been found that the Young’s modulus of CH is between 21.33 and 92.23 GPa in function of the crystal plane orientation. The results of Voigt and Reuss approximations show upper and lower limits in Figure 5; the difference arises due to different elastic constants of single CH crystals.

Mechanical moduli of CH polycrystalline by Reuss-Voigt-Hill method are shown in Table 3.

Pressure (GPa)	Bv (GPa)	Br (GPa)	Gv (GPa)	Gr (GPa)	B (GPa)	G (GPa)	E (GPa)	μ
Ref. [33]	36.2300	26.6316	22.6510	13.9196	31.4308	18.2853	45.9460	0.2564
Ref. [34]	37.3(±0.4)	26.0(±0.3)	24.4(±0.3)	17.5(±0.4)	31.6(±0.3)	20.9(±0.3)	51.4(±1.0)	0.23(±0.02)
Ref. [35]	45.9778	22.3513	26.4167	8.4413	34.1646	17.4290	44.6878	0.2820
Ref. [36]	40.2422	28.0648	21.2959	15.3032	34.1535	18.2996	46.5796	0.2727
0 (GGA)	36.6395	28.1827	27.0233	21.0793	32.4111	24.0513	58.0794	0.2074
0 (LDA)	30.8885	19.5826	22.9999	14.8900	25.2355	18.9450	45.45917	0.199767
0.5 (LDA)	35.8308	27.1002	24.0796	18.2410	31.4655	21.1603	51.85653	0.225326
1.0 (LDA)	40.5673	35.0865	24.6871	17.7577	37.8269	21.2224	53.63647	0.263676
1.5 (LDA)	42.7098	35.6819	26.0083	18.9866	39.1959	22.4975	56.65329	0.259102
2.0 (LDA)	44.8341	37.4192	28.4868	25.8209	41.1267	27.1539	66.76733	0.229424
2.5 (LDA)	47.3443	41.2657	29.6026	26.3588	44.3050	27.9807	69.34408	0.239141
3.0 (LDA)	51.2366	46.0601	27.9731	23.6034	48.6484	25.7882	65.74721	0.274754
3.5 (LDA)	54.4701	49.1269	27.9306	23.5463	51.7985	25.7384	66.24323	0.286856
4.0 (LDA)	56.4897	51.9437	28.8208	25.0783	54.2167	26.9495	69.35677	0.286791
4.5 (LDA)	55.9515	51.4301	31.3143	28.8379	53.6908	30.0761	76.03140	0.263984
5.0 (LDA)	65.1811	61.3509	30.6193	26.9974	63.2660	28.8084	75.03591	0.302327

Table 3.

Mechanical moduli of CH polycrystalline by different methods.

From Table 3, based on elastic constants by GGA method, elastic moduli of CH by Jia [13] are as: B = 32.41 GPa and G = 24.05 GPa. Thus, Young’s modulus is averaged to 58.08 GPa by RVH method. We can see that elastic moduli by LDA method are close to other Refs. [33, 34, 35, 36]. However, there is significant deviation with the results obtained by GGA method. Moreover, Y parameter is consistent with the results of Voigt and Reuss models. Based on elastic constants using LDA method, by Ref. [13], elastic moduli at 0 GPa are verified and averaged as: G_v = 30.889 GPa, B_v = 22.999 GPa, G_r = 19.583 GPa, B_r = 14.890 GPa, B = 25.236 GPa, G = 18.945 GPa, and μ = 0.200 in Table 3. Young’s modulus of CH is about 45.459 GPa by Reuss-Voigt-Hill estimation.

5.2 Elastic modulus of hexagonal calcite structure

For calcite single crystal, elastic constants have been calculated by GGA method as follows: c₁₁ = 143.938 GPa, c₁₂ = 55.579 GPa, c₁₃ = 50.970 Pa, c₃₃ = 81.175 GPa, and c₄₄ = 32.658 GPa. The results of Voigt and Reuss approximations showing upper and lower limits are given in Table 4; the difference arises due to the use of different elastic constants of single crystals.

Pressure (GPa)	Bv (GPa)	Br (GPa)	Gv (GPa)	Gr (GPa)	B (GPa)	G (GPa)	G/B	E (GPa)	μ
Ref. [37]	80.6433	74.6632	38.9690	36.5788	77.6533	37.7739	0.486443	97.5106	0.2907
Ref. [38]	77.8078	71.5529	36.7167	35.0768	74.6803	35.8967	0.480672	92.8185	0.2929
Ref. [39]	89.7778	83.8382	39.1667	37.4911	86.8080	38.3289	0.441536	100.2343	0.3076
Ref. [40]	79.3111	72.9228	37.3967	35.1543	76.1169	36.2755	0.476576	93.9083	0.2944
Ref. [41]	78.0556	72.5157	36.5933	34.4239	75.2856	35.5086	0.471652	92.0535	0.2962
0.0 (LDA)	78.4495	70.9425	38.9105	36.1199	74.6960	37.5152	0.5022	96.4061	0.285
0.0 (GGA)	76.0099	69.6255	36.0014	28.7137	72.8177	32.3576	0.4444	84.5492	0.306
0.5 (GGA)	78.8573	72.0996	36.9119	29.1686	75.4785	33.0402	0.4377	86.4992	0.309
1.0 (GGA)	80.8540	73.9455	37.0892	29.6910	77.3998	33.3901	0.4314	87.5769	0.311
1.5 (GGA)	82.7850	76.1680	37.4184	29.7426	79.4765	33.5805	0.4225	88.3047	0.315
2.0 (GGA)	86.2955	78.5694	38.3801	30.3575	82.4325	34.3688	0.4169	90.5253	0.317
2.5 (GGA)	89.5705	80.5595	36.2223	27.4455	85.065	31.8339	0.3740	84.9097	0.334
3.0 (GGA)	90.8761	82.5141	38.9433	28.9708	86.6951	33.9571	0.3917	90.1068	0.327
3.5 (GGA)	93.6468	84.6268	39.4405	29.4758	89.1368	34.4582	0.3866	91.5743	0.329
4.0 (GGA)	95.8936	86.6201	39.1752	29.1386	91.2568	34.1569	0.3743	91.1041	0.334
4.5 (GGA)	98.4435	88.9200	39.3356	28.2723	93.6818	33.8139	0.3609	90.5476	0.306
5.0 (GGA)	100.9450	90.7940	39.7741	27.9952	95.8695	33.8846	0.3534	90.9398	0.309

Table 4.

Mechanical moduli of calcite polycrystalline by different methods.

From Table 4, the results with Y parameters are consistent with the results of Voigt and Reuss models. It can also be found that E and G are satisfied with the relation G = E/[2(1 + μ)].

Poisson ratio, Young’s, and shear moduli of calcite structure by Y parameter are shown in Figure 6.

Figure 6.

Poisson ratio, Young’s, and shear moduli of calcite by Y parameter.

From Figure 6(a)–(c) of calcite polycrystals, Y parameters of Reuss and Voigt models are very close to each other based on the elastic constants by GGA method. Besides, the Poisson ratio, Young’s, and shear moduli by GGA method have a very good agreement with other Refs. [37, 38, 39, 40, 41]. Moreover, by Reuss and Voigt bounds of Y parameter, shear modulus is between 15.08 and 46.13 GPa and Young’s modulus is between 42.28 and 116.65 GPa.

Based on elastic constants by Ref. [13], elastic moduli of calcite at 0 GPa are averaged as: G_v = 36.001 GPa, B_v = 76.010 GPa, G_r = 28.714 GPa, B_r = 69.626 GPa, B = 72.818 GPa, G = 32.358 GPa, and μ = 0.306 in Figure 6. Young’s modulus is about 84.549 GPa by RVH estimation.

Moreover, for the homogenized moduli of monoclinic polycrystals (gypsum and 11 Å tobermorite) at larger scale, the content is detailed in our previous chapter work [42]. Based on the elastic constants of 11 Å tobermorite (Ca/Si = 0.67) calculated by DFT, Young’s modulus is homogenized about 79.512 GPa, which is close to the simulation result of Pellenq et al. [43] (89 GPa) and Shahsavari et al. [44] (78.939 GPa). However, Young’s moduli considering the ordered Si-chain at long-range are far away from the nanoindentation test by Ulm [45, 46]. That confirms another time that the absence of order at long range in this phase and that the upscaling to polycrystals cannot be done with the tobermorite model. This finding is in agreement with the results of Manzano et al. [47].

6. Conclusions

The elastic properties of typical crystals are investigated and C_ij determination is given by DFT method. Y parameters have then been determined for various crystal structures and can be seen as an intermediate step in the homogenization. By means of the Y parameter, we can obtain the Young’s modulus in function of the orientation of the crystal plane. Contrary to Hill approach to obtain the isotropic elastic properties of polycrystals, the Y-parameter method enables to study the anisotropic behavior of a single crystal. By means of the Y parameter, which is a function of the stiffness coefficient C_ij (or the compliance coefficient S_ij) and the crystal plane orientation, the Young?s modulus in function of the orientation of the crystal plane may be obtained. Contrary to Hill approach which is used to obtain the isotropic elastic properties of a polycrystals, the Y parameter method enables to study the anisotropic behavior of a single unit cell. When $\overline{\mathrm{\Gamma }}=1/5$ or w² = 1/3, the result is equal to that of Hill model. Results are as follows:

1. For cubic CaO structure, Young’s modulus at 0 GPa is about 175.758 GPa by elastic constants and RVH estimation. Elastic moduli are as: B = 104.894 GPa, G = 74.134 GPa, and μ = 0.2209. Young’s moduli in directions of [100], [110], and [111] are separately determined by elastic constants as: E_[100] = 175.758 GPa, E_[110] = 181.037 GPa, and E_[111] = 182.868 GPa.

2. For cubic MgO structure, Young’s modulus of at 0 GPa is about 293.172 GPa by elastic constants and RVH estimation. Elastic moduli are homogenized as: B = 185.395 GPa, G = 131.449 GPa, and μ = 0.2364. Young’s moduli at 0 GPa in directions of [100], [110], and [111] are separately determined as: E_[100] = 293.173 GPa, E_[110] = 324.938 GPa, and E_[111] = 337.114 GPa.

3. For hexagonal CH structure, Young’s modulus of CH at 0 GPa is 45.459 GPa by LDA and RVH methods. Other elastic moduli are as: G_v = 30.889 GPa, B_v = 23.000 GPa, G_r = 19.583 GPa, B_r = 14.890 GPa, B = 25.236 GPa, G = 18.945 GPa, and μ = 0.200. On the other hand, by GGA and RVH methods, Young’s modulus is 58.07 GPa. Other elastic moduli are as: G_v = 27.023 GPa, B_v = 36.640 GPa, G_r = 21.079 GPa, B_r = 28.183 GPa, B = 32.41 GPa, G = 24.05 GPa, and μ = 0.207.

4. For hexagonal calcite structure, Young’s modulus at 0 GPa is about 84.542 GPa. Elastic moduli are homogenized as: G_v = 36.001 GPa, B_v = 76.010 GPa, G_r = 28.714 GPa, B_r = 69.626 GPa, B = 72.818 GPa, G = 32.358 GPa, and μ = 0.306. Young’s modulus is 84.549 GPa by RVH method.

It is worth to mention that the crystalline tobermorite structure cannot accurately represent the layered C▬S▬H structure at upper scale. So, modeling of a represented amorphous C▬S▬H structure-considered disordered Si chain, instead of the crystalline one, will be realized and further discussed based on MD method in the following Chapter 5.

Acknowledgements

The author acknowledges the financial support provided by the China Scholarship Council and the startup foundation of Xi’an Shiyou University. Thanks Qiufeng Wang for her proofreading, and very sincere thankness to prof. Fabrice Bernard and prof. Siham Kamali-Bernard for the their enthusiastic and responsible guideness.

Classical Molecular Dynamics (MD): Atomistic Simulation of Typical C▬S▬H Structures

Abstract

C▬S▬H has a great influence on the mechanical properties of cement and it is the most important binding phase of cement paste. In the light of recent computational material technology, a promising way to establish a larger microscale structure to obtain the desired C▬S▬H phase structure with the consideration of porosity is proposed, using a mechanical property to verify the rationality of the designed structure. Here, we discuss the basic atomic unit at nanoscale and the inverse approach, which propose strategies for the design of a possible C▬S▬H structure over nano- to macroscales of an achievement of nearly experimental phase. At the nanoscale, a 11 Å tobermorite monoclinic crystal with the ordered Si-chain is used to enlarge the scale to the monolithic “globule” C▬S▬H structure about 5.5 nm³ using molecular dynamics simulation. However, the inverse approach with full structural C▬S▬H information (atoms and their positions) from the macroscale by the neutral scattering tests to nanoscale by the reorganized amorphous C▬S▬H cell is unlinked nowadays. Our contribution is to find the interlinkages between the “globule” at nanoscale with LD/HD C▬S▬H at microscale, thus to seek experimental validation of the C▬S▬H hydrated phases.

Keywords: nanoscale, tobermorite, molecular dynamics, DFT, anisotropic elasticity, Young’s moduli

1. Introduction

Concrete is the most used building material in the world, because of its appreciable advantages such as load-bearing capacity, durability, maneuverability, flowability. Recently, a platform called multiscale modeling of computational concrete (MuMoCC) using various modeling methods has been developed. In these previous works, cementitious materials have been studied at four different scales:

1. The nanoscale that considers the element cells of the various components of the hardened cement paste [1, 2, 3]. The anisotropic elastic properties of different hydrated phases are determined thanks to various modeling methods [4]: density functional theory, molecular dynamics, atomic finite element method.

2. The microscale that considers the hardened cement paste as a heterogeneous and complex porous material in which the main solid phases are calcium silicate hydrates (C▬S▬H), portlandite (calcium hydroxides or CH), and hydrated aluminates or sulfoaluminates (ettringite or Aft, monosulfoaluminate or AFm) phases. The solid phase is itself in equilibrium with the internal solution, which fills the porosity [5].

3. The mesoscale that is in fact divided into two scales: the submesoscale, firstly, is the scale of mortar, that is, the smallest fraction of the particle size distribution bonded by homogeneous cement paste with porosity and air voids; then the mesoscale itself is the scale of concrete: the largest aggregates embedded in a matrix of mortar. These two scales introduce also the presence of the interfacial transition zone between particles and matrix, which modifies crack patterns and diffusion paths but not necessarily homogeneous effective physical and mechanical properties at the various scales [6, 7, 8].

In such a hierarchical approach, the adequate outputs at a given scale are used as input data at the next higher level.

C▬S▬H as the main and typical constituents of hydrated cementitious systems separately represents approximately 70% of the volume fraction of portland cement paste [9], which has a great influence on the mechanical properties of the cement paste, and it is the most important binding phase of cement paste. Its Young’s moduli are needed in the modeling of cement at the macro- and microscales [10]. The traditional continuum models [11] and nonlocal continuum theory [12] are not adequate in modeling [13] of these materials [14]. As C▬S▬H theoretical mechanism is very complex [15] and size effect [16, 17] are not commonly considered, the simulation is used for comparison of mechanical properties of C▬S▬H according to loading/unloading curves between load and displacement in nanoindentation technique. Although there are some relative reports on experimental analysis and the numerical calculation of C▬S▬H structure, its nanostructure has not yet been revealed.

The crystal structure of C▬S▬H structure is basically known and commonly modeled as tobermorite-like (i.e., tobermorite 9, 11, 14 Å) and jennite-like systems [18] and/or with distorted semicrystalline variations of them [19]. Tobermorite as one of the earliest models proposed by Taylor and Howison [20] is thought as its hydration degree, which may describe the relative C▬S▬H nanostructure. From the view of Ca/Si ratio, there are two types of C▬S▬H structures: the C▬S▬H(I) (with Ca/Si of 0.6–1.5) and C▬S▬H(II) structures [21] (with Ca/Si of 1.5–2.0), of which the later is close to the experimentally confirmed structure C_1.7▬S▬H_1.8 with Ca/Si of 1.7 in the C▬S▬H gel [22]. From the view of several nanometers length scale, C▬S▬H gel has a nanogranular aspect composed of monolithic C▬S▬H (full-dense C▬S▬H) and porosity. At this scale, C▬S▬H gel “globule” model exists in two forms: low-density C▬S▬H (LD C▬S▬H) and high-density C▬S▬H (HD C▬S▬H). Recently, Fabrice et al. have investigated the specific heat of cementitious materials using multiscale modeling approach [23]. Hou et al. [24] have investigated an amorphous C▬S▬H structure on the elastic properties of the layered C▬S▬H based on the CSH-FF field.

In recent years, some scholars have successfully simulated C▬S▬H structure by MD simulations. Pellenq et al. [25] have simulated several tobermorite structures by using the General Utility Lattice Program (GULP), where GULP uses the core-shell potential. Kalinichev et al. [26] have simulated the tobermorite structure using ClayFF force field [27], where water molecules are modeled using simple point charge (SPC) model, and it is indicated that the water shows very strong binding force in the solid surface. Pellenq et al. [28] have developed a new modeling method to simulate C▬S▬H gel structure, of which it has a shorter [SiO₄]⁴⁻ chain under the system equilibrium structure C_1.65▬S▬H_1.75, comparing with experimental values of C_1.7▬S▬H_1.8 [22, 29]. Shahsavari et al. [30] proposed a new CSH-FF field with a higher accuracy and more efficient than the traditional core-shell potential. Ji et al. [31] have simulated the influencing role of a variety of water molecules (usually optional SPC model) on C▬S▬H by using GROMACS [32] program and mainly used a simple harmonic force potential description of C▬S▬H between atoms and ions. Zaoui [33] has simulated tobermorite structure of different Ca/Si ratios under various pressures, of which Buckingham potential and core-shell potential are used. This author found that the elastic modulus of different Ca/Si ratio tobermorite structures tends to be similar with increasing pressure. Some scholars [34] develop the DL_POLY program to simulate the tobermorite by MD, finding that amorphous form of C▬S▬H(I) is a short-range structure and long-range disordered structure, which is more consistent with the experimental parameters. Dai et al. [35] have compared two C▬S▬H structures of different Ca/Si ratios and discussed MD modeling parameters, of which the COMPASS force field is used. According to the NMR measurement of C▬S▬H [19], mechanical properties of this amorphous structure can be simulated by MD method to study the deformation mechanism at atomic scale.

Moreover, microporomechanics technique has been used to calculate elastic properties of LD C▬S▬H and HD C▬S▬H gels [36]. For the elastic properties of LD and HD C▬S▬H, Constantinides and Ulm [37] have obtained elastic moduli by nanoindentation experiment. The porosity is also a factor of influencing elastic modulus of cement paste, and it has been investigated based on the backscattered electron image analysis and the HYMOSTRUC model [38]. The impact of the porosity on the matrix Young’s modulus is that the Young’s modulus decreases highly with the increasing of the porosity [39]. According to LD and HD C▬S▬H models described by Jennings [22], the gel porosity of LD C▬S▬H solid phase is 35–37%, while the gel porosity of HD C▬S▬H solid phase is 24%. However, the relationship between the C▬S▬H structure with the size of about 5 nm and C▬S▬H phases (LD and HD C▬S▬H) has not been revealed yet. In this chapter, the porosity and homogenization are used to explain the difference, which is meaningful for the verification of both nanoindentation simulation [40] and nanoindentation experiment [41, 42] at nanoscale in Chapter 6.

2. Molecular dynamics modeling of typical C▬S▬H structures by LAMMPS

Molecular dynamics is a deterministic method that offers the possibility of a microscopic description of a physical system in consideration of all the interactions involved. The main advantage of this method is that it gives the information on the evolution of the system over time and this by numerically solving Hamilton equations of motion, Lagrange or Newton.

2.1 Principle of MD method and ClayFF field used in LAMMPS

2.1.1 Introduction and principle of molecular dynamics method

The applied numerical integration for Newton’s equation of motion by MD [43] is described as:

$$H=\frac{1}{2}\sum _{i=1}^N\frac{p_i^2}{m_i}+\sum _{i=1}^{N-1}\sum _{j=i+1}^{N}U\left(r_{\mathit{ij}}\right)\hspace{1em}\mathrm{(1)}$$

$${\mathbf{p}}_i=m_i\frac{d{\mathbf{r}}_i}{\mathit{dt}}=m_i{\mathbf{v}}_i\hspace{1em}\mathrm{(2)}$$

$$\frac{d{\mathbf{p}}_i}{\mathit{dt}}=\sum _{i=1}^{N-1}\sum _{j=i+1}^N\mathbf{F}\left(r_{\mathit{ij}}\right)=-\sum _{i=1}^{N-1}\sum _{j=i+1}^N\frac{\partial U\left(r_{\mathit{ij}}\right)}{\partial {\mathbf{r}}_{\mathbf{ij}}}\hspace{1em}\mathrm{(3)}$$

where p is momentum that depends on the mass m and the acceleration a, F is the net force vector on the ith molecule, r is the position vector of the ith particle. The force on each molecule is obtained from the interaction potential, that is, the Lennard-Jones pair potential.

Positions and velocities of the particles are computed by integrating the equations of motion using finite difference algorithm at equal time intervals. One of the most used algorithms are the ones developed by Verlet that is the velocity Verlet leapfrog algorithm, which makes use of the half-time step velocities to calculate positions and velocities of particles [44]:

$${\mathbf{r}}_{\mathbf{i}}\left(\mathrm{t}+\mathrm{\Delta t}\right)={\mathbf{r}}_{\mathbf{i}}\left(\mathrm{t}\right)+v_{\mathbf{i}}\left(\mathrm{t}+\mathrm{\Delta t}/2\right)\mathrm{\Delta t}\hspace{1em}\mathrm{(4)}$$

$${\mathbf{v}}_{\mathbf{i}}\left(\mathrm{t}+\mathrm{\Delta t}/2\right)={\mathbf{v}}_{\mathbf{i}}\left(\mathrm{t}‐\mathrm{\Delta t}/2\right)+\frac{d^2{\mathbf{r}}_i\left(t\right)}{{}^{}}$$