Open access peer-reviewed compact

Atomistic Simulation of Anistropic Crystal Structures at Nanoscale

By Jia Fu

Reviewed: January 21st 2019Published: May 10th 2019

DOI: 10.5772/intechopen.84597

Downloaded: 318

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–10−7 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−4 m), the mesolevel (about 10−6–10−4 m), the microlevel (about 10−7–10−6 m), and atomic/nanolevel (about 10−10–10−7 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 methodsLength scales (m)Time scales (s)
Quantum mechanics method<10−810−15–10−12
Molecular dynamics method10−8–10−610−12–10−9
Micromechanics methods10−6–10−410−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 E8 model [18]. Garrett Lisi believes that all known physics particles and elementary various forces can be accommodated at the E8 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=Ur+Uθ+Uϕ+Uω+Uvdw+Ucoul(1)

where Ur is bond stretching, is bond angle bending, is dihedral angle torsion, is improper torsion, Uvdw is nonbonded van der Waals interaction, and Ucoul 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−1 m) down to the C▬S▬H solid phase (10−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 (C3S), dicalcium silicate (C2S), tricalcium aluminate (C3A), tetracalcium aluminoferrite (C4AF), 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 C4(A,F)H13. 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 C3S are noted. Aluminates are as a “stick” structure between the crystals C3S calcium silicates and/or C2S. Alkali sulfates are condensed on the surfaces of the silicate crystals during solidification process. The hydration of C2S and C3S 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), SiO2 (S), Al2O3 (A), Fe2O3 (F), H2O (H), SO3 (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.

σ1σ2σ3σ4σ5σ6=c11c12c13c14c15c16c21c22c23c24c25c26c31c32c33c34c35c36c41c42c43c44c45c46c51c52c53c54c55c56c61c62c63c64c65c66ε1ε2ε3ε4ε5ε6(2)

Elastic constant matrix Cijkl 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 typeNumber of CijIndependent elastic coefficient
Isotropic2c11, c12
a = b = cα = β = γ = 90°Cubic3c11, c12, c44
a = b ≠ cα = β = 90°, γ = 120°Hexagonal5c11, c12, c13, c33, c44
a = b ≠ cα = β = γ = 90°Tetragonal6c11, c12, c13, c33, c44, c66
a ≠ b ≠ cα = β = 90°, γ = 120Trigonal6c11, c12, c13, c33, c14, c44
a ≠ b ≠ cα = β = γ = 90°Orthorhombic9c11, c22, c33, c12, c13, c23, c44, c55, c66
a ≠ b ≠ cα = γ = 90°, β ≠ 90°Monoclinic13c11, c22, c33, cl2, c13, c23, c44, c55, c66, cl6, c26, c36, c45
a ≠ b ≠ cα ≠ γ ≠ β ≠ 90°Triclinic21c11, 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 typeLattice relationsCij number and elastic coefficientMethodObject
Cubica = b = c; α = β = γ = 90°3c11, c12, c44DFTCaO, MgO
Hexagonala = b ≠ c; α = β = 90°, γ = 120°5c11, c12, c13, c33, c44DFTCalcite, portlandite
AFEMPortlandite, graphene
Monoclinica ≠ b ≠ c; α = γ = 90°, β ≠ 90°13c11, c22, c33, c12, c13, c23, c44, c55, c66, cl5, c25, c35, c46DFTTobermorite, gypsum
MDTobermorite
OthersAFEMCNT
MDC▬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 (ε1, ε2, ε3, ε4, ε5, ε6) nonzero each time. Using the average stress tensor,<σij>, one can then define the effective elastic tensor by <σij>=<Cijkl>ε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.

MethodThe method of measurement or calculationMeasured or calculated parameters
ExperimentStatic, quasistatic methodStrain, stress
Indentation methodStrain, stress
Pulse-echo methodVelocity
Resonance ultrasonic spectrum (RUS) methodFrequency
X-ray or neutron methodDiffraction or scattering spectra
Brillouin scattering methodDiffraction or scattering spectra
CalculationPending elasticity parameter methodParameters from experiment
Classical molecular dynamicsEnergy, stress, strain
Monte Carlo methodEnergy, stress, strain
First-principles molecular dynamicsEnergy, stress, strain
DFT-QHA methodPhonon density of states, Helmholtz free energy
(Present)Atomic finite element methodStrain, 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.

MethodMeasurement or calculation methodsMeasured or calculated parameters
ExperimentIndentation methodStrain, stress
Classical numerical calculationClassical molecular dynamicsEnergy, stress, strain
First-principles molecular dynamicsEnergy, stress, strain
Method developedAFEMStress, 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 Ur, , , and 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 Uvdw 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 Cijkl in general. One can analytically average this tensor, using the full Cijkl 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 Cij 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.

YearAuthorMethodt (nm)E (TPa)E.t
1996Yakobson et al. [5]MD0.0665.50.3630
1999Goze et al. [9]Tight binding0.341.260.4284
2000Popov et al. [10]FCM0.331.020.3368
2000Zhou et.al [7]Tight-binding MD0.0745.10.3774
2001Kudin et al. [11]Ab initio0.0893.8590.3450
2002Tu and Ou-Yang [12]LDA0.0754.70.3525
2003Li and Chou [8]Structure mechanics0.341.020.3463
2003Chang and Gao [13]Molecular mechanics0.360.9970.3592
2004Sears and Batra [14]Molecular mechanics0.1342.520.3377
2004Pantano et al. [15]Continuum shell0.0754.840.3630
2005Zhang et al. [16]C-B rule0.3351.080.3618
2006To[17]FEM0.341.0240.3482
2007Chandraseker and Mukherjee [18]Ab initio0.3350.990.3317
2008Wang and Zhang [19]Continuum shell0.13.50.350
2008Jinan Lu and Chen [20]FEM0.341.2740.4332
2010Mahmood et al. [4]Continuum mechanics0.341.033–1.0420.3512–0.3543
2012Xiaoxing and Zhong Hu [21]Equivalent-continuum0.331.0580.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=Ur+Uθ+Uϕ+Uω+Uvdw(1)

where Ur is for a bond stretch interaction, for a bond angle bending, for a dihedral angle torsion, for an improper (out-of-plane) torsion, Uvdw 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:

Ur=12krrr02=12krΔr2(2)
Uθ=12kθθθ02=12kθΔθ2(3)
Uτ=Uϕ+Uω=12kτΔϕ2(4)

where kr, , and 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 Ix=Iy=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:

EAL=kr,EIL=kθ,GJL=kτ(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, kr is the stretching constant, is the bending constant, and 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 kr, , and 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: Fxi, Fyi, Mi for node i and Fxj, Fyj, Mj for node j. There are three nodal displacements for each node: ui, vi, θi for node i and uj, vj, θj for node j.

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

ΔFxi=krΔuiΔuj(6)
ΔFxj=krΔui+Δuj(7)
ΔMi=2kθiΔθLΔθi(8)
ΔMj=2kθjΔθjΔθL(9)

where ΔθL=ΔvjΔviL, kr is the stretching stiffness, kyi and kyj 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: kr=2.78aJ/Å2, kθ=0.498aJ/rad, ε=0.585×103aJ, r0=1.47×1010m, θ0=1.88rad, σ=3.53×1010m, kr=938kcal×mol1×Å2, kθ=126kcalmol1rad2, kτ=40kcalmol1rad2, aCC=L=0.142nm. For the round section of C▬C bonds, when t = 0.34 nm, thus A = 0.090792 (nm2), E = 1019.259 (GPa), and I = 1.21959 × 10−4 (nm4).

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 Ch=na1+ma2, a1=a2=3acc, so Ch=3accn2+nm+m2; chiral angle cosθ=2n+m2n2+nm+m2or θ=tg13mm+2n; diameter d=Chπ=3πaccn2+nm+m2.

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:

Etotalx=Utotalxi=1NF¯ixi(10)

where xiis the coordinate of atom i vector and F¯iis the external force exerted on the atom i vector. The minimum energy state is:

Etotalxx=0(11)

For Ftotalx, the initial state x0can be expanded by the developed Taylor, we obtain:

EtotalxEtotalx0+Etotalxxx=x0xx0+12xx0T2Etotalxxxx=x0xx0(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:

ψxn=dEtotaldxn+F¯(13)
ψxn+1=ψxn+dxnΔxn=0(14)

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

Δxn=dx1ψxn(15)
xn+1=xn+Δxn(16)

where F¯is the external force vector, ψxnis imbalance force vector after n iterations, xnis 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=dx=d2Etotaldx2(17)
F=dEtotaldx+F¯.(18)

The above equation can be written as:

Ku=P(19)

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

K=2Etotalxxx=x0=2Utotalxxx=x0(20)
P=Etotalxx=x0=F¯Utotalxx=x0(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:

KA=2EtotalaAaA2×2122EtotalaAai2×18122EtotalaiaA18×2018×18(22)

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

FA=F¯AEtotalaA2×1018×1(23)

where aArepresents the coordinates of atomic vector, and in the two-dimensional case, aA=xAyAT; the range of i is B-J, and the sum of the order number of KAand FAis related to the interaction between atoms A and B.

For FA, the component EtotalaA=EtotalxAEtotalyAT; For KA, the component 2EtotalaAaA.

For component FA:

UrALxA=UrALrALrALxA(24)
UrALyA=UrALrALrALyA(25)

For component KA:

U2rALxAxA=2UrALrAL2rALxA2+UrALrAL2rALxA2(26)
U2rALyAyA=2UrALrAL2rALyA2+UrALrAL2rALyA2(27)
U2rALxAyA=U2rALyAxA=2UrALrAL2rALxArALyA+UrALrAL2rALxAyA(28)

As is known to us, UVdWr=4εσrij12σrij6, so:

UrALrAL=4ε12×σ12rAL13+6×σ6rAL7(29)
2UrALrAL2=4ε156×σ12rAL1442×σ6rAL8(30)
rALxA=xAxLrAL(31)
2rALxA2=1rALxAxLrAL2rALxA(32)
rALyA=yAyLrAL(33)
2rALyA2=1rALyAyLrAL2rALyA(34)
2rALxAyA=2rALyAxA=xAxLyAyLrAL2(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=Ur+Uθ+Uϕ+Uω+Uvdw+Ucoul(36)

where Ur is for a bond stretch interaction, for a bond angle bending, for dihedral angle torsion, for an improper (out-of-plane) torsion, Uvdw for a nonbonded van der Waals interaction, and Ucoul for Coulombic (electrostatic) interactions. Equivalent beam elements [8] can be used to represent covalent bond to calculate Ur, , , and 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 Uvdw and Ucoul of ionic bond. Around the equilibrium position, van der Waals force is developed by the first-order Taylor:

FVdWri=FVdWr0+dFr0dririr0=25εσ2σr013σr0724εσ26σr0147σr08rir0(37)

where r0 (unit: Å) is the distance between two atoms around the equilibrium position without external force, ri 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:

FVdWatomAatomBr=24εABσAB2σABrABr0AB13σABrABr0AB7FVdWatomAatomBr0AB(38)

For Uvdw, 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:

UijCoulomb=qiqjrirj/4πε0rirj2(39)

where qi and qj are partial charges, ε0 is the dielectric permittivity of vacuum (8.85419 × 10−12 F/m).

Nonlinear spring element is used to represent interatomic force of Lennard-Jones potential. For Uvdw and Ucoul represented by nonlinear spring element, since transforming from xy1to xy1by the matrix [T], the expression is as:

xy1=xy1T=xy11,0,00,1,0r0,F0,0=x+r0y+F01(40)

where r0 and F0 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]:

UijVdWr=Urepuls if+Uattract if=ij4εijσijrij12σijrij6(41)

where rij 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, Rcut, usually set to Rcut = 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 rij. 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]:

[σ11σ22σ33σ12σ13σ23]=[c11c12c13000c11c13000c33000symc4400c4401/2×(c11c12)][ε11ε22ε33γ12γ13γ23]=[D1111D1122D1133000D2222D2233000D3333000symD121200D13130D2323][ε11ε22ε33γ12γ13γ23](42)

where σ represents the normal stress and shear stress in each direction (unit: nN/nm2); ε and γ are the normal strain and shear strain in each direction. 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 Cij 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=ε11,B=ε22,C=ε33,D=ε31,E=ε11+ε21(here εijindependent strain of each unit in all directions), as elastic constants Cij = Del, Eq. (42) can be changed as:

ABC00BAC0000A+BC000ED00000γ12c11c12c13c33c44=σ11σ22σ330σ12c11=D1111=D2222=Aσ11Bσ22A2B2CDσ33DA+B2CA+BEc12=D1122=Aσ22Bσ11A2B2CDσ33DA+B2CA+BEc13=D1133=D2233=Dσ33DA+BCE;c33=D3333=Eσ33DA+BCE;c44=D1212=D1313=σ12γ12;c66=12c11c12=D2323(43)

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

0=c11ε13+c12ε23+c13ε330=c12ε12+c11ε23+c13ε330=c13ε13+c13ε23+c33E3ε33c11c12c13c12c11c13c13c13c33E3=0(44)

Eq. (44) can be simplified as:

|c11c12c13c12c11c13c13c13c33E3|R2R1__|c11c12c13c12c11c11c120c13c13c33E3|C2+C1__|c11c11+c12c13c12c1100c132c13c33E3|=|c11+c12c132c13c33E3|=0.

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

E3=c332c132c11+c12(45)

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

E2=c11c121+c12c33c132c11c33c132(46)

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

E1=E2=c11c121+c12c33c132/c11c33c132(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 cij 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 ε1=εa,ε3=εc, so εbcan be easily obtained by the relative displacement in b-direction.

bcos60°sin60°0=b1/23/20,εb=1/23/20c11c12c13c12c11c13c13c13c331/23/20(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 (ε12=ε13=ε23=0), so

εb=14ε11+34ε22,ε22=43εb13ε11(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:

ECNTzigzag=43KC9C+4Kl02λz12+2λz22(50)

where K, C, respectively, correspond to constants of bond angle changes and bond stretching; λz1 and λz2 are structure parameters of zigzag nanotube, λz1=343cos2π2ncosπ2n/8323cos2π2n, and λz2=129cos2π2n/16343cos2π2n. 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:

ESLGSszigzag=83KC18C+Kl02(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=σε=F/AΔl/l=FlπDtΔl(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: rmax = 0.425 nm, Fmax = 7.0174 nN. When the distance r = 0.3833208 nm, the van der Waals force is zero, whereas the equilibrium distance r0 is 0.3415 nm, and the F0 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 θ0 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 ESLGSs and ECNT in Eq. (52), the results are as: ESLGSs = 0.867 TPa, ECNT = 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 (ESLGSs = 1.06 TPa) and literature results (ECNT = 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 {kr, , }. 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/m2, 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: kr = 669 nN/(nm2),  = 0.876 nN.nm, and  = 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]: Uvdw=15.2/r6+24.1×103/r12,r0=26σ=0.383nm. 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 θ0 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: h1 = 0.992σC-C = 0.3373 nm, h2-h1 = 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 anisotropyStrain-stress by beam-spring element
σij (nN/nm2)ε11ε22ε33σij (nN/nm2)ε11ε22ε33
σ11 = 902.320.879739−0.145212−0.301956σ11 = 957.7690.8797−0.15202−0.112044
σ22 = 889.883−0.1430920.8676238−0.297795σ22 = 944.594−0.149930.86762−0.112044
σ33 = 31.8732−0.0106694−0.01066940.882005σ33 = 15.641−5.278e-3−5.403e-30.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, Cij 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 = (Ex + Ey)/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.

ParametersC11C12C13C33C44C66BvBrGvGrEμ
FEM1057.21183.4915.0036.494.86436.84286.4335.75218.4811.00507.270.212
AFEM1123.02201.335.5035.722.56460.82300.7134.39231.156.07525.010.214
Equation-AFEM1123.10201.395.1238.242.61460.86300.8636.65231.406.20526.390.215
Equation-FEM1057.51183.8313.2945.244.86436.83286.7843.64219.3011.19511.530.217
Graphene [45]1060.00180.0015.0036.504.5440.00286.2835.76219.5710.26507.860.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 Cij 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 C33 larger than C11. 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 Ca2+ 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 O2− 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
SymbolWeight [u]Charges [e]BondDij (eV)Rij (nm)
Oh16.00−1.40O▬O0.155400.35532
Ho1.000.40O▬Ca3.83368 × 10−50.48980
Ca39.902.00Ca▬Ca2.18105 × 10−70.62428
Ohw16.00−0.95O▬H6.73854 × 10−30.35532
How1.000.42H▬H0.020000.20000

Table 5.

CH parameters of atomic charges and atomic force field employed (1e = 1.6021892 × 10−19 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 LO▬H = 0.986 nm, potential cutoff radius is 0.8 nm, EO▬H = 4.611 × 103 nN/nm2, dO▬H = 0.175 nm, GO▬H = 1.921 × 103 nN/nm2. 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 Cij 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 Cij 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 anisotropyStrain-stress by beam-spring element
σij (nN/nm2)ε11ε22ε33σij (nN/nm2)ε11ε22ε33
σ11 = 74.38320.834493−0.249651−0.118612σ11 = 41.01210.453152−0.064053−0.052899
σ22 = 59.7722−0.2006120.670576−0.095314σ22 = 32.3671−0.0571480.337453−0.052226
σ33 = 36.1151−0.057590−0.0575901.018542σ33 = 11.3929−0.041024−0.0410240.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.

MethodC11C12C13C33C44C66BvBrGvGrμ
FEM99.3830.787.3636.297.8834.3136.2326.6322.6513.920.2
AFEM94.0225.516.0232.317.09934.2532.8323.6721.8812.750.246
Equation-AFEM95.80625.2944.62635.6717.09835.25632.9324.5922.7412.990.243
Equation-FEM100.2231.6116.38937.4307.87834.30536.2926.7922.9114.000.255
Brillouin [56]102.0032.008.4033.6012.0035.0037.2426.0224.3918.180.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 Ex,y,z and cij of hexagonal crystal is derived by crystal elastic theory. Besides, a reliable equation to calculate these Cij 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 E1, y-axial modulus E2, and z-axial modulus E3 of hexagonal crystal as a function of elastic constants Cij are as:

      E3=c332c132c11+c12;E1=E2=c11c121+c12c33c132/c11c33c132

  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=ε11,B=ε22,C=ε33,D=ε31,E=ε11+ε21(εijare independent strains in all directions).

  • 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/m2. Young’s modulus of zigzag SWCNT is 1.05 Tpa and E.t is close to be 0.3604 kJ/m2. Besides, the calculated E of TLGSs model is close to the measurement of Berryman.

  • 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 cij are in an acceptable range with experimental values.

  • CH crystal has shown anisotropy. Based on the stress-strain data of AFEM modeling, elastic constant cij 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 σiis positioned R=Rii=1Nnand r=riσii=1Ne, respectively.

    Rr=Rr(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 wi, the total Helmholtz free energy is:

    F=E+iFith=U+i12wi+kBTiIn1ewikBT(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]:

    Ĥ=2meii2+iVexri+12ije2rirj(3)
    2m2+Veffrφir=εiφirVeffr=Vr+Veer+VexrVxer=δExeρδrρr=iϕir2(4)

    where Vex(r) is an external potential; ri and rj are the nucleus position vector; while , 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:

    ρ0FηijT=ρ0FηijT+12ijklcijklTηijηkl++1n!ijklcijklTηijηkl(5)

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

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

    ΔE=V2i=16j=16cijeiej(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.

    ValueCoarseMediumFineUltrafine
    SCF tolerance (eV/atom)1.0 × 10−52.0 × 10−61.0 × 10−65.0 × 10−7
    k-Points separation (1/Å)0.070.050.040.04
    Energy cutoffCoarseMedium340 eV380 eV
    Energy tolerance (eV/atom)5.0 × 10−52.0 × 10−51.0 × 10−55.0 × 10−6
    Max. force tolerance (eV/Å)0.10.050.030.01
    Max. stress tolerance (GPa)0.20.10.050.02
    Max. displacement tolerance (Å)5.0 × 10−32.0 × 10−31.0 × 10−35.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: ψRr=χRφRrand assume that the electron mobility (φ) 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.

    Ĥ=T̂e+Ûne+Ûee+Ûnn(7)
    ĤRφR0r=EBORφR0rr(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:

    EBORφR0ĤφR0=minφR0ĤφR0(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:

    ĤBO=T̂n+EBOR(10)
    ĤBOχR=R(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 Vexr, where the Hamiltonian is written as:

    Ĥ=2meii2+iVexri+12ije2rirj(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 Vexriand 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.

    EHKn=Tn+Eintn+d3rVexr+EnmRFHKn+d3rVexr+EnnR(13)

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

    EBOR=minERnr(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 nrd3r=Ne. In this stage, the DFT can reformulate the problem rather than solve an uncertain functional FHKn.

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

    EKS=Fn+d3rVexr=TSn+EHn+EXCn+d3rVexr(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,

    TSn=ψNIT̂eψNI=i=1Neφi122φi(16)

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

    EHartreen=12d3rd3rnrnrrr(17)
    nr=i=1Neφir2(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 TS 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:

    δEKSδφir=δTSδφir+δEexδnr+δEHartreeδnr+δExcδnrδnrδφir=0(19)

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

    ĤKSεiφir=0(20)

    εirepresents the eigenvalues, and ĤKSis the effective Hamiltonian H,

    ĤKSr=122+VKSr(21)
    VKSr=Vexr+δEHartreeδnr+δExcδnr(22)

    Eqs. (20)–(22) are known equations of Kohn Sham, the density n(r), and the resulting total energy EKS. These equations are independent of any approximation on the functional EXCn, resolution provides the exact values of the density and the energy of the ground state of the interacting system, provided that EXCnis known. The latter can be described in terms of Hohenberg Kohn function (12)

    Excn=FHKnTSn+EHartreen(23)
    Excn=T̂TSn+V̂intEHartreen(24)

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

    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 c11c12c12000c12c11c12000c12c12c11000000c44000000c44000000c44, 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 δ.

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

    e=e1e2e3e4e5e6(25)

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

    ΔE=V2(c11e1e1+c11e2e2+c11e3e3+c12e1e2+c12e1e3+c12e2e1+c12e2e3+c12e3e1+c12e3e2+c44e4e4+c44e5e5+c44e6e6)(26)

    Firstly, considering the simplest constant C44, i and j are equal to 4, and by the cubic matrix elements, C44 = C55 = C66. By setting e=000δδδof Eq. (25) and submitting into Eq. (5), we get:

    ΔE=V02c44e4e4+c44e5e5+c44e6e6(27)

    Therefore,

    ΔEV0=32c44δ2(28)

    Similarly, setting e=δδ0000, we obtain:

    ΔE=V02c11e1e1+c11e2e2+c12e1e2+c12e2e1(29)

    Therefore,

    ΔEV0=c11+c12δ2(30)

    On this basis, if a strain e=δδδ000is imposed, we have:

    ΔE=V02c11e1e1+c11e2e2+c11e3e3+c12e1e2+c12e1e3+c12e2e1+c12e2e3+c12e3e1+c12e3e2(31)

    Therefore,

    ΔEV0=32c11+2c12δ2(32)

    After fitting the three equations, we can obtain three independent elastic constants c11, c12, and c44.

    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 δ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 tensorLCEC
    e=δδ00002c11 + 2c12
    e=δδ00002c11−2c12
    e=00000δ3c11 + 6c12
    e=δδδ0003c11 + 6c12
    e=δ2/1δ2002δ004c44
    e=000δδδ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]:

    C11>0,C44>0,C11>C12,C11+2C12>0(33)

    4.2 Elastic strain energy formulas of hexagonal crystal structures

    Matrix element of hexagonal crystal structure is as: c11c12c13000c12c11c13000c13c13c33000000c44000000c4400000012c11c12, system has a total of five independent components, of which independent elastic constants are c11, c12, c13, c33, and c44.

    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:

    ρ0F(ηJ,T)=ρ0F(ηJ,T)+12{c11Tη1η1+c11Tη2η2+c33Tη3η3+c44Tη4η4+c44Tη5η5+c11Tc12T2η6η6+2c12Tη1η2+2c13Tη1η3+2c13Tη2η3}++1n!JPcJPTηJηP(34)

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

    R=a1a2a3=32a12a032a12a000c(35)

    where a and c are the lattice constants of crystals.

    We can apply strain e=δδ0000to calculate the value of c11 + c12:

    ΔEV0=c11+c12δ2(36)

    Similarly, we can apply strain e=00000δto calculate the value of c11−c12:

    ΔEV0=14c11c12δ2(37)

    Also for c33, we can apply strain e=00δ000:

    ΔEV0=12c33δ2(38)

    Then, for c44, we can apply strain e=000δδ0:

    ΔEV0=c44δ2(39)

    Finally, we can apply strain e=δδδ000to calculate the values of c11, c12, c13, and c33:

    ΔEV0=c11+c12+2c13+12c33δ2(40)

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

    c44>0,c11>c12,c11+2c12c33>2c132(41)

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

    Distortion matrixEnergy (Taylor’s expansion)1V02E(V,δ)δ2|δ=0
    1+δ0001+δ0001EVδ=EV00+V0τ1+τ2δ+c11+c12δ22c11+c12
    1+δ0001δ0001EVδ=EV00+V0τ1τ2δ+c11c12δ22c11c12
    100010001+δEVδ=EV00+V0τ3δ+1/2c33δ2c33
    10001δ0δ1EVδ=EV00+V0τ5δ+2c44δ24c44
    1+δ00010001+δEVδ=EV00+V0τ1+τ3δ+1/2c11+2c13+c33δ2c11+2c13+c33
    1+δ00010001δEVδ=EV00+V0τ1+τ3δ+1/2c112c13+c33δ2c112c13+c33
    1+δ0001+δ0001+δEVδ=EV00+V0τ1+τ2+τ3δ+ 1/22c11+2c12+4c13+c33δ22c11+2c12+c33+4c13

    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 Cij 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 B1–B2 transition at 61 GPa [17], and phase transformation from B1 to B2 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.

    AtomxyzOccupancy rateUiso or Ueq
    Ca0.0000.0000.0001.001.00
    O0.5000.5000.5001.001.00

    Table 4.

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

    The elastic properties of the B1-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 B1-CaO at zero pressure are given. In addition, the elastic constants, shear modulus, and bulk modulus of the B1-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. Mg2+ ions will fill into the O−2 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.

    AtomxyzOccupancy rateUiso or Ueq
    Mg0.0000.0000.0001.001.00
    O0.5000.5000.5001.001.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−6 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−6 eV per atom.

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

    aC11C12C44BG
    Present (LDA)4.7124289.49771.43276.652144.1288.295
    Present (GGA)4.8339200.87956.90275.600104.8974.134
    Ref. [21]4.722745453128
    Ref. [22]4.820206506610271
    Ref. [23]4.83823951.677.411783.6
    Ref. [24]4.840223538410985
    Exp. [37]221.8957.8180.3280.59
    Exp. [38]220.5357.6780.0380.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.

    PaV(Å3)c11c12c44
    0.0 (LDA)4.7124104.647289.49771.43276.652
    0.0 (GGA)4.8339112.951200.87956.90275.600
    0.1 (GGA)4.8324112.847201.54856.74275.621
    0.2 (GGA)4.8310112.748202.53856.98675.639
    0.3 (GGA)4.8299112.665203.06256.89375.659
    0.4 (GGA)4.8279112.534204.50657.36175.683
    0.5 (GGA)4.8264112.428205.66757.73675.717
    0.6 (GGA)4.8251112.336206.10557.44075.717
    0.7 (GGA)4.8238112.249207.00157.66475.753
    0.8 (GGA)4.8219112.110208.20257.76275.784
    0.9 (GGA)4.8204112.008209.41558.19275.890
    1.0 (GGA)4.8193111.934209.83058.05975.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, Cij, B, and G of the B1-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.

    aC11C12C44BGA
    Present (LDA)4.2130298.50892.228141.197160.988124.506
    Present (GGA)4.2130343.468106.358140.823185.395131.449
    Ref. [19]4.2130167.600
    Ref. [22]4.30270731271391150.21
    Ref. [24]4.259031887144164.0001320.18
    Ref. [29]298.9696.42157.13163.930131.80
    Ref. [33]4.221030093.6147162.4127.60.29
    Ref. [35]4.251029191139158.00121.8000.27
    Ref. [39]4.16–4.24286–35291–108158–188157.0–198.0
    Ref. [40]4.212

    Table 8.

    Lattice a, elastic constants Cij, 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 cij of CaO at zero temperature and different pressures. It is found that c11 varies largely under the effect of pressure, compared with the variations in c12 and c44. The change of elastic constant c11 represents a produced change with a longitudinal strain in length elasticity. The elastic constants c12 and c44 are related to the elasticity in shape, which is a shear constant. Besides, c12 and c44 are less sensitive of pressure as compared with c11. As pressure increases, c11 and c12 increase monotonically, but c44 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−6 eV per atom.

    As can be seen in Table 8, the obtained a, Cij, 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 c12 and c44 for MgO are positive indicating that the MgO is stable at zero pressure (Table 9).

    PaV3)c11c12c44
    0.0 (GGA)4.213074.778343.468106.358140.823
    0.0 (LDA)4.213074.778298.50892.228141.197
    0.1 (GGA)4.203874.289309.89295.878142.181
    0.2 (GGA)4.203074.246308.98994.200142.178
    0.3 (GGA)4.202274.202310.99495.635142.300
    0.4 (GGA)4.201274.153312.88596.635142.371
    0.5 (GGA)4.200574.115313.05096.187142.484
    0.6 (GGA)4.199674.069313.56596.092142.519
    0.7 (GGA)4.198874.025314.76696.493142.597
    0.8 (GGA)4.197973.980314.61295.599142.677
    0.9 (GGA)4.197173.937316.35396.652142.857
    1.0 (GGA)4.196373.893317.06796.738142.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 c11, c12, and c44. 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 c11 varies largely under the effect of pressure as compared with the variations in c12 and c44.

    5.2 Elastic coefficients Cij 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, Ca2+ 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. Ca2+ is octahedrally coordinated by oxygen and each O2− 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 O2− and H+ is 2.757 Å (Table 10).

    Atomic speciesXYZOccupancy rateUiso or Ueq
    Ca0.00000.00000.00001.001.00
    O0.33330.66670.23461.001.00
    H0.33330.66670.42801.001.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) CO3 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 (CO3−2), form together an equilateral triangle, including the central carbon atom and each vertex totalizing three oxygen atoms. (2) CO3 group interatomic distances are as: C▬O is 1.2815 Å, Ca▬O1 is 3.4589 Å, O▬O1 is 2.2197 Å, O▬O2 is 4.7911 Å, O▬O3 is 3.4107 Å, and O▬O4 is 4.2459 Å (Table 11).

    Atomic speciesXYZOccupancy rateUiso or Ueq
    Ca0.00000.00000.00001.001.00
    C0.00000.00000.25001.001.00
    O0.256420.00000.25001.001.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−6 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.

    Pressurea0(Å)c0(Å)V03)C11C12C13C14C33C44C66
    Cal. [43]03.59304.909054.883099.3930.787.36036.297.8834.305
    Exp. [44]03.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]03.65804.793055.5000135.145.66.6026.04.144.7
    Exp. [59]03.58904.911054.800099.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)
    Present0.0 (GGA)3.59184.906354.8162103.94830.9232.5644.33642.67014.14036.513
    0.0 (LDA)3.62335.008656.944595.41727.1621.3130.24227.5858.99834.127
    0.5 (GGA)3.61184.938355.789997.14027.0429.5530.50335.89912.00335.049
    1.0 (GGA)3.59994.8671654.6274101.46333.61310.1932.72554.18210.90433.925
    1.5 (GGA)3.59034.820753.8137107.68533.88912.5933.39650.86712.04536.898
    2.0 (GGA)3.58154.785553.1590107.23634.68117.6500.28549.07220.81836.278
    2.5 (GGA)3.57334.75552.5737112.24937.19017.1392.66158.66619.95937.529
    3.0 (GGA)3.56514.724552.0042114.38839.66422.1112.37264.58216.34037.362
    3.5 (GGA)3.55724.696151.4622117.50541.66126.4852.84665.96016.47637.922
    4.0 (GGA)3.54924.671850.9649119.69644.86926.4980.64873.28517.54337.414
    4.5 (GGA)3.54154.646150.4651120.40044.26125.0892.85173.88422.54338.070
    5.0 (GGA)3.53424.625850.0389128.86549.44536.7024.78183.20320.34639.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: c11 = 103.948 GPa, c12 = 30.923 GPa, c13 = 2.564 GPa, c33 = 42.670 GPa, and c44 = 14.140 GPa. C14 is turned out not to be zero, which is different to the result of Laugesen. We can see that all elastic coefficients (c14 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−6 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: c11 = 143.938 GPa, c12 = 55.579 GPa, c13 = 50.970 GPa, c33 = 81.175 GPa, and c44 = 32.658 GPa. The elastic coefficients of calcite under pressure 0–5 GPa show slight growing trends except for c11 and c12 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/a0V03)c11c12c13c14c33c44c66
    Ref. [46]146.2659.7050.76−20.7685.3134.0543.28
    Ref. [48]4.98917.044367.39
    Ref. [49]4.98115.902341.676
    Ref. [51]4.99017.061367.8
    Ref. [54]4.98817.068367.761159706319.5983944.5
    Ref. [55]149.457.953.52085.234.145.75
    Ref. [56]145.755.953.520.585.333.444.9
    Ref. [8]5.06117.097379.279152.357.0554.8317.1487.7736.0047.625
    0.0 (LDA)4.971816.52963.3247353.526153.926155.806450.969813.654182.701533.945159.9905
    0.0 (GGA)5.054817.25513.4136381.818143.938355.579150.969816.867881.175332.658344.1796
    0.5 (GGA)5.047817.18083.4036379.127148.47957.757453.496717.685583.256533.688745.3608
    1.0 (GGA)5.041117.10893.3939376.536150.162960.223255.528417.357384.800534.597144.9699
    1.5 (GGA)5.035717.03383.3826374.074152.058460.937758.00817.841287.041235.065245.5604
    2.0 (GGA)5.029316.96423.3731371.609159.588663.108160.753317.881788.25334.694248.2403
    2.5 (GGA)5.024516.89393.3623369.350158.223975.89661.965919.568790.031335.531941.1640
    3.0 (GGA)5.019416.82423.3518367.090164.144968.867365.108521.300591.426636.766947.6388
    3.5 (GGA)5.014416.75693.3418364.885167.898472.126567.452321.418492.962537.703647.8860
    4.0 (GGA)5.010116.68873.3310362.789170.100974.842869.61821.453494.68337.322447.6291
    4.5 (GGA)5.004916.62853.3224360.724172.889576.984472.447522.881496.453937.687147.9526
    5.0 (GGA)5.000916.56203.3118358.704177.115778.618674.845423.657197.655137.54849.2486

    Table 13.

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

    From Table 13, elastic constants have been calculated by GGA method as follows: c11 = 143.938 GPa, c12 = 55.579 GPa, c13 = 50.970 Pa, c33 = 81.175 GPa, and c44 = 32.658 GPa. By LDA method, elastic constants have been calculated as follows: c11 = 153.926 GPa, c12 = 55.806 GPa, c13 = 50.970 Pa, c33 = 82.702 GPa, and c44 = 59.991 GPa.

    Similarly, for elastic constants Cij 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: c11 = 200.879 GPa, c12 = 56.902 GPa, and c44 = 75.600 GPa.

      2. Elastic constants of our MgO single crystal are as: c11 = 298.508 GPa, c12 = 92.228 GPa, and c44 = 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: c11 = 103.948 GPa, c12 = 30.923 GPa, c13 = 2.564 GPa, c14 = 4.336 GPa, c33 = 42.670 GPa, and c44 = 14.140 GPa. Besides, elastic constants based on LDA method are as: c11 = 95.417 GPa, c12 = 27.162 GPa, c13 = 1.313 GPa, c14 = 0.242 GPa, c33 = 27.585 GPa, and c44 = 8.998 GPa. We can see that C14 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: c11 = 143.938 GPa, c12 = 55.579 GPa, c13 = 50.970 GPa, c33 = 81.175 GPa, and c44 = 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 CaCO3) 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 CaCO3. 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 SHill=12SReuss+CVoigt1and CHill=12CVoigt+SReuss1, it can be seen that the stiffness coefficient C and the compliance coefficient S do not satisfy the relation CS=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 L3 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 L3 within system [S] relative to l, m, n, one can get [3]:

    l=sinψcosϕm=sinψsinϕn=cosψ(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 εϕψat the direction L3 is as [22, 23]:

    εϕψ=l2ε11+m2ε22+n2ε33+lmγ12+mnγ23+nlγ31(2)

    Therefore, whole stress σϕψat the direction L3 is as [22, 23]:

    σϕψ=l2σ11+m2σ22+n2σ33+2lmτ12+mnτ23+nlτ31(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: λ=2μG12μ, G=E21+μ; here, E, μ, 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 τ12=τ23=τ31=0; then, we can see that γ12=γ23=γ31=0. If suppose ϕ=0°, Eq. (3) can be simplified as [3]:

    σϕψ=l2σ11+m2σ22+n2σ33=λε11+ε22+ε33+2Gε11ε33sin2ψ+ε33(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., σ110, other stress components are zero) and under the conditions of ϕ=0°, the ratio expression between the whole strain and the uniaxial stress from Eq. (4) will be eventually developed as [3]:

    εϕψσ11=1+μEsin2ψμE(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 C11, C12, and C44 can be calculated as follows:

    E100=c112+c11c122c122c11+c12(6)
    E110=4c112+c11c122c122c44c112+c11c12+2c11c442c122(7)
    E111=3c11+2c12c44c11+2c12+c44(8)

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

    GYuvw=12s11s1232s11s12s44Γ(9)
    EYuvw=1s112s11s12s44Γ(10)
    μYuvw=2s12+2s11s12s44Γ2s112s11s12s44Γ(11)

    where Γ¯=u2v2+v2w2+w2u2. 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 L3 (u, v, w) relative to the cubic crystal coordinates [K], the relationship between the macroscopic stress σL3 in L3 direction and the stress component (σx, σy, σz, τyz, τzx, τxy) of the coordinate system [K] is as [22, 23]:

    σL3=u2σx+v2σy+w2σz+2vwτyz+wuτzx+uvτxy(12)

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

    σxσyσzτyzτzxτxy=c11c12c12000c12c11c12000c12c12c11000000c44000000c44000000c44εxεyεzγyzγzxγxy(13)

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

    σL3=c11u2εx+v2εy+w2εz+c12v2+w2εx+u2+w2εy+u2+v2εz+2c44vwγyz+uwγzx+uvγxy(14)

    With respect to the coordinate system [K], assuming the three direction cosines of coordinate axes [S] (respectively S1, S2, S3) are S1 (l1, m1, n1), S2 (l2, m2, n2), S3 (l3, m3, n3), and considering relationship between strain components in coordinates system [S] and strain component of coordinates system [K], we can get relation γ12=γ23=γ31=0when the polycrystalline material is limited only by the principal stress (τ12=τ23=τ31=0). So the following formula of strain holds [22]:

    εx=l12ε11+l22ε22+l32ε33εy=m12ε11+m22ε22+m32ε33εz=n12ε11+n22ε22+n32ε33γyz=2m1n1ε11+m2n2ε22+m3n3ε33γzx=2n1l1ε11+n2l2ε22+n3l3ε33γxy=2l1m1ε11+l2m2ε22+l3m3ε33(15)

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

    σL3=c11[u2(l12ε11+l22ε22+l32ε33)+v2(m12ε11+m22ε22+m32ε33)+w2(n12ε11+n22ε22+n32ε33)]+c12[(v2+w2)(l12ε11+l22ε22+l32ε33)+(u2+w2)(m12ε11+m22ε22+m32ε33)+(u2+v2)(n12ε11+n22ε22+n32ε33)]+4c44[vw(m1n1ε11+m2n2ε22+m3n3ε33)+uw(n1l1ε11+n2l2ε22+n3l3ε33)+uv(l1m1ε11+l2m2ε22+l3m3ε33)](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 σL3¯so as to determine its macroscopic stress. Here, if ϕ=0°, we can obtain the mean stress σYL3¯by averaging the orientation in the crystal plane:

    σYL3¯=02πσL302π(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):

    σYL3¯=c11c123c11c122c44Γε11ε33sin2ψ+c11c123c11c122c44Γε33+ε11+ε22+ε33c12+c11c122c44Γ(18)

    where Γ=u2v2+v2w2+u2w2,uvwis 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]:

    GYuvw=c11c1223c11c122c44Γ2(19)
    λYuvw=c12+c11c122c44Γ(20)

    As λ=2μG12μ, G=E21+μ, formulas can be developed as:

    1+μEY=1c11c12c11c122c443Γ(21)
    μEY=c12+c11c122c44Γc11+2c12c11c123c11c122c44Γ(22)
    EYuvw=c11+2c12c11c123c11c122c44Γc11+c12c11c122c44Γ(23)
    μYuvw=c12+c11c122c44Γc11+c12c11c122c44Γ(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 Γ¯=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]:

    GM=c11c12+3c445,λM=c11+4c122c445(25)
    1+μEM=52c11c12+3c44(26)
    μEM=c11+4c122c442c11+2c12c11c12+3c44(27)
    EM=c11+2c12c11c12+3c442c11+3c12+c44(28)
    μM=c11+4c122c4422c11+3c12+c44(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., σ110, other stress components are zero) and under the conditions of ϕ=0°, the ratio expression between the whole strain and the uniaxial stress from Eq. (4) will be eventually developed as [4]:

    σϕψ=λε11+ε22+ε33+2Gε11ε33sin2ψ+ε33(30)

    3.1 Reuss bounds of hexagonal polycrystal structure

    As in Figure 1, supposing that the direction cosines are L3 (u, v, w) relative to the hexagonal crystal coordinates [K], the relationship between the macroscopic strain εL3(equivalent to εϕψ) in L3 direction and the strain component (εx, εy, εz, γxy, γyz, γzx) of the coordinate system [K] is as [22, 23]:

    εL3=u2εx+v2εy+w2εz+uvγxy+vwγyz+wuγzx(31)

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

    εxεyεzγyzγzxγxy=s11s12s13000s12s11s13000s13s13s33000000s44000000s440000002s11s12σxσyσzτyzτzxτxy(32)

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

    εL3=s11u2σx+v2σy+2uvτxy+s12u2σy+v2σx2uvτxy+s13u2+v2σz+w2σ+xσy+s33w2σz+s44wvτyz+uτzx(33)

    With respect to the coordinate system [K], assuming the three direction cosines of coordinate axes [S] (respectively S1, S2, S3) are S1 (l1, m1, n1), S2 (l2, m2, n2), S3 (l3, m3, n3), the relationship between the stress components in the coordinate system [K] and the stress component of the coordinate system [S] is as follows [26]:

    σx=liljσij(34)
    σy=mimjσij(35)
    σz=ninjσij(36)
    τxy=limjσij(37)
    τyz=minjσij(38)
    τzx=niljσij(39)

    When the material is only under the uniaxial stress σ110and other stress components are zero, under these conditions, substituting Eqs. (34)–(39) into Eq. (54), one can get:

    εL3=s11u2l12+v2m12+2uvl1m1+s12u2m12+v2l122uvl1m1+s13u2+v2n12+w2l12+m12+s33w2n12+s44wvm1n1+ul1n1σ11(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 εL3¯so as to determine its macroscopic strain. Here, if ϕ=0°, we can obtain the mean strain εYL3¯by averaging the orientation in the crystal plane:

    εYL3¯=02πεL302π(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):

    εYL3¯σ11=122s11s12s13w225s11s125s13+s333s44+3w42s112s13+s33s44×sin2ψ+12s12+s13+w22s11s12s13+s33s44w42s112s13+s33s44(42)

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

    1+μEY=122s11s12s13w225s11s125s13+s333s44+3w42s112s13+s33s44(43)
    μEY=12s12+s13+w22s11s12s13+s33s44w42s112s13+s33s44(44)

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

    EYuvw=1s112s11s13s44w2+s112s13+s33s44w4(45)
    μYuvw=s12+s13+s11s12s13+s33s44w2s112s13+s33s44w42s112s11s13s44w2+s112s13+s33s44w4(46)
    GYuvw=12s11s12s135s11s125s13+s333s44w2+3s112s13+s33s44w4(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 w2¯=13and w4¯=15[25], the mechanical properties of polycrystalline, such as Young’s modulus E, Poisson’s ratio μ, and the shear modulus G, are as [4]:

    1+μEM=7s115s124s13+2s33+3s4415(48)
    μEM=s11+5s12+8s13+s33s4415(49)
    E=1524s11+2s13+s44+3s33(50)
    μ=s11+5s12+8s13+s33s4424s11+2s13+s44+3s33(51)
    G=1527s115s124s13+2s33+3s44(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:

    σL3=c11u2εx+v2εy+uvγxy+c12u2εy+v2εxuvγxy+c13u2+v2εz+w2ε+xεy+c33w2εz+2c44wvγyz+uγzx(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:

    σYL3¯=122c11c12c13w25c11c125c13+c3312c44+3w4c112c13+c334c44ε11ε33sin2ψ+122c11c12c13w25c11c125c13+c3312c44+3w4c112c13+c334c44ε33+12ε11+ε22+ε33c12+c13+w2c11c12c13+c334c44w4c112c13+c334c44(54)

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

    GYuvw=142c11c12c13w25c11c125c13+c3312c44+3w4c112c13+c334c44(55)
    λYuvw=14c12+c13+w2c11c12c13+c334c44w4c112c13+c334c44(56)

    As G=E21+μ, λ=2μG12μ, formulas can be developed as follows [4]:

    1+μEY=22c11c12c13w25c11c125c13+c3312c44+3w4c112c13+c334c44(57)
    μEY=c12+c13+w2c11c12c13+c334c44w4c112c13+c334c44c11+c12+c13w2c11+c12c13c332c11c12c13w25c11c125c13+c3312c44+3w4c112c13+c334c44(58)
    EYuvw=c11+c12+c13w2c11+c12c13c332c11c12c13w25c11c125c13+c3312c44+3w4c112c13+c334c442c11+c12+c13w23c11+c123c13c334c44+w4c112c13+c334c44(59)
    μYuvw=c12+c13+w2c11c12c13+c334c44w4c112c13+c334c442c11+c12+c13w23c11+c123c13c334c44+w4c112c13+c334c44(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 w2¯=13and w4¯=15[25], the mechanical properties of polycrystalline, such as Young’s modulus E, Poisson’s ratio μ, and shear modulus G, are as [4]:

    GM=7c115c124c13+2c33+12c4430(61)
    λM=c11+5c12+8c13+c334c4415(62)
    1+μEM=157c115c124c13+2c33+12c44(63)
    μEM=3c11+5c12+8c13+c334c442c11+c12+2c13+c337c115c124c13+2c33+12c44(64)
    E=2c11+c12+2c13+c337c115c124c13+2c33+12c4439c11+5c12+43c13+c33+c44(65)
    μ=c11+5c12+8c13+c334c449c11+5c12+43c13+c33+c44(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 and shear modulus ) 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 in either direction of crystal c-axis with the angle Ф can be described as [1]:

    1Eϕ=S11l14+2S12l12l22+2S13l12l32+2S15l13l3+S22l24+2S23l22l32+2S25l1l22l3+S33l34+2S35l1l33+S44l22l32+2S46l1l22l3+S55l12l32+S66l12l22(67)

    where l1l2l3=sinθcosϕsinθsinϕcosθ, Sij is the single crystal compliance tensor in collapsed matrix notation and l1, l2, and l3 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 Ehkl 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 Cij and the elastic compliance constants Sij 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: c11 = 200.879 GPa, c12 = 56.902 GPa, and c44 = 75.600 GPa. By LDA, the Cij are as: c11 = 289.497 GPa, c12 = 71.432 GPa, and c44 = 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.

    PBvBrGvGrBGEμ
    Ref. [27]127.3333127.333375.800066.8578127.333371.3289180.31700.2640
    Ref. [28]10210270.800070.327910270.5639172.02310.2189
    Ref. [29]114.0667114.066783.920083.1886114.066783.5543201.47030.2056
    Ref. [30]109.6667109.666784.400084.3972109.666784.3986201.50380.1938
    0 (GGA)104.894104.89474.15574.113104.89474.134179.99780.2140
    0 (LDA)144.120144.12089.60486.985144.12088.295219.96400.2456
    0.1 (GGA)105.011105.01174.33474.300105.01174.317176.6170.2197
    0.2 (GGA)105.503105.50374.49474.467105.50374.481177.5120.2196
    0.3 (GGA)105.616105.61674.62974.608105.61674.619178.1590.2189
    0.4 (GGA)106.409106.40974.83974.825106.40974.832179.3770.2190
    0.5 (GGA)107.046107.04675.01675.006107.04675.011180.3570.2192
    0.6 (GGA)106.995106.99575.16375.157106.99575.160181.0670.2180
    0.7 (GGA)107.443107.44375.31975.315107.44375.317181.8730.2179
    0.8 (GGA)107.908107.90875.55875.557107.90875.558183.1130.2172
    0.9 (GGA)108.600108.60075.73075.730108.60075.730184.1080.2175
    1.0 (GGA)108.650108.65075.84875.848108.65075.848184.6640.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 Cij. 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.

    PBvBrGvGrBGEμ
    Ref. [28]138.6667138.6667115.6000113.8262138.6667114.7131269.75390.1758
    Ref. [30]164.0000164.0000132.6000131.0638164.000131.8319311.91720.1830
    Ref. [31]163.9333163.9333134.7860128.7278163.9333131.7569311.75040.1831
    Ref. [32]162.4000162.4000129.4800125.6660162.4000127.5730303.30010.1887
    0 (GGA)185.395185.395131.9160130.9822185.395131.449293.1720.2364
    0 (LDA)160.988160.988125.974123.037160.988124.506254.9690.2360
    0.1 (GGA)167.217167.217128.111125.659167.217126.885264.5820.2363
    0.2 (GGA)165.796165.796128.265125.871165.796127.068264.9720.2336
    0.3 (GGA)167.421167.421128.452126.085167.421127.268266.0090.2352
    0.4 (GGA)168.718168.718128.673126.362168.718127.517267.2790.2360
    0.5 (GGA)168.474168.474128.863126.583168.474127.723267.8350.2350
    0.6 (GGA)168.583168.583129.006126.766168.583127.886268.4840.2346
    0.7 (GGA)169.250169.250129.213127.020169.250128.116269.4860.2346
    0.8 (GGA)168.603168.603129.409127.258168.603128.334270.0540.2330
    0.9 (GGA)169.886169.886129.654127.530169.886128.592271.1160.2340
    1.0 (GGA)170.181170.181129.830127.738170.181128.784271.8360.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 cij, 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: c11 = 343.468 GPa, c12 = 106.358 GPa, and c44 = 140.823 GPa. By LDA, the Cij are as: c11 = 298.508 GPa, c12 = 92.228 GPa, and c44 = 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.230026.631622.651013.919631.430818.285345.94600.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.977822.351326.41678.441334.164617.429044.68780.2820
    Ref. [36]40.242228.064821.295915.303234.153518.299646.57960.2727
    0 (GGA)36.639528.182727.023321.079332.411124.051358.07940.2074
    0 (LDA)30.888519.582622.999914.890025.235518.945045.459170.199767
    0.5 (LDA)35.830827.100224.079618.241031.465521.160351.856530.225326
    1.0 (LDA)40.567335.086524.687117.757737.826921.222453.636470.263676
    1.5 (LDA)42.709835.681926.008318.986639.195922.497556.653290.259102
    2.0 (LDA)44.834137.419228.486825.820941.126727.153966.767330.229424
    2.5 (LDA)47.344341.265729.602626.358844.305027.980769.344080.239141
    3.0 (LDA)51.236646.060127.973123.603448.648425.788265.747210.274754
    3.5 (LDA)54.470149.126927.930623.546351.798525.738466.243230.286856
    4.0 (LDA)56.489751.943728.820825.078354.216726.949569.356770.286791
    4.5 (LDA)55.951551.430131.314328.837953.690830.076176.031400.263984
    5.0 (LDA)65.181161.350930.619326.997463.266028.808475.035910.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: Gv = 30.889 GPa, Bv = 22.999 GPa, Gr = 19.583 GPa, Br = 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: c11 = 143.938 GPa, c12 = 55.579 GPa, c13 = 50.970 Pa, c33 = 81.175 GPa, and c44 = 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/BE (GPa)μ
    Ref. [37]80.643374.663238.969036.578877.653337.77390.48644397.51060.2907
    Ref. [38]77.807871.552936.716735.076874.680335.89670.48067292.81850.2929
    Ref. [39]89.777883.838239.166737.491186.808038.32890.441536100.23430.3076
    Ref. [40]79.311172.922837.396735.154376.116936.27550.47657693.90830.2944
    Ref. [41]78.055672.515736.593334.423975.285635.50860.47165292.05350.2962
    0.0 (LDA)78.449570.942538.910536.119974.696037.51520.502296.40610.285
    0.0 (GGA)76.009969.625536.001428.713772.817732.35760.444484.54920.306
    0.5 (GGA)78.857372.099636.911929.168675.478533.04020.437786.49920.309
    1.0 (GGA)80.854073.945537.089229.691077.399833.39010.431487.57690.311
    1.5 (GGA)82.785076.168037.418429.742679.476533.58050.422588.30470.315
    2.0 (GGA)86.295578.569438.380130.357582.432534.36880.416990.52530.317
    2.5 (GGA)89.570580.559536.222327.445585.06531.83390.374084.90970.334
    3.0 (GGA)90.876182.514138.943328.970886.695133.95710.391790.10680.327
    3.5 (GGA)93.646884.626839.440529.475889.136834.45820.386691.57430.329
    4.0 (GGA)95.893686.620139.175229.138691.256834.15690.374391.10410.334
    4.5 (GGA)98.443588.920039.335628.272393.681833.81390.360990.54760.306
    5.0 (GGA)100.945090.794039.774127.995295.869533.88460.353490.93980.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: Gv = 36.001 GPa, Bv = 76.010 GPa, Gr = 28.714 GPa, Br = 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 Cij 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 Cij (or the compliance coefficient Sij) 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 Γ¯=1/5or w2 = 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: Gv = 30.889 GPa, Bv = 23.000 GPa, Gr = 19.583 GPa, Br = 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: Gv = 27.023 GPa, Bv = 36.640 GPa, Gr = 21.079 GPa, Br = 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: Gv = 36.001 GPa, Bv = 76.010 GPa, Gr = 28.714 GPa, Br = 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 nm3 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 C1.7▬S▬H1.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 [SiO4]4− chain under the system equilibrium structure C1.65▬S▬H1.75, comparing with experimental values of C1.7▬S▬H1.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=12i=1Npi2mi+i=1N1j=i+1NUrij(1)
    pi=midridt=mivi(2)
    dpidt=i=1N1j=i+1NFrij=i=1N1j=i+1NUrijrij(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]:

    rit+Δt=rit+vit+Δt/2Δt(4)
    vit+Δt/2=vitΔt/2+d2ritdt2Δt+ΟΔt3(5)
    d2ritdt2=Fimi=1miUr(6)

    where F is the net force vector on the ith molecule, m is the mass, and r is the position vector of the ith particle. For Verlet algorithm, particle acceleration: ait=Fitmi. Initial conditions of Verlet algorithm are as: rit=0=ri0, dridtt=0=vi0.

    There are a few methods to control temperature, the simplest method involves velocity scaling or coupling with heat bath, and other more complicated approaches are the Andersen and Nose-Hoover thermostats [45]. In this research, a special form of Nose-Hoover thermostats developed by Berendsen is used, in which the friction coefficient λ rather than its time derivative varies according to the following equation rather than its derivative [46]:

    λ=i=1Npi22migkT2/Q(7)

    where Q = gkTτ2, g is the degree of freedom in the system (g = 3N − 4) for a thermostated monatomic system with zero momenta, τ is thermostat relaxation time typically about 0.5.

    From Figure 1, the domain reduction only performs analysis of a representative substructure. During MD simulation, only one box is modeled explicitly. Each particle interacts with other particles in the box and with their images in nearby boxes. The mean square displacement (MSD) is a measure of the average distance a molecule travels, it represents the statistical average of the particle trajectories change over time. By solving differential motion equations, which is as [44]:

    MSDt=1ni=1nrntrno2(8)

    Figure 1.

    Diagram from the original box to the conversion box in statistical mechanics.

    where rn(t) is the atom position at time t.

    The diffusion coefficient is one-sixth of the slope of the mean square displacement, in Figure 2.

    Figure 2.

    Description of the probability to find a particle within the shell r to r + Δr [44].

    As shown in Figure 2, the probability to find a pair of a distance r apart, relative to a uniform random distribution of particles at the same density, can be defined as g(r). For each particle, the probability to find another particle within the shell r to r + Δr is measured by g(r). Radial distribution function (RDF) is a measure of arrangement of atoms in a structure, and the MD is a deterministic method that offers the possibility of a microscopic description of a physical system in consideration of all the interactions involved. In simulation, an input file was generated for LAMMPS code to calculate the response of C▬S▬H structures with a large number of particles (atoms) to uniaxial loading, which leads to strain-stress data to calculate elastic modulus.

    2.1.2 Parameters of ClayFF field and programming code by LAMMPS

    The total potential energy of the classical molecular force field expands to:

    u=bondsus+anglesub+torsionsut+oopsuo+improoersui+invuinv+crossucross+vdWuvdW+clucl(9)

    The total energy is the sum of Coulombic (electrostatic) interactions, short-range interactions (van der Waals), and bonded (stretching/angular) interactions [27]:

    Utotal=Uvdw+Ucoul+Ubond+Uangle(10)

    For bonded intramolecular interactions, the Coulombic and VDW interactions are excluded. The Coulombic energy is represented as:

    Ucoul=e24πε0i<jqiqjrij(11)

    where qi and qj are partial charges, e is the charge of the electron, while ε0 is the dielectric permittivity of vacuum (8.85419 × 10−12 F/m). Note that the Coulombic interaction is long range and requires techniques such as the Ewald sums to be properly calculated. The VDW interactions are represented with the conventional 12-6 Lennard-Jones function that includes the short-range repulsion and the attractive dispersion energy:

    Uvdw=i<jDijRijrij122Rijrij6(12)

    where D and Rij are empirical parameters derived from the fitting of the ClayFF model to a number of observed structural property data for oxides, hydroxides, and oxyhydroxides.

    The interaction parameters between the unlike atoms are calculated by the arithmetic mean rule for the distance parameter Rij, the geometric mean rule for the energy parameter Dij:

    Rij=12Ri+RjandDij=DiDj(13)

    Bond stretching energy is considered between O and H of either a hydroxyl or a water molecule and is described by a simple harmonic term as:

    Ubond=k1rijr02(14)

    where k1 is a force constant and r0 represents the equilibrium bond length, both values taken from the flexible version of the SPC water model. For the description of the vibrational motion of hydroxyl groups, a bending (three-body) term is introduced in form of a harmonic relationship:

    Uangle=k2θijkθ02(15)

    where k2 is a force constant of Eq. (15), θijk is the bond angle for the oxygen-hydrogen, and θ0 refers to the equilibrium bond angle.

    2.1.3 Parameters of ClayFF field used in C▬S▬H structures

    In this study, bonded interactions include bond stretching potential energy function and bond angle bending potential energy function, which can be described by a harmonic potential function [47]. For nonbonded parameters of ClayFF force field, Lennard-Jones (12-6) potential is mainly considered. The empirical interatomic potentials in ClayFF field have been described above, and the corresponding parameters used for C▬S▬H structures are listed in Tables 1–3.

    (a) Bond stretching potential parameters
    Species iSpecies jK1 (kcal mol−1 Å−2)R0 (Å)
    owhw554.13491.0000
    ohho554.13491.0000
    ohsho554.13491.0000
    (b) Bond angle bending potential parameters
    Species iSpecies jSpecies kK1 (kcal mol−1 rad−2)θ0 (degree)
    hwowhw45.7696109.4700

    Table 1.

    Bonded parameters of ClayFF force field [27].

    TypeR0 (Å)D0 (kcal/mol)ElementCharge (ē)Atom type specification
    h*4.57750.0000H0.410Water hydrogen
    ho4.57750.0000H0.425Hydroxyl hydrogen
    o*3.55320.1554O−0.820Water oxygen
    oh3.55320.1554O−0.950Hydroxyl oxygen
    ob3.55320.1554O−1.050Bridging oxygen
    St3.70641.8405E−06Si2.100Tetrahedral silicon
    cao6.24845.0298E−06Ca1.360Octahedral calcium
    obss3.55320.1554O−1.300Bridging oxygen with double substitution
    obts3.55320,1554O−1.169Bridging oxygen with tetrahedral substitution
    obos35,5320.1554O−1.181Bridging oxygen with octahedral substitution
    ohs3.55320.1554O−1.081Hydroxyl oxygen with substitution
    cah6.24285.0298E−06Ca1.050Hydroxide calcium
    oh-3.55320.1554O−1.410Oxygen in hydroxyl anion
    Ca3.22370.1000Ca2.000Aqueous calcium ion

    Table 2.

    The nonbonded parameters of ClayFF force field [27].

    In fact, development of new parameters is very difficult, we mainly look for force field and select suitable parameters to be programmed in LAMMPS software, as many existed and complex force fields at present are difficult to be fully taken into account in this developing software.

    2.2 Modeling of typical C▬S▬H structures considering the Ca/Si ratio

    2.2.1 Description of C▬S▬H structures with various Ca/Si and Qn ratios

    In neat Portland cement, only the C▬S▬H with the highest Ca/Si ratio (>∼1.5) is observed, whereas C▬S▬H with the whole compositional range may exist in cement pastes containing fly ash, metakaolin, or silica fume [48]. 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 his scale, C▬S▬H gel exists in two forms: low-density C▬S▬H and high-density C▬S▬H. Two typical C▬S▬H structures of LD and HD C▬S▬H have been provided by Jennings [49] and Masoero et al. [50] in Figure 3.

    Figure 3.

    Structures at nanoscale and shear moduli of two typical C▬S▬H [49, 50]. (a) Two typical C▬S▬H structures at nanoscale [49] and (b) Shear moduli of typical LD and HD C▬S▬H structures [50].

    (a) Charge and species of atoms of C▬S▬H (C/S = 1.67) structure used [27]
    SpeciesCharge (ē)D0 (kcal/mol)R0 (Å)
    Water hydrogen (Hw)0.41
    Hydroxyl hydrogen (Ho)0.42
    Water oxygen (Ow)−0.820.15543.5532
    Hydroxyl oxygen (Oh)−0.950.15543.5532
    Bridging oxygen (O)−1.050.15543.5532
    Silicon (Si)2.11.84E−063.7064
    Calcium (Ca)1.055.03E−066.2428
    (b) The nonbonded Lennard-Jones (12-6) of C▬S▬H (C/S = 1.67) structure [27]
    Species iSpecies jDijRijSpecies iSpecies jDijRij
    CaCa5.0298E−066.2484SzOw0.0005348023.6298
    CaCw5.0298E−066.2456SzSz1.8405E−063.7064
    CaSz3.04259E−064.9774SzO0.0005348023.6298
    CaOss0.0008840994.9008SzOss0.0005348023.6298
    CaO0.0008840994.9008OssO0.15543.5532
    CaHw03.1242OssOss0.15543.5532
    CaOw0.0008840994.9008OssHw01.7766
    CwHw03.1214OssOw0.15543.5532
    CwOw0.0008840994.898OO0.15543.5532
    CwSz3.04259E−064.9746OOw0.15543.5532
    CwCw5.0298E−066.2428OHw01.7766
    CwO0.0008840994.898OwOw0.15543.5532
    CwOss0.0008840994.898OwHw01.7766
    SzHw01.8532HwHw00

    Table 3.

    Parameters for ClayFF field of C▬S▬H (C/S = 1.67) structure used [27].

    From Figure 3, from SEM image of C▬S▬H, two typical structures at nanoscale are provided, of which the LD C▬S▬H solid phase with the diameter of C▬S▬H gel is less than 16.6 nm, the gel porosity is about 35–37%, while the gel porosity of HD C▬S▬H solid phase is 24%. The solid porosity manifests itself (about 18%) at a scale smaller than the characteristic solid dimension of 2.2 nm [49, 50], irrespective of the type of C▬S▬H. It should be noted that, beyond this scale of nanoporosity, there is another type of porosity, which was found to differ from two C▬S▬H structures. This solid phase is named “globules,” which was found to have a characteristic size of 5.6 nm.

    Relation of C▬S▬H phase at microscale and “globule” C▬S▬H at nanoscale is linked, in Figure 4.

    Figure 4.

    Relation of C▬S▬H phase at microscale and “globule” C▬S▬H at nanoscale.

    From Figure 4, the porosity is an important factor to reveal the differences between “globule” C▬S▬H at nanoscale and LD/HD C▬S▬H phase at microscale. Then elastic modulus of the “globule” C▬S▬H about 5.5 nm3 can be adopted to assess the elastic moduli of LD and HD C▬S▬H phases, where the values are verified by the outer and inner C▬S▬H distinguished in nanoindentation experiment.

    It is generally considered that C▬S▬H solid solution structure is a kind of tobermorite structure between CH structures to form a sandwich structure [51]. Calcium-rich model and silicon-rich structural models: depending on the Ca/Si ratio, C▬S▬H can be divided into silicon-rich C▬S▬H (Ca/Si = 0.65–1.0) and calcium-rich C▬S▬H (Ca/Si = 1.1–1.7) [52]. Structure and defect model of C▬S▬H [53] and the NMR measurement of Grutzeck et al. [19, 52], Q1/Q2 ∼ C/(S + A) are as in Figure 5.

    Figure 5.

    Structure and defect model of C▬S▬H and the NMR measurement data of Q1/Q2 ratio. (a) Structure and defect model of C▬S▬H [53] and (b) Q1/Q2 ratio of various C▬S▬H [19, 52].

    From Figure 5(a), we can see that such nanostructures and intermediary structural models consider that C▬S▬H composition is uniform within the scope of nanocrystals (about 5 nm). However, in the range of short-range order of 1 nm, C▬S▬H composition and structure can vary, especially in the area of matrix composition of the amorphous C▬S▬H. From Figure 5(b), the results show that the measured Ca/Si ratio has a mutation in both sides of 1.0–1.1 region, which proved that C▬S▬H between 1.0 and 1.1 is vacant. However, Q1/Q2 in calcium-rich C▬S▬H is fluctuated from 1.0 to 1.5, which can be considered about 50% dimers and 50% repeated silicon-oxygen tetrahedron chain ternary.

    There, a Qn nomenclature is used in general for the peaks. Qn is the chemical shift of a silicon atom, which is bound to n bridging oxygens. The Qn distribution and the mean chain lengths (MCL) of C▬S▬H structures with various Ca/Si ratios are given by Tajuelo et al. [54], shown in Table 4.

    Ca/SiQ1 (%)Q2 (%)Q3 (%)MCL
    0.759.579.810.718.8
    0.8319.177.13.810.1
    135.762.02.35.5
    1.2574.026.002.7
    1.3379.920.102.5
    1.587.712.302.3

    Table 4.

    The Qn distribution and mean chain lengths (MCL) of C▬S▬H with various Ca/Si ratios [54].

    From Table 4, when the C▬S▬H structure with a C/S ratio higher than 1.0, the number of Q3 is zero, which means there is no silicate tetrahedral with three bridging oxygen atoms in its structure.

    2.2.2 Structural, bonding, and interaction of C▬S▬H structures

    For modeling of a 14Å tobermorite: (C/S = 0.83), we chose Bonaccorsi model [55] with the empirical formula of Ca5Si6O16(OH)2·7H2O (Ca/Si = 0.83), and crystallographic parameters are as follows: (1) crystal system: triclinic; (2) the cell parameters: a = 6.735, b = 7.425, c = 27.987 Å, α = 90°, β = 90°, γ = 123.25°; (3) space group: B11b. Figure 6 shows the structure and atomic connections of 14 Å tobermorite.

    Figure 6.

    The structure and atomic connections of 14 Å tobermorite structure [55]. (a) The real cell of 14Å tobermorite by Bonaccorsi [55]. (b) Silicon chain links and atomic connections.

    From Figure 6, 14 Å tobermorite structure is basically a layered structure. Besides, the central part is a Ca▬O sheet (with an empirical formula: CaO2, although the chemical aspect CaO should always be CaO, which implies that the oxygen in CaO2 also includes that of the silicate tetrahedron part). Moreover, silicate chains envelope the Ca▬O sheet on both sides, with the characteristics mentioned later. Last but not least, between individual layers, Ca2+ and H2O are filled into the space to balance the charges and determine the layer distance, respectively.

    For modeling of the jennite structure [56] (C/S = 1.5), the ordered jennite is triclinic, space group P-1, with the unit cell dimensions of a = 10.576 Å, b = 7.265 Å, c = 10.931 Å, α = 101.3°, β = 96.98°, and γ = 109.65° [57]. The chemical constitutional formula is Ca9Si6O18(OH)6·8H2O and Ca/Si ratio is 1.5. Figure 7 shows the structure and atomic connections of jennite unit cell.

    Figure 7.

    Modeling structure and atomic connections of jennite crystal structure. (a) The topology structure of jennite crystal [56] and (b) the real structure of jennite crystal [57] (atom types: green, Ca; yellow, Si; red, oxygen; white, hydrogen).

    We can see in Figure 7 that jennite has the same link of the silicate tetrahedral chains similar to 14 Å tobermorite [55]. The difference is that only half of the oxygen atoms of the Ca▬O layer in jennite structure are connected to the silicate tetrahedral chains, while the other half of the oxygen atoms are linked to the hydroxyl groups. Moreover, the bridging tetrahedra of the silicate chains are connected to the Ca▬O layer.

    For modeling of a 11 Å tobermorite: (C/S = 0.67), the morphology of 11 Å tobermorite is monoclinic and the initial lattice is as follows: a = 6.69 Å, b = 7.39 Å, c = 22.779 Å, α = β = 90°, γ = 123.49°, its structure is monoclinic with space group P21 [58]. Modeling of 11 Å tobermorite is shown in Figure 8.

    Figure 8.

    Modeling of 11 Å tobermorite crystal. (a) 11 Å tobermorite by Hamid [58]; (b) silicon chain links and atomic connections; and (c) the real unit cell.

    In Figure 8(a), the 11 Å tobermorite crystal can be summarized as follows: (1) the structure is basically a layered structure with the central part as a Ca▬O sheet (with an empirical formula of CaO2, of which the oxygen in CaO2 also includes that of the silicate tetrahedron part). (2) Silicate chains envelop the Ca▬O sheet on both sides in Figure 8(b). (3) Ca2+ and H2O are filled between individual layers to balance the charges. The unit cell is in Figure 8(c), the infinite layers of calcium polyhedra, which is parallel to (001), with tetrahedral chains of wollastonite-type along b and the composite layers stacked along c and connected through formation of double tetrahedral chains [59].

    2.2.3 Radial distribution function (RDF) of C▬S▬H structures with different Ca/Si ratios

    It should be pointed out that the description for C▬S▬H structures by 11 Å tobermorite structure is not fixed, and the position and the appearance of certain atoms can be changed and adjusted, thus leading to various Ca/Si ratios. There are three different structures with various Ca/Si forms based on Hamid model [60], of which 11 Å tobermorite is denoted as Ca8Si12O28(OH)8·4H2O (Ca/Si = 0.67), Ca10Si12O32(OH)4·4H2O (Ca/Si = 0.83), and Ca12Si12O36·4H2O (Ca/Si = 1.00). The jennite (Ca9Si6O18(OH)6·8H2O, Ca/Si = 1.50) is also given [56]. Then a total radial distribution function (RDF) and atomic partial RDF curves of C▬S▬H structures with different C/S ratios (0.67–1.5) can be calculated, shown in Figure 9.

    Figure 9.

    The total radial distribution function of C▬S▬H structures with different Ca/Si ratios. (a) A three-dimensional view of tobermorite structure [58] and (b) the total RDF of C▬S▬H with different Ca/Si ratios.

    It can be seen in Figure 9(a) that, for the system of Ca/Si = 0.66 (or Water/Si = 0.67), four sites of O8, O9, O14, O18 are protonated caused by the lack of Ca5 and Ca6 atoms; For the system of Ca/Si = 0.83 (or Water/Si = 0.50), two sites of O8, O9, O14, O18 are protonated caused by the lack of Ca6 atoms. When Ca/Si = 1.0 (or water/Si = 0.33), there are only two water molecules of W1 and W2, all sites of O8, O9, O14, O18 are not protonated because both Ca5 and Ca6 atoms are present. In Figure 9(b), there are three kinds of atomic partial RDF curves, of which the chemical bonds include Ca▬O, Si▬O, Si▬Si. The first peak is OH, the second peak is Si▬O, the third peak Ca▬O, the fourth peak is Si▬Si, the respective location is 1.28, 1.56, 2.54, and 3.10 Å, for the interatomic distance between corresponding pairs of atoms.

    3. Construction of “glouble” C▬S▬H by 11 Å tobermorite with various C/S ratios

    As the initial configuration is monoclinic, in order to study the deformation of independence in a particular direction, the orthogonal box conversion of the C▬S▬H structure is needed from the initial monoclinic one, shown in Figure 10.

    Figure 10.

    Schematic conversion diagram of simulation box in x-y plane. (a) Original box of a monoclinic cell and (b) transformed box of a orthogonal cell.

    From Figure 1(a)), it is noteworthy that the relative atomic position and model size are unchanged in C▬S▬H model after box transformation in Figure 1(b) with the conversion modes of orthogonal box by LAMMPS. Here the conversion is done to obtain the orthogonal structure as an analog initial structure.

    3.1 Construction of crystalline 11 Å tobermorite structures (C/S = 0.67)

    In order to obtain an amorphous structure in the long-range disorder and short-range order, silicon-oxygen tetrahedra chain is rotated and twisted. Above all, the system was a preequilibrated ensemble of NVT system at 3000 K, and the temperature of the system was gradually lowered to 300 K. After 300 K rebalance, the output is written as an input of the NVE heald, then we obtain stable amorphous C▬S▬H. Finally, the simulation results of the C▬S▬H at 300 K in terms of structural properties and dynamical behavior are analyzed.

    Construction of a monolithic C▬S▬H cell structure corresponding to “globule” C▬S▬H is in Figure 11.

    Figure 11.

    Nanoscale modeling of C▬S▬H with a size of a “globule” by Jennings. (a) 11 Å tobermorite, (b) model of 2.5 nm3, (c) amorphous model of 5.5 nm3 (atom types: red-Ow + white-Hw → water molecule (H2O); red, O; green, Ca; purple, Ca; yellow, Si).

    As in Figure 11, based on single crystal of 11 Å tobermorite (C/S = 0.67) in Figure 3(a) discussed in Chapter 3, two supercells 4a × 3b × 1c in Figure 3(b) and 4a × 3b × 2c of 11 Å tobermorite model are established first. Then the amorphous model of 5.5 nm3 is obtained after MD simulation considering an annealing treatment, shown in Figure 3(c). MD simulation of monolithic C▬S▬H is seen as the “glouble” C▬S▬H.

    The 11 Å tobermorite structure is transformed from crystalline to amorphous by means of an annealing process simulated by molecular dynamics simulation. The modeling procedure is deeply presented in a recent publication of the authors [2], that is why only a brief description is developed here.

    An assembly of a supercell with the size of 8a × 6b × 2c (5.352 × 4.434 × 4.554 nm3) of single 11 Å tobermorite cells. These dimensions are chosen to be coherent with the model of Jennings, which describes the C▬S▬H gel as a packing, more or less dense, of monolithic C▬S▬H (with a nanogranular aspect of about 5 nm diameter) [49].

    Molecular dynamics is a computational simulation method, which aims to study the physical movements of atoms and molecules, by resolving Newton’s equations. LAMMPS software is used to perform this study. ClayFF potential is retained to simulate the atomic force field between the various atoms of the supercell.

    The various steps to simulate the annealing process and to obtain an amorphous structure are:

    • energy minimization (conjugate gradient method) and balance of the structure at 300 K separately under NVT and NVE systems;

    • increase of the temperature until 3000 K under NVT system, then balance of the structure at this high temperature until NVT and NVE systems;

    • decrease of the temperature from 3000 to 300 K under NVT system;

    • balance of the structure at 300 K under NVT and NVE systems.

    After this annealing treatment, the ordered silicon chain becomes disordered and the system is more stable in comparison with the initial structure.

    3.2 Construction of amorphous 11 Å tobermorite structure (C/S = 1.67)

    Modeling of an amorphous C▬S▬H structure (C/S = 1.67) with various sizes is in Figure 12.

    Figure 12.

    Modeling of an amorphous C▬S▬H structure (C/S = 1.67) with various sizes. (a) A 2.676 × 2.217 × 2.278 nm3 model, (b) a 5.352 × 4.434 × 4.556 nm3 model and (c) amorphous 5.5 nm3 model (atom types: light blue-Ow + white-Hw → water molecule (H2O); red, O; green, Ca; purple, Ca; yellow, Si).

    In Figure 12(a), a 3D image of a C▬S▬H structure with a C/S ratio equal to 1.67 and with a distribution of Q0 = 13%, Q1 = 67%, and Q2 = 20% is presented. Amorphous C▬S▬H (C/S = 1.67) models in Figure 12(c) are constructed by periodic extension in x, y and z directions of the primary cell of 5.352 × 4.434 × 4.556 nm3 with 8328 atoms presented in Figure 12(b), with an initial density of 2.257 g/cm3.

    Based on the layered supercell of 11 Å tobermorite (Ca6[Si6O18]·2H2O, Ca/Si = 1.0), an amorphous C▬S▬H(II) (C/S = 1.67) supercell (4a × 3b × 1c) is constructed, as shown in Figure 13.

    Figure 13.

    Construction of an amorphous C▬S▬H (C/S = 1.67) supercell of 4a × 3b × 1c 11 Å Tob. (a) Initial 2.676 × 2.217 × 2.278 nm3 model; (b) removal of light silicate tetrahedral along [010]; (c) along [001] (atom types: red-Ow + white-Hw → water molecule (H2O); red, O; green, Ca; purple, Ca; yellow, Si). (d) GCMC water adsorption (initial); (e) GCMC water adsorption (half); (f) final amorphous C▬S▬H (atom types: light blue-Ow + white-Hw → water molecule (H2O); red, O; green, Ca; yellow, Si).

    Figure 13 records the C▬S▬H model construction procedures, which are as follows: above all, from Figure 13(a), a supercell of 4a × 3b × 1c 11 Å tobermorite is periodically established by Materials Studio software. Then a cleave surface is built and makes it possible to change the box into an orthogonal type [61], which can be used to better understand mechanical bahaviors of cement paste using MD method by Murray et al. [62]. The simplified formula is (CaO)1.67(SiO2)(H2O)1.75, which is close to the reference [28]. As presented in Figure 13(b) and Figure 13(c), the light silicate tetrahedrals are broken along [010] and [001] and the layered supercell of 11 Å tobermorite without water is taken as the initial configuration. Secondly, Figure 13(d) describes a dry C▬S▬H in which the continuous silicate chains were broken to achieve Q0, Q1, and Q2 percentages of 13, 67, and 20%, respectively, without considering the O▬H bonds.

    Finally, Grand Canonical Monte Carlo (GCMC) simulation of the water adsorption [63] by LAMMPS using the ClayFF parameters (in Tables 1–3) is carried out to obtain the disordered C▬S▬H structure in Figure 13(e) and (f)), which is performed on the initially dry C▬S▬H structure at 300 K. The number of atoms is about 4164, with the initial density of 2.56 g/cm3. The hydrated phase is presented at various stages of the simulation in Figure 13(d)(f)). At the equilibrium state, the C▬S▬H reached saturation shown in Figure 13(f)), which is also the same with the disordered C▬S▬H structure shown in Figure 13(c)). The final chemical composition and density of the amorphous disordered C▬S▬H are (CaO)1.67(SiO2)(H2O)1.75 and 2.257 g/cm3, which are in a good agreement with the results obtained from the neutral scattering tests [22]. As for a C/S ratio of 0.67, this supercell represents the monolithic full-dense C▬S▬H, in other words, the “globules” by Jennings [49].

    4. Simulation and properties of mechanical two kinds of C▬S▬H (I) (C/S = 0.67)

    In order to get the stress-strain relation of C▬S▬H, the model was subjected to uniaxial tensile loading through gradual elongation with a strain rate of 0.001 ps−1. In the whole simulation process, NPT ensembles are set for the system. Nosé Hoover thermostat temperature control method and the Nosé Hoover barostat pressure method were used for NPT calculations. A time step of 1 femtosecond (fs) was chosen since 1 fs is suitable for most purposes. The dynamic time ranges from 100 to approximately 400 picoseconds (ps) depending on the size of the simulation cell, which ensures thermodynamic equilibrium and the convergence of energy with a reasonable computational time. For example, the procedures in the uniaxial tension test in the x direction are as follows. Firstly, the supercells are relaxed at 300 K and coupled to zero external pressure in the x, y, and z dimensions for 50 ps. Then, after the pressures in three directions reach equilibrium, the C▬S▬H structure is elongated in the x direction. Meanwhile, the pressure in y and z direction should be kept at zero. Pressure evolution in the x direction is taken as the internal stress σxx. The pressure perpendicular to the tensile direction is set to zero, which allows the normal direction to relax anisotropically without any restriction. The setting, considering the Poisson’s ratio, can eliminate the artificial constraint for the deformation. Additionally, in the simulation process, 10,000 configurations are recorded for structural analysis.

    4.1 Tensional simulation of the crystaline and amorphous C▬S▬H(I)

    4.1.1 Tensional simulation of the crystalline C▬S▬H(I)

    The crystalline 11 Å tobermorite structure is obtained first and then the tensional deformation under a strain rate of 10−3 ps−1 is simulated. We assume that, for a tensile strain rate of 10−3 ps−1, quasistatic conditions are obtained. The crystaline C▬S▬H(I) structure (4a × 3b × 1c) by 11 Å tobermorite is obtained first, and then the tensional deformation under the strain rate 1.0 × 10−3 ps−1 is simulated. The tensile deformation process of C▬S▬H(I) structure with the size of 2.676 × 2.217 × 2.278 nm3 are calculated and analyzed. Take the z direction, for example, the structure evolution of the difference strain under the strain rate of 1 × 10−3 ps−1 is in Figure 14.

    Figure 14.

    Results of a crystalline C▬S▬H(I) structure for various strains under 10−3 ps−1 strain rate.

    From Figure 14, it is indicated that the water molecule has great impact on the mechanical properties of C▬S▬H(I) during tensile deformation. Moreover, for crystalline C▬S▬H(I) by 11 Å tobermorite structure (4a × 3b × 1c) with the size of 2.676 × 2.217 × 2.278 nm3, it is likely to be broken in the z direction, while it is difficult to break in the x and y directions.

    During the stretching process, the aqueous layer is easily broken along the z direction, and the stratified phenomenon becomes apparent at both sides. A comparative study of the behaviors of crystalline and amorphous C▬S▬H(I) is presented in this section. In order to maintain comparative structure, use the same initial structure with the size of 2.676 × 2.217 × 2.278 nm3.

    4.1.2 Tensional simulation of an amorphous C▬S▬H(I) (C/S = 0.67)

    After annealing treatment on Lammps, the ordered silicon chain becomes disordered and the system becomes more stable compared with the initial structure, which is consistent with the real structure of C▬S▬H(I). Then the stretching process of an amorphous C▬S▬H (the ratio of Ca/Si is 0.67) at a certain strain rate 1 × 10−3 ps−1 is simulated. Take the z direction, for example, the nanostructure evolution of the difference strain under the 1 × 10−3 ps−1 is shown in Figure 15.

    Figure 15.

    Aggregation of water molecules during tensile deformation under different strains. Deformation during a tensional test of a 4a × 3b × 1c supercell model after annealing processes (atom types: light blue-Ow + white-Hw → water molecule (H2O); red, O; green, Ca; yellow, Si).

    From Figure 15, it is indicated that the water molecule has a great impact on the mechanical properties of C▬S▬H(I) during tensile deformation. With an increasing strain, the water molecules in an aqueous layer have a tendency to be pulled off, which may explain the softening behaviors of C▬S▬H(I) during tensile process. In the structure after annealing, the distribution of water molecules changes a lot, and water molecules become dispersed disorderly. During stretching process, the aqueous layer is only formed in the upper layer, and the deformation resistance is much greater than that in z direction in the crystalline structure.

    4.2 Elastic moduli results of the crystalline and amorphous C▬S▬H(I)

    4.2.1 Elastic moduli for a crystalline C▬S▬H(I) (C/S = 0.67) in various directions

    The tensional process of an amorphous C▬S▬H(I) is calculated and analyzed. The corresponding tensile and compressive stress-strain curves are in Figure 16.

    Figure 16.

    Results of a crystalline C▬S▬H(I) structure under strain rate of 10−3.

    As is shown in Figure 16, for the crystaline C▬S▬H(I) structure (the supercell 4a × 3b × 1c structure by 11 Å tobermorite) using the ClayFF parameters, Young’s modulus in three directions is averaged to about 79.95GPa, which is close to the value of 79.51 GPa obtained by DFT seen in Chapter 3. Similarly, by using BHM potential, the influencing rule of loading direction on strain-stress curves during the elastic linear strain is close to the results by using ClayFF parameters. The elastic moduli in the corresponding x, y, and z directions are 65.61, 82.97, and 85.92 GPa, respectively. Moreover, parameters of ClayFF field seem more stable than those of BHM field (not shown here).

    4.2.2 Elastic moduli for an amorphous C▬S▬H(I) (C/S = 0.67) in various directions

    With the same force field, the tensile test is then simulated in the three directions. The stretch and compression process of an amorphous C▬S▬H(I) structure (4a × 3b × 1c supercell of 11 Å tobermorite, 5.352 × 4.434 × 4.556 nm3) are calculated and analyzed. The corresponding tensile stress-strain curves are in Figure 17.

    Figure 17.

    Results of an amorphous C▬S▬H(I) structure after annealing by ClayFF. (a) Supercell structure of 5.352 × 4.434 × 4.556 nm3 and (b) Comparison of C▬S▬H(I) with different sizes.

    From Figure 17, the three curves are close to each other, showing the isotropy of the supercell as well as its amorphous structure. From these results, by determining the slope at the beginning of these curves (in between 10 and 50% of the maximal stress), the homogeneous Young’s modulus may be determined: it is found to be equal to 59.4 GPa. The annealing process leads to a decrease of this elastic property from 79.5 to 59.4 GPa; this result is in agreement with previous results of Vandamme et al. [65]. This value will be used as homogenized Young’s modulus to calculate the specific heat of monolithic C▬S▬H.

    In a word, in the disordered structure obtained by annealing treatment, the anisotropy in three directions significantly decreased. It can be inferred that, if the disordered structure becomes much larger, it leads to the isotropic behaviors in a certain degree.

    4.2.3 Comparison of crystalline and amorphous C▬S▬H(I) during tensional process

    We can see from Figures 16 and Figures 17 that, under the strain rate of 10−3 ps−1, peak stress does not change obviously before and after annealing. Result comparison of crystalline and amorphous 11 Å tobermorite supercells during tensional process is shown in Table 5.

    Crystal/structureInitial size ( nm3)Cell/box typeMethodExEyEzE = (Ex + Ey + Ez)/3
    11Å tobermorite0.669 × 0.739 × 2.278Single crystalDFT79.51279.512
    11Å tobermorite [60]0.669 × 0.739 × 2.278Single crystalDFT78.93978.939
    C▬S▬H(I) (C/S = 0.67)2.676 × 2.217 × 2.278CrystallineClayFF (MD)55.8286.6397.4079.95
    C▬S▬H(I) (C/S = 0.67)2.676 × 2.217 × 2.278DisorderedClayFF (MD)67.8984.6378.6277.05
    C▬S▬H(I) (C/S = 0.67)5.352 × 4.434 × 4.556DisorderedClayFF (MD)57.8566.1354.1459.37

    Table 5.

    Elastic modulus comparison of C▬S▬H(I) structures after optimization.

    From Table 5, we can see that after annealing process, Young’s modulus of C▬S▬H(I) (C/S = 0.67) structure with a size of 2.676 × 2.217 × 2.278 nm3 using BHM potential slightly decreased from 78.17 to 68.49 GPa. However, with the ClayFF potential, a decrease of the Young’s modulus after annealing process under the strain rate of 10−3 has not been pointed out.

    Similarly, by using the ClayFF potential, the influencing role of loading direction on strain-stress curves during the elastic linear strain is very close to the results by using BHM potential. After annealing treatment on Lammps, the elastic modulus using the ClayFF field in the corresponding x, y, and z directions are averaged to be 67.89, 84.63, and 70.53 GPa, respectively. Moreover, for C▬S▬H(I) structure with various sizes at nanoscale in Table 5, it seems that the size has a significant influence on the Young’s modulus of the amorphous material.

    In all, as a kind of C▬S▬H(I), 11 Å tobermorite supercell (C/S = 0.67) becomes a denser and more stable configuration after annealing treatment, which is in agreement with the previously proposed quantitative “colloid” model of C▬S▬H gel, resulting in an improved understanding of the microstructural changes associated with drying and heat curing [64].

    5. Simulation and mechanical properties of an amorphous C▬S▬H(II) (C/S = 1.67)

    Before simulation, it must be also noted that the total energy obtained by this method (called lattice energy for crystalline systems) has no physical sense for most force fields. After repeated tests, ClayFF force field is used to complete the MD simulation.

    5.1 Tensional simulation of amorphous C▬S▬H(II) structures

    The stretching of the C▬S▬H structure of 5.352 × 4.434 × 4.556 nm3 at the strain rate of 10−3 ps−1 under various directions is shown in Figure 18.

    Figure 18.

    The stretching C▬S▬H structure at the strain rate of 10−3 ps−1 under various directions.

    From Figure 18, the aggregation of water molecules changes in different directions during the stretching process. The layers of water molecules in the x direction are extended and the width of the water molecules layer becomes larger, while the water molecule layer in the y direction is narrowed. The water molecules gathered toward the center of the structure in the z direction; the different states of aggregation form of water molecules in various directions cause the anisotropic properties of C▬S▬H structure with size of 5.352 × 4.434 × 4.556 nm3.

    5.2 Elastic moduli of the amorphous C▬S▬H(II) determination

    In order to better understand the influence of size on the performance of C▬S▬H(II) structures (C/S = 1.67), the stretching process is simulated to obtain stress-strain curves and then to find the influencing rules of these parameters on the properties of C▬S▬H(II) structures (C/S = 1.67).

    5.2.1 Effect of loading direction on strain-stress curves of size 5.352 × 4.434 × 4.556 nm3

    Figure 19 shows the tensional stress-strain curves in x, y, and z directions of an amorphous C▬S▬H structure of size 5.352 × 4.434 × 4.556 nm3 under the strain rate of 10−3 ps−1.

    Figure 19.

    Comparison of supercell structures under various directions by MD simulation. (a) Supercell structure of 5.352 × 4.434 × 4.556 nm3 and (b) comparison of two supercell structures.

    Figure 19(a) shows the strain-stress curves of a supercell structure with size of 5.352 × 4.434 × 4.556 nm3 under various directions. From Figure 19(b), we can see that the maximum stress of tensional stress-strain curve with the size of 2.676 × 2.217 × 2.278 nm3 is higher than the other size of 5.352 × 4.434 × 4.556 nm3. With the increase of size, the peak stress becomes lower. It can also be seen that the presence of the maximum stress in stress-strain curves greatly varies due to the scale effect between atomic interactions. Young’s modulus with the size of 5.352 × 4.434 × 4.556 nm3 is lower than the other size of 2.676 × 2.217 × 2.278 nm3, which may be due to the size effect at nanoscale. However, there is little change of stress-strain curve slope in the initial elastic stage. Moreover, the averaged modulus in three directions is 60.95 GPa under strain rate of 10−3 at 300 K, which is close to the other results of the C▬S▬H structure with the value of 63.50 [25] and 69.01 GPa [55]. A supercell of 5.352 × 4.434 × 4.556 nm3 as a monolithic full-dense C▬S▬H (globules) appears to be sufficient to determine Young’s modulus, and it can be used to further estimate Young’s moduli of LD and HD C▬S▬H, which is the aim of this chapter.

    Elastic modulus comparison of C▬S▬H constructed by 11 Å tobermorite is in Table 6.

    Crystal/structureInitial size (nm3)Box typeMethodExEyEz(Ex + Ey + Ez)/3
    C▬S▬H(I) (C/S = 0.67)2.676 × 2.217 × 2.278CrystallineMD55.8286.6397.4079.95
    C▬S▬H(I) (C/S = 0.67)5.352 × 4.434 × 4.556DisorderedMD57.8566.1354.1459.37
    C▬S▬H(II) (C/S = 1.67)5.352 × 4.434 × 4.556DisorderedMD55.8258.8968.1360.95

    Table 6.

    Elastic modulus comparison of C▬S▬H structures constructed by 11 Å tobermorite.

    From Table 6, we can see that the deformation direction has a little influence on the mechanical properties of the amorphous C▬S▬H(II) (C/S = 1.67) under the strain rate of 10−3 ps−1. Meanwhile, for the structure of C▬S▬H(II) (C/S = 1.67) with size of 5.352 × 4.434 × 4.556 nm3, the corresponding averaged elastic modulus in y and z directions is higher than that in x direction by the stress-strain curve slope within about 0.002–0.055 strain. Moreover, there is only a slight difference between elastic moduli of C▬S▬H(I) and C▬S▬H(II). For all results of the C▬S▬H(I) structures in three directions mentioned in Table 5, the averaged value is very close to each other by using ClayFF parameters. Conversely, the size has a significant influence on the Young’s modulus of the amorphous C▬S▬H(II) structure.

    5.2.2 Assessment of the elastic properties of amorphous C▬S▬H structures by MD

    As the solid phase described by Jennings is named by “globules,” which was found to have a characteristic size of about 5.6 nm [49]. Here we can use the elastic modulus of the CSH(II) structure with a size of 5.352 × 4.434 × 4.556 nm3 to approximate the modulus of “globules” so as to investigate the LD C▬S▬H and HD C▬S▬H.

    A self-consistent scheme may be used [65]:

    Kks=4G1ϕs4G+3ksϕs(16)
    Ggs=1254ϕs316rs3ϕs+11664320ϕs+400ϕs2+rs144168ϕs120ϕs2+rs28154ϕs+9ϕs2(17)

    where rs=ksgs=21+μs312μs, gs is shear modulus of the monolithic full-dense C▬S▬H (here gs = Es/[2(1 + μs)] with Es = 60.95GPa and μs = 0.25 for C▬S▬H(II) found in this work). ks is bulk modulus of the monolithic full-dense C▬S▬H (here ks = Es/[3(1 − 2μs)] with Es = 60.95GPa and μs = 0.25 for C▬S▬H(II) structure).

    The Mori-Tanaka scheme may also be used [65]:

    K=rfrkr1+αskrks11rfr1+αskrks11(18)
    G=rfrkr1+βsgrgs11rfr1+βsgrgs11(19)
    αs=3ks3ks+4gs(20)
    βs=6ks+2gs53ks+4gs(21)

    where αs and βs are coefficients, fr is the volume fraction, k is the bulk modulus, g is the shear modulus. The subscripts “o” stands for the monolithic C▬S▬H.

    Figure 20 presents the evolution of the Young’s modulus according to the porosity.

    Figure 20.

    Young’s modulus evolution of C▬S▬H(II) (C/S = 1.67) according to the porosity.

    As is shown in Figure 20, Young’s modulus of LD C▬S▬H and HD C▬S▬H can be estimated by a self-consistent scheme and the Mori-Tanaka scheme. For example, with the self-consistent scheme in Eqs. (16) and (17), and assuming a porosity of, respectively, 0.35 and 0.24, the Young’s moduli of LD C▬S▬H and HD C▬S▬H are found to be equal to 18.11 and 31.45 GPa, which is in agreement with the results of Constantinides (2006). Similarly, by using the Mori-Tanaka scheme in Eqs. (18)–(21), the values of Young’s moduli of LD C▬S▬H and HD C▬S▬H are, respectively, 29.78 and 37.71 GPa.

    6. Conclusions

    In this chapter, the stretch process of a represented amorphous full-dense C▬S▬H(I) and C▬S▬H(II) with the Ca/Si ratios of 0.67 and 1.67 is simulated by LAMMPS program after setting parameters of related potential functions. Some conclusions can be drawn.

    1. For two typical crystalline and amorphous C▬S▬H(I) structures (formula: Ca4Si6O14(OH)4·2H2O), the stretch process at a certain strain rate is calculated, and the corresponding tensile stress-strain curves are analyzed. The results are as follows:

      • Before annealing process, the crystalline supercell (4a × 3b × 1c) of C▬S▬H(I) at nanoscale has the averaged elastic modulus value of 79.95 GPa by MD simulation, which is similar to the homogenized Young’s modulus value of 79.512 GPa by elastic constants of single 11Å tobermorite model based on DFT.

      • Through annealing process of the C▬S▬H(I) with a size of 2.676 × 2.217 × 2.278 nm3, using ClayFF field parameters, the difference of Young’s modulus in three directions becomes smaller during the tensional deformation, with an averaged value of 77.05 GPa. For the C▬S▬H(I) structure with a size of about 5 nm, using MD simulation, the averaged Young’s moduli are separately about 59.37 GPa.

      • By comparison of both the crystalline and amorphous C▬S▬H(I) structures, it is found that the water molecule in both structures has a great influence on tensile strength.

      • After annealing simulation, the annealed structure becomes more uniform in three directions. Young’s modulus is very slightly reduced; however, we think that with a larger supercell corresponding to globule C▬S▬H, such a decrease would be obtained, indicating that the annealed supercell with disordered silicon chain is more softened and more consistent in x, y and z directions.

    2. A typical C▬S▬H(II) structure (formula: (CaO)1.67(SiO2)(H2O)1.75) is simulated by using the ClayFF field. Above all, a layered supercell of 11Å tobermorite of Hamid model without water is firstly considered as the initial configuration. Then, SiO2 groups in silica tetrahedra have been removed, guided by the NMR results to achieve realistic percentages of Q0, Q1, and Q2. Finally, Grand Canonical Monte Carlo simulation of water absorption on the dry C▬S▬H structure at 300 K is performed to obtain a disordered C▬S▬H structure.

      • The water molecule in the structure has a great influence on tensile strength.

      • It is difficult to find a RVE to assess the tensile strength and the postpeak part of the tensile curve. However, Young’s modulus of C▬S▬H(II) structure with a size of about 5 nm has been identified to be between 55 and 68 GPa (the averaged modulus in three directions is 60.95 GPa under strain rate of 10−3 ps−1 at 300 K), which is consistent with previous literature results. The Young’s modulus in x, y, and z directions is rather the same, which suggests that the model is indeed amorphous.

      • The effect of size on the response of the amorphous C▬S▬H structure has been also investigated; by the variation of stress-strain curve due to the atomic interactions, the presence of scale effect is observed (this means that, with the increase of size, the peak stress becomes lower). Moreover, Young’s modulus with the size of 5.352 × 4.434 × 4.556 nm3 is lower than the other size of 2.676 × 2.217 × 2.278 nm3.

      • Water layer is an important factor, which affects the mechanical properties, showing the weakest loading resistance perpendicular to the aqueous layer.

      Overall, we mainly considered the C/S influencing factor, which has an effect on the stiffness and strength of C▬S▬H gel. The presence of impurities in C▬S▬H gel (i.e., sulfate or aluminate), porosity, or a specified density should be also considered in future work.

    3. The values we have found for the Young’s moduli of amorphous C▬S▬H(I) and C▬S▬H(II) are those of the monolithic C▬S▬H (globules). At the upper scale, we can deduce the Young’s moduli of LD C▬S▬H and HD C▬S▬H assuming a porosity for these two kinds of C▬S▬H. This study is limited to the C▬S▬H(II), the kind of C▬S▬H found in the hardened cement pastes experimentally tested and simulation results in Chapter 6. We can see that there is a significant difference between the two homogenization methods, as previously pointed out by Constantinides.

    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

    Nanoindentation Investigation of Typical C▬S▬H Structures: FEM Simulation

    Abstract

    This work focuses on elastic modulus of two main constituents of cement-based materials: portlandite (CH) and calcium silicate hydrate (C▬S▬H). The large-scale models using QC methods are considered trustworthy and on the ABAQUS software, where, on behalf of the purely elastic scene, the unloading curves of the three (CH, LD and HD CSH, CSH) models are simulated. Based on the elastic modulus of the monolithic C▬S▬H structure by MD simulations, the assessment results on LD C▬S▬H and HD C▬S▬H after homogenization are very close to nanoindentation simulation data. The findings are as follows: (1) the homogenized elastic properties of polycrystals can be obtained by elastic constants of single crystal (using DFT and RVH estimation) and thus can be used to explain the relationship between structure and mechanical properties of CH from nano- to microscale. (2) It is found that the results of comparison of the simulated and experimental unloading curves of CH and LD/CH C▬S▬H are essentially coincident and have a very good agreement, showing the feasibility and the rationality of the simulation methods above. Results enable to provide useful parameters for composite cement system modeling and a method to calculate elastic modulus of other similar structures.

    Keywords: calcium silicate hydrate, elastic modulus, nanoindentation, FEM

    1. Introduction

    Multiscale simulations of atomistic/continuum coupling in computational materials science have become a hot spot for recent years, and the different scales will be considered to reflect the atomic interaction [1]. Multiscale modeling of structural concrete performance is presented as a systematic knowledge base of coupled cementitious composites and structural mechanics [2]. At the scale of micromechanical investigations, the interaction of the various phases is taken into account. In this way, for the greater length scale, a unique constitutive behavior is extracted that typically cannot be captured fully by standard closed-form continuum models [3]. An illustrative way for denoting these different length scales is the micro-, meso-, and macroscale [4]. Among these, the level of components or structural parts is called macroscale. At present, there are two major directions in multiscale modeling [5]: extending the atomistic length scales (spatial multiscale methods) and extending the atomistic time scales (temporal multiscale methods). For the length scale extension, atomistic-continuum coupling may be applied: atomistic method can be used in the region of localized deformation where resolution down to atomic scale is desired (fine scale domain), and continuum method (e.g., FEM) everywhere else (coarse scale domain).

    In theory, the modulus estimation of deformation involves different types of processes, including elastic, plastic, and fracture theories, and its theoretical mechanism is very complex [6]. Therefore, the prediction of hardness and mechanical modulus of deformation has a lot of difficulties in theory. So, the experimental measurements to calculate elastic modulus become the easiest way, mainly a variety of indentation experiments [7]. That is why, the nanoindentation experiment has been widely used to characterize mineral and metal materials [8].

    Nanoindentation technology has become an important means to measure material properties and is widely used in the fields of medical devices, microelectro mechanical systems (MEMS), bioengineering, civil engineering, aircraft, and film material research. Indentation technology is the expansion of conventional hardness test to microscale (no more than submillimeter) [9], the principle of it is to get parameters of microindentation hardness and elasticity modulus, etc., according to loading/unloading curves between load and displacement [10, 11]. The disadvantage of the classical plasticity theory is that it has not considered size effect [12, 13]. Models commonly used during micro- and nanoindentation process mainly include elastic recovery model, Kick model [14], proportioned sample model (PSR), Hays-Kendall model, and Taylor dislocation model [15, 16]. In general, the elastic moduli of covalent crystals are intrinsic and equivalent to the sum of resistance of each bond per unit area to indenter. In this chapter, analytical method is used to solve the nanoindentation simulation with indentation depth during nanoindentation process and fitting methods are adopted to establish the finite-element constitutive model.

    Using backscattered images and point analysis, outer and inner C▬S▬H can be distinguished. Figure 1 shows SEM BSE micrograph of different phases in ordinary Portland cement (OPC) paste at 7 days of hydration [17] and point’s position of the paste with w/c = 0.40 after 28 days of hydration for the EDX-SEM analysis [17].

    Figure 1.

    BSE micrograph of typical phases and the C▬S▬H points in cement paste. (a) SEM BSE micrograph of an OPC paste [17]. (b) BSE micrograph map of the inner and outer C▬S▬H points [17]. (c) SEM BSE micrograph of the OPC paste [17]. (d) Regional distribution of chemical elements [17].

    From Figure 1(a, b), it is indicated that the outer C▬S▬H has a brighter gray level, while inner C▬S▬H is mainly localized within the original cement grains. As is shown in Figure 1(c, d), there are mainly four regions of hydrated products, such as: voids region (the cumulative hydration voids), C▬S▬H region (including low density and high density), CH phase region, and RC phases region (the residual cement in a cement paste).

    Elastic moduli and hardness can be used to understand the extent, to which a solid resists both elastic and plastic deformation. C▬S▬H gel has high surface areas and shows excellent adhesive characteristics. Nanoindentation can be used to understand differences in hydration products that result from water availability for clinker during hydration. Ye et al. [18] have simulated the microstructure of cementitious materials by using the cement hydration model-HYMOSTRUC. Ulm and coworkers [19, 20] have provided the experimental techniques and analysis tools to provide reliable, mechanical technical tools to identify the relative volume fractions and spatial distributions of cement hydration products in cement paste samples produced with conventional cements. Besides, they have illustrated that there are three types of C▬S▬H (LD C▬S▬H, HD C▬S▬H, C▬S▬H/CH nanocomposite) in hardened cement paste. The average nanoscale material responses (e.g., indentation modulus and indentation hardness) for each of these three C▬S▬H phases are virtually the same for all Portland cement pastes tested to date [21, 22].

    Nanoindentation method may provide a wide variety of mechanical characteristics without substrate effects to obtain information on the mechanical properties of hard, soft, brittle, and ductile materials at the nanometer scale. The techniques have been widely used in the characterization of small-scale mechanical behavior of materials, which provides an ideal approach for measuring phases at nanoscale upon a substrate [23, 24]. Although there are many reports on nanoindentation, however, these studies are almost based on experimental analysis and its numerical calculation at nanoscale has not yet been revealed. Multiscale modeling and simulation of atomistic/continuum coupling in computational CH, LD C▬S▬H, and HD C▬S▬H structures with the size of nanometer unit during nonlinear indentation is considered and simulated to compare the results of Keinde [25] and Constantinides [26]. In this chapter, nanoindentation simulations are carried out to determine elastic moduli of several phases in cement pastes, compared with results of nanoindentation experiments. As the volume fraction of main components of CH and C▬S▬H in hydrated Portland cement paste separately reach more than 14 and 40%, our main task is to predict the Young’s modulus via multiscale simulation and nanoindentation experiment.

    2. Background of theory, principle, and measurement of nanoindentation

    2.1 Numerical description of nanoindentation procedure

    During indentation, the load and displacement were continuously monitored by a three-plate capacitive force/displacement transducer. Nanoindentation hardness is defined as the maximum indentation load Pmax divided by the contact area Ac:

    H=PmaxAchc(1)

    The indenter shape function A (nm2) can be seen as a function of the indentation depth hc:

    A=fhc(2)

    Here, area function selected depends on the function of Oliver and Pharr proposed and modified considering loading system flexibility and curvature radius of the indenter tip correction. So, the projected contact area calculation is as:

    Ap=C0hc2+C1hc+C2hc1/3+C3hc1/7+(3)

    where C1, C2Cn are coefficients describing geometrical shape of indenter, of which high major is used to describe deviation degree of relative ideal shape at indenter tip.

    The unloading stiffness S is then established by differentiating P at the maximum depth of penetration. S is the initial unloading contact stiffness, which is represented by the slope of the initial portion of the unloading curve and described as:

    S=dPdhh=hmax(4)

    The contact area as a function of contact depth hc is calibrated by indenting on a standard cement specimen. The contact depth hc can be estimated from the load-displacement data using Sneddon’s equation [27] and assuming that pile-up is negligible:

    hc=hmaxεPmaxS(5)

    where ε is a constant that depends on the indenter geometry; ε = 0.75 for the Berkovich indenter, ε = 0.72 for the conical indenter, and ε = 1 for the flat punch. Pmax is the maximum load and S is the slope of unloading curve.

    The indentation depth hmaxat any time during loading can be described as follows:

    hmax=hc+hs=hc+εPmaxS(6)

    where hs is the displacement (nm) over the initial surface nearby the section between indenter and tested sample. The reduced modulus of material is defined as:

    Er=π2βSA(7)

    which can be further expressed as:

    1Er=1ν2E+1νi2Ei(8)

    where hc is the contact depth, Es and vs are the elastic modulus and Poisson’s ratio of the specimen, and Ei and μi are the elastic modulus and Poisson’s ratio of the indenter. For a diamond indenter, Ei = 1140 GPa and μi = 0.07. The unloading curve defined by Oliver and Pharr is described by a power law relation:

    P=Ahhfm(9)

    where A and m are fitting parameters empirically determined.

    It can also be noted that Young’s modulus E and Poisson coefficient μ are linked in Eq. (8), and without any assumption on the value of one of these characteristics, the other one cannot be determined by the nanoindentation test.

    2.2 Principle of the finite-element analytical method

    Guillonneau [28] has discussed the process of nanoindentation in detail both the measurement and the numerical simulation, where FEM method is used to know the deformation distribution in the contact area.

    Figure 2 is a schematic diagram of rigid conical indenter and its geometrically necessary dislocations underneath. Assuming that the indenter is rigid, microplastic deformation caused by geometric necessary dislocations underneath illustrated in Figure 2 can be expressed as follows:

    tanθ=h/a=b/s,s=ba/h(10)

    Figure 2.

    Schematic diagram of rigid conical indenter and its dislocations region underneath [15]. (a) Different schemas and different types of indentation points [28]. (b) Schematic diagram of rigid conical indenter and its dislocations region underneath [15].

    where θ is the angle between indenter surface and sample surface, a is indenter radius contacted, h is indentation depth, b stands for Burgers vector, and s represents distance within a slipped step of indentation.

    Huang et al. [29] have developed the model of Nix and Gao [30] and proposed a new theoretical model based on the maximum geometric necessary dislocation density allowed. Assuming that λ is the total length of dislocation loops, integration from r to r + dr within dislocation loops is carried out in the whole contacted radius a of indenter and λ can be described as follows:

    λ=0a2πrsdr=0ahba2πrdr=hbaπr20a=habπ(11)

    Considering size effect, the intrinsic material length is introduced. Assuming that flow stress is under control of microscale dislocation motion and compliances with Taylor hardening relation, the equivalent stress in microscale can be expressed as: σ¯=σyf2ε¯+. The model of Nix and Gao et al. [30] on the basis of continuum mechanics can be described as follows:

    σ=Kεp2n+=Kεp2n+lρGb(12)

    where σrepresents stress, K is strengthening coefficient, n is strain hardening exponent, b is Burgers vector, ρGis geometric dislocation density, εpis the peak of effective strain, ηis equivalent strain gradient, expressed as η=ρGb, and l is intrinsic material length (microns) and can be expressed as: l=Mαμ/σref2b.

    Assuming that all dislocation loops are within the volume V of hemisphere with its radius a, the geometrically necessary dislocation density caused by the ideal conical rigid indenter with the sample being pressed can be described as follows:

    ρG=λV=ahπ/b2a3π/3=3h2a2b=32bhtan2θ(13)

    where h is the indenter depth, and θ is the angle between indenter surface and sample surface with its value 19.7°.

    Assuming that conical indenter is rigid and the contact surface with indenter is under friction-free state, the schematic diagrams of load state and displacement state brought by conical indenter with material pressed can be expressed in Figure 3.

    Figure 3.

    Schematic diagrams of displacement and load during microindentation. (a) Conical indenter. (b) Spherical indenter.

    The indentation schematic diagrams using Berkovich indenter are shown in Figure 3(a). Take the spherical indenter for example, as is illustrated in Figure 3(b), the displacement boundary on contact sphere and the contact boundary of interfacial nodes at the top of indenter must keep the sliding condition at spherical part. Moreover, the displacement boundary on conical surface and the contact boundary of interfacial nodes at conical part of indenter must keep the sliding condition too. So, the displacement boundary of node i at spherical part and the incremental displacement of node j at conical part under the incremental calculation step n during microindentation process can be described as follows:

    ria2ri2Δui+Δvi=δn(14)
    Δuj/tgθ+Δvj=δn(15)

    where r and z are cylindrical coordinates, (ri, zi) stands for coordinate of node i; Δui and Δvi are incremental node displacements in r and z directions, a is the hemispherical radius of indenter, δnis incremental displacement caused by indenter with material being pressed, θis half-cone angle, and the value is 70.3° for Berkovich indenter.

    Assuming that conical indenter is fully rigid and there is no friction between indenter and material [31], as the contact surface stress of sample pressed by indenter is at normal direction, so the load of node can be adopted approximately at normal direction. The nodal force at node N can be divided into FNr and FNz at directions of r and z, so the incremental loads at directions of r and z can be separately expressed as ΔFNr and ΔFNz. The incremental loads of spherical contact node i and conical contact node j can be separately described as follows:

    ΔFircosα+ΔFizcosα=0(16)
    ΔFjr+ΔFjztgθ=0(17)

    The modified incremental load ΔFemust be added in the direction of radius to keep the deformation load along the direction of normal at contact surface required, the expression of ΔFecan be described as follows:

    ΔFe=FirFiztgα=0(18)

    where exists the equation: α=arctgri/a2ri2.

    For the regular resolution of contacted fine-grid section, the incremental finite element equation can be described as:

    KΔd=ΔF(19)

    where [K] is total stiffness matrix, Δd and ΔF, respectively, represent incremental displacement and load. Insert load boundary conditions and displacement boundary conditions of corresponding contact surface nodes into Eq. (19) and Eqs. (14), (16), and (18) are inserted for the ball contact surface node i while the proportionally conical surface contact node j with Eqs. (15) and (17) inserted, so the microindentation hardness obtained ultimately can be resolved by the simultaneous equations described as:

    k2i1,1+k2i,1tgαk2i1,2i1+k2i,2i1tgαk2i1,2i+k2i,2itgα0ri/a2ri21k2j1,1+k2j,1/tgθk2j1,2j1+k2j,2j1/tgθk2j1,2j+k2j,2j/tgθ01/tgθ1×ΔuiΔviΔujΔvj=FirFiztgαδn0δn(20)

    Solving equations of the contact nodes in Eq. (14) as well as noncontact nodes treated by conventional computing method, displacement distribution during unloading from Pmax can be finally obtained.

    Assuming that the material follow Mises flow law and the Tabor factor [32] is equal to 3, the microindentation hardness of material can be converted by equivalent stress described by shear strength with the dislocation density [32] according to the model of Nix and Gao [30] described the Taylor hardening relationship by dislocation density, so the equation of microindentation hardness can be described as follows:

    H=3σ¯=33τ=33αμbρS+ρG(21)

    where ρSis the statistics stored dislocation density, μis the shear modulus, and αis an empirical constant, whose magnitude order value is one.

    As experimental microindentation hardness values are related with maximum load Pmax (mN) and indenter shape function A (nm2), the relationship between relative microindentation hardness values and indentation depth considering the intrinsic material length l without considering the effect of indentation elastic recovery can be expressed as follows:

    HH0=1+hh=1+khd(22)

    where H is indentation hardness (GPa), H0 is indentation hardness without considering strain gradient effect, h is indentation depth, d is the reciprocal of indentation diagonal length, k is identified as 4tanθ/2, and A is area function.

    The expressions of A, H, H0, and h were separately described as:

    H0=33αμbρS(23)
    h=812bα2tan2αμH0(24)
    H=Pmax/A(25)
    A=Pmax/n=0SCnhc1/2n1(26)

    In case of analytical extraction of elastic properties through nanoindentation data, it is possible to assess roughly the yield stress as being directly correlated with hardness value according to Tabor [33], and given by:

    σyH/C(27)

    where H is indentation hardness (GPa), and C is the contact creep modulus (GPa); the coefficient can change in the range of 2.6–3.0.

    3. Nanoindentation modeling and simulation at nano-/microscale

    3.1 Nanoindentation modeling using the QC method

    3.1.1 Representative volume element (RVE)

    In general, 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 [34]. 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 [35] and MuMoCC platform [35, 36] are developed. The representative volume element primarily uses numerical approximations of the effective properties of composite. Square or cubic RVEs are used for most numerical approximations because of the ease of numerically solving boundary value problems with these geometries. The difficulties involved in generating statistical information about particle distributions and concentrations lead to difficulties in the rigorous determination of RVE sizes. Sab [37] has shown that if an RVE exists for a random composite material, the homogenized properties of the material can be calculated by the simulation of one single realization of the medium. The “ergodic” hypothesis assumes that the ensemble average is equal to the volume average. The ensemble average is the mean of a large number of realizations of the microstructure. The volume average, on the other hand, is the average as the RVE volume becomes infinitely large compared to the volume of a particle.

    The MuMoCC platform operates in two main steps. The first one consists of generating realistic 3-D representative volume elements of concrete and their three different scales (cement paste, mortar, and concrete). This hierarchical modeling is based on the assumption of scale separation which is for the moment necessary. Indeed, because of the wide range of feature size in concrete, it is impossible to represent concurrently all these structural features in a single mesomodel and then to represent the complete system as one. The NIST 3-D model, CEMHYD3D, is then used for this purpose. The MuMoCC platform enables one to couple this well-known model to the FE software ABAQUS. A set of algorithms has been developed in order to convert the voxelized images obtained in this first step by CEMHYD3D to FE meshes compatible with ABAQUS. Figure 4 shows a 2-D cross section of the generated submesostructure (mortar) and the modeled microstructure (hydrated cement paste).

    Figure 4.

    Structural modeling of the meso- and microstructure models [38, 39]. (a) Mesostructure model of mortar (125 × 125 μm); (b) microstructure model of the cement paste (1 × 1 μm).

    From Figure 4(a), we can see that, at the mesoscale of mortar, large aggregates with interfacial transition zone (ITZ) embedded in a matrix of mortar can typically be observed. From Figure 4(b), at microscale, the hydrated cement paste corresponding to the considered mortar is modeled using NIST’s CEMHYD3D code. As also seen before, this paste is a heterogeneous and complex porous medium in which the main solid phases are calcium silicate hydrate (C▬S▬H), Portlandite (calcium hydroxides) and hydrated aluminates, and sulfoaluminate phases. The solid phases are in chemical equilibrium with an interstitial solution that partially or totally fills the porosity. The second step of the method consists of imposing various boundary conditions on the meshed RVE. The failure in mortar is supposed to be the consequence of crack initiation and propagation at highly localized regions under large tensile stress concentration. Failure in compression is supposed to be the consequence of tensile stresses resulting from restrained deformations due to the heterogeneous mesostructure of mortar. Thus, the tensile behavior of the hardened cement paste is needed as input data to bridge the scales.

    Thus, the quasicontinuum (QC) method is a mixed continuum and atomistic approach for simulating the mechanical response of polycrystalline materials. For nanoindentation modeling using the QC method, the representative volume element (RVE) is used for nanoindentation simulation.

    3.1.2 The 2-D/3-D modeling and axisymmetric finite element simulations

    The main eight steps of nanoindentation simulation using ABAQUS finite element analysis software are mainly divided as:

    1. ABAQUS CAE component module: A solid model composed by indenter and the sample is established first. The sample is established as a 2-D/3-D axisymmetric deformable part, and the indenter is thought as rigid since the hardness and elastic modulus of the indenter material is far greater than the sample material. For Berkovich indenter model, the half-cone model with a degree of 70.32° is used, thus has the same relation between area and penetration of nanoindentation experiment. In order to eliminate the influence of outside boundary conditions on result of the indentation region, the sample length of each side is set 40 times the indentation depth.

    2. Characteristics module: Material properties (elastic constants, plastic, etc.) are defined.

    3. Assembly module: Model instances of indenter and sample are defined and assembled.

    4. Grid module: The samples are meshed and defined by certain element types. In order to better simulate the deformation behavior of the indentation affected zone, setting with a bias grid, unit type is quadrilateral reduced integration unit with the type of CPE4R element.

    5. Analysis step module: The unloading step, time step, and analysis step are defined.

    6. Interaction module: Related properties can be defined. In this chapter, indenter surface is set as the main surface, and the surface of sample is subjected to the main surface. The friction coefficient between contact surfaces may be defined.

    7. Load module: The bottom surface of the sample is fixed, the left axis is of symmetry, and boundary conditions are applied to the contacted indenter point.

    8. Analysis job module: The job is established and then submitted to analysis solver.

    The nanoindentation constitutive model with embedded programs was used to simulate nanoindentation process by finite element analysis. The meshing diagram is shown in Figure 5.

    Figure 5.

    Indentation simulation on 2-D/3-D models of typical cement paste phases. (a) 2-D geometrical details and finite element mesh. (b) 3-D geometrical details and finite element mesh.

    Based on the steps mentioned above, 2-D/3-D axisymmetric finite element simulations are performed to investigate the effect of friction on the indentation response of phases, e.g., C▬S▬H [38]. The indenter was modeled as a rigid cone with half apex angle of 70.32°. Details of the model geometry are shown in Figure 5. This conical angle ensures the same contact depth-projected area of contact relation as in Vickers and Berkovich pyramidal indenters.

    3.2 Typical P-h curve analysis of the elastic modulus

    The penetration depth is measured in nanometers continuously during the loading process in nanoindentation simulation. Apart from the scale movement, the distinguishing feature of most nanoindentation testing is the indirect measurement of the contact area (real area where contact is made indenter). In conventional indentation test, the contact area is calculated by measuring the size of the residual footprint on the surface of the sample after unloading. In the nanoindentation test, the contact surface is determined so indirectly from the geometry of the indenter and the depth of penetration; meanwhile, the force and the displacement of the indentor in the material are continuously recorded and shown in a force-depth curve, as in Figure 6.

    Figure 6.

    Typical load-depth curves in nanoindentation experiment.

    As is shown in Figure 6, S is the initial unloading contact stiffness, hp is the residual depth during plastic deformation, hc is the contact depth, Pmax is the maximum of loading stress, and hmax is the maximum of indentation depth. The typical indentation curve in Figure 6 represents a loading-unloading cycle, not superimposed when the test leaves a permanent imprint on the surface. The first part is the loading curve corresponding to the penetration of the indentor. This ends when the indenter reached its maximum load Pmax or the maximum penetration hmax. The second part represents the discharge curve. It corresponds to the withdrawal of the indenter. The analysis of a cycle gives a lot of information by reading the indentation curve total penetration depth hmax and depth of the residual depth hp. The area enclosed by load curve, unload curve, and horizontal axis represents plastic deformation work, while the area enclosed by unload curve, horizontal axis, and the vertical ordinate of maximum load represents elastic deformation work.

    During the nanoindentation simulation, the movement can apply forces from micronewtons up to several millinewtons and measurement depths down to a few nanometers. To achieve such good resolution at nanometer length, the model should be generally placed in a soundproof enclosure on an antivibration table to eliminate the negative influence of vibrations and noise so as to obtain the correspondingly accurate load-depth curves which we need.

    4. Nanoindentation technique for CH and C▬S▬H structures by FEM

    Nanoindentation method is a relatively simple and effective technique to evaluate the mechanical properties of material, and it can not only obtain the relevant performance parameters of the material but also reflect the mechanism of elastic-plastic behaviors and reveal relationship between microstructure and macroscopic mechanical properties. However, it has a complex problem of contact, which is affected by many factors, such as nanoindentation test, the surface roughness, the substrate effect, the grain boundary effects, indenter geometry, lattice anisotropy, size effect, etc. Even if the nanoindentation process is tested in the same device and under the same test conditions, the result cannot be guaranteed repeatability. That is why, the investigation of nanoindentation process research using numerical simulation method is required.

    4.1 Models of nanoindentation on 2-D/3-D models

    The nanoindentation constitutive model with embedded programs was used to simulate nanoindentation process by FE analysis. At the mesoscale, subdomains are regarded with locally varying loads, and interaction of the constituents is accounted for at the microscale [39]. For simplicity purposes, the Berkovich indenter is commonly modeled as a conical indenter with a semiapex angle of 70.32°, where the conical cross-sectional area closest to the actual situation [40]. According to literatures [41, 42], 3-D models can be equivalent to an axisymmetric 2-D model and the results of two models are consistent with each other. Thus, the nano-indentation 3-D model is treated as an axial symmetrical 2-D model, shown in Figure 7.

    Figure 7.

    Indentation simulation on 3-D and 2-D models of typical cement paste phases. (a) 3-D geometrical model. (b) 2-D model with three mesh regions. (c) Partial 2-D model with fine mesh.

    From Figure 7(a), the 3-D model has a size of 15.0 × 15.0 × 7.5 μm, which has the sufficient space so as to compare the experimental results. The CPE4R elements with higher accuracy and axisymmetric model are adopted to simulate nanoindentation process. Berkovich indenter has an angle of 70.32°. To simplify the calculation, Berkovich indenter is assumed to be rigid and the deformation is subjected to the Mises yield criterion and isotropic hardening criterion. The size of microindentation model regarded as the semiinfinite solid in Figure 7(b) is 12.0 × 12.0 μm. This continuum space is discretized using four-noded axisymmetric, isoparametric element (CPE4R element). A mesh sensitivity analysis was performed to ensure that the simulation results were insensitive to the mesh size in the indenter region. To save the computational time, the transitional gird division method is adopted. The element size was continuously refined as approaching the indenter contact region for greater accuracy, a refined mesh in Figure 7(c) was near the tip area (contact zone is 3.0 × 3.0 μm, the number of grid is 300 × 300) in order to minimize the computing resource and to maintain the reasonable accuracy. The minimum length of elements within the final domain is about 10 nm.

    4.2 Parameters of typical structures in simulation

    Correspondingly, the indentation is simulated by ABAQUS software at the same depth from the simulation. The applied displacement force is added. Mechanical parameters of typical phases in simulation are found in the previous chapters and listed in Tables 1 and 2.

    Phases/parametersLD C▬S▬HHD C▬S▬HCH
    E (GPa)17.3532.7039.88
    Friction coefficient f0.60.60.6
    μ0.2–0.30.2–0.30.2–0.3

    Table 1.

    Mechanical properties and parameters of typical phases and structures used in simulation: parameters of the representative experimental CH and C▬S▬H phases.

    Phases/parametersLD C▬S▬HHD C▬S▬HCH
    S-C schemeM-T schemeS-C schemeM-T schemeAFEMDFT + RVH
    Elastic modulus E (GPa)18.2229.7831.6337.7143.1345.46
    Poisson’s ratio μ0.250.250.250.250.250.20
    Friction coefficient f0.60.60.60.60.60.6

    Table 2.

    Mechanical properties and parameters of typical phases and structures used in simulation: parameters of the CH, LD C▬S▬H, and HD C▬S▬H structures using various simulation techniques.

    In the following simulations, as in the reference of Asroun et al. [43], an elastic perfectly plastic material model with Von Mises yield criterion σy (GPa) is assumed for the indented materials. The yield stresses are considered to be equal to 0.198 GPa for LD C▬S▬H, 0.347 GPa for HD C▬S▬H, and 0.440 GPa for CH.

    Frictional effects in the indenter-material interface were included in the analysis through an isotropic model. We assume that the loading rate is slow enough such as static friction can securely model the interface response. Simulations are proceeded in two steps. The indenter was first subjected to a ramped vertical displacement, followed by an indenter retraction to the original position which corresponded to complete unloading at zero load. Moreover, only half of the material is chosen to simulate nanoindentation with fraction coefficients between indenter and material adopted. Axisymmetric boundary conditions were used along the symmetry axis beneath the indenter. During simulation, five degrees of freedom are limited except y-axis and the bottom plane (constrained vertically) and planes around are fixed, then material model is defined. The contact mode is used during loading and unloading processes.

    For comparative analysis, the FEM simulation of nanoindentation on the LD C▬S▬H, HD C▬S▬H, and CH structures using parameters in Tables 1 and 2 are carried out using ABAQUS.

    5. Discussion on elastic moduli by both experiment and simulation

    5.1 Experiment and discussion of elastic moduli of CH and C▬S▬H

    5.1.1 Nanoindentation experiment and regional indentation analysis

    After the sample preparation and chemical composition, nanoindentation experiment is done on the instrument of INSA de Rennes, and CSM instrument and sample are shown in Figure 8.

    Figure 8.

    CSM instrument and cement paste used in experiment. (a) Nanohardness platform. (b) Berkovich indenter. (c) Cement paste sample.

    As in Figure 8(a), CSM instruments include: ultrananoindentation tester head, nanoindentation tester head, atomic force microscopy, and optical video microscope. Besides, Berkovich indenter is shown in Figure 8(b). Samples with the dimension 1.1 × 1.1 × 1.1 cm underwent mechanical polishing and electrolytic polishing separately and lapped to a mirror, as shown in Figure 8(c). Parameters are as: the maximum load is 1.5 mN; velocity after contact with the sample surface is 10 m/s; continuous maximum load is 5 s; unloading time is 10 s. The speed of loading/unloading was 9.6841 mN/s and the P-h data were recorded. The test is repeated several times to average.

    5.1.2 Elastic moduli of typical phases in various samples

    Nanoindentation experiments were carried out on CSM microindentation tester with conical Berkovich indenter. The process is mainly described as follows: samples with the dimensions 1.1 × 1.1 × 1.1 cm underwent mechanical polishing and electrolytic polishing separately and lapped to a mirror. Then, experimental maximum loads of 1.5 and 2 mN were mainly selected by loading/unloading with the speed of 9.6841 mN/s and the data of load and depth were recorded. Besides, indentation diagonal lengths were measured by data collection software on computer. Under each parameter, experiments were repeatedly done five times and then averaged to decrease the relative error. By experimentally investigating the mechanical properties of cement paste at different length scales provides a means of correlating such nanoscale properties to nanoscale applications. The identification of the parameters is provided in Table 3.

    Experimental parametersUnitCEMICEMIIICEMI_MKCEMI_Sed
    Maximum load PmaxmN1.5–2.01.51.51.5–2.0
    Strain rate P*/Pmaxs−10.010.010.010.01
    Loading and unloading speedmN/min3–43–42–52–5
    Acquisition frequency rateHz10101010
    Delta slope (%)80808080
    Poisson’s ratio0.20.20.20.2
    Total number of indents830363334296

    Table 3.

    Characteristics of experimental parameters used in loading programs for samples.

    For the tested curves of CH phase collected from different samples, the parameters of CH phase obtained microindentation experiment are shown in Table 4.

    Experimental parametersLD C▬S▬HHD C▬S▬HCHCalciteClinker
    Maximum depth hmax (nm)441.79278.08212.9897.1199.54
    The projected contact area Ap (nm2)4577259.841761371.10977898.05137943.39267721.30
    The unloading stiffness S (mN/nm)0.05320.05590.04580.03630.0588
    The fitting parameter m4.766.014.121.291.66
    Indent number (e.g., CEMIII, etc.)30225613

    Table 4.

    Experimental parameters and corresponding characteristics of typical phases.

    The representative experimental curves of typical phases characterized by individual elastic modulus are averaged. Using the indents of LD C▬S▬H phase under different samples, the corresponding P-h curves are found. The representative experimental curves of outer LD C▬S▬H phase are given in Figure 9.

    Figure 9.

    The representative experimental curves of outer LD C▬S▬H phase. (a) Typical curves of outer LD C▬S▬H phase. (b) The representative curve of outer LD C▬S▬H phase.

    Figure 9 shows the representative experimental curve of LD C▬S▬H phase, with elastic modulus of 17.35 GPa. Similarly, representative experimental curves of HD C▬S▬H phase is shown in Figure 10.

    Figure 10.

    The representative experimental curves of inner HD C▬S▬H phase. (a) Typical curves of inner HD C▬S▬H phase. (b) The representative curve of inner HD C▬S▬H phase.

    Figure 10 shows the representative experimental curve of HD C▬S▬H phase, with elastic modulus of 32.70 GPa. Similarly, the representative corresponding P-h curve of CH phases is given in Figure 11.

    Figure 11.

    The representative experimental curves of CH phase. (a) Typical curves of CH phase. (b) The representative curve of CH phase.

    Figure 11(b) shows the representative experimental curve of CH phase, with its corresponding elastic modulus of 39.88 GPa.

    Overall, nanoindentation experiments show that the range of elastic modulus corresponding to low-density C▬S▬H is 15–26 GPa and the range of elastic moduli corresponding to high-density C▬S▬H is 26–39 GPa. The inverse method described above is applied to each matrix phase that is what has been called previously: LD C▬S▬H, HD C▬S▬H, portlandite (CH), and clinker.

    5.2 Simulation and discussion of elastic moduli of CH and C▬S▬H

    There are two methods of applying load in simulation: forced loading and forced displacement. Comparison of the load-displacement curves using the forced loading and the forced displacement is the same; here, we use the forced displacement method.

    5.2.1 Simulation and discussion of elastic modulus of CH phase

    From the load-depth curve of CH phase in nanoindentation experiment, the pressure is about 1.5 mN and the maximum indent depth of CH phase is about 212.98 nm, the representative curve of CH is as in Ref. [44]. CH simulation parameters were separately set as: elastic moduli E are separately 45.459 GPa (DFT + RVH method, μ = 0.2), 43.13 GPa (AFEM method, μ = 0.25), and 39.880 GPa (nanoindentation experiment). The region range of Poisson ratio μ is 0.25–0.30. The yield stress can be averaged by the experimental data [21] and the simulated data [38], which is 0.440 GPa by the Eq. (27). The contact friction is set about 0.6 [26].

    Stress distribution and deformed displacement zones of CH are shown in Figure 12.

    Figure 12.

    Stress distribution and displacement zones of CH under various frictions. (a) Stress distribution zone (f = 0.6). (b) Deformed displacement zone (f = 0.6).

    As is shown in Figure 12, the maximum residual stress after unloading on the ideal contact condition distributes just below the indenter, while the maximum residual stress after unloading on the friction condition distributes around the top of indenter and the residual stress below the indenter is not too large.

    Using ABAQUS software, under the loading force of 1.5 mN, the comparison of experimental and simulated unloading P-h curves of CH phase is shown in Figure 13.

    Figure 13.

    Comparison of experimental and simulated P-h and unloading curves of CH phase. (a) Comparison of experiment and AFEM and DFT methods. (b) Unloading curves under various Poisson’s ratio.

    As is shown in Figure 13, we can see that the simulated unloading part is close to the representative experimental part. It can be seen that the smaller Poisson’s ratio is, the smaller the slope of unloading curve will be. Comparison of purely elastic unloading curve of CH phase is shown in Figure 13; we can see that DFT methods (45.459 GPa, by LDA method) and AFEM method (43.13 GPa) are very consistent with experimental results (39.880 GPa), which verifies the reliability of simulation method.

    For CH structure in this book, the elastic constants have been calculated using the AFEM and DFT methods. Then, the homogenized elastic modulus of CH structure is estimated by the Y parameter, elastic moduli by DFT method (45.459 GPa) and by AFEM (43.14 GPa) are very close to the results of Constantinides. In order to verify the simulation results of the CH structure, nanoindentation experiments of four cement samples are carried out. The averaged Young’s modulus of CH phase is 39.880 GPa by nanoindentation experiment, which is close to the results of 40.3 ± 4 GPa by Constantinides. Comparison of purely elastic unloading curve of CH phase shows that DFT and AFEM methods are very consistent with experimental results on average, which verifies the reliability of both simulation methods. It has demonstrated the reasonableness of simulation conclusions by the reliable experiment verification.

    The simulated elastic modulus is close to the average by experiment, which has demonstrated the reasonableness of simulation conclusions. Although the mechanism of the complex plastic deformation is still unknown, the yield flow stress parameter is still turned to be of certain rationality. The nanoindentation simulation of CH phase is done on ABAQUS using the corresponding parameters, thus reversely verify the reliability of the initial parameters. The large-scale models using QC methods are considered and verified on ABAQUS software, where, on behalf of the purely elastic stage, the unloading curves of the CH model with various Poisson’s ratios are simulated. It can be seen that the smaller Poisson’s ratio is, the smaller the linear slope of unloading curve will be.

    5.2.2 Simulation and discussion of elastic modulus of LD C▬S▬H

    Firstly, elastic constants of the monoclinic 11 Å tobermorite are calculated by DFT. Then, the homogenized elastic modulus is determined to be 79.51 GPa, which is close to the other references. However, the ordered 11 Å tobermorite model cannot reflect the C▬S▬H structure, and the crystalline C▬S▬H structures—C▬S▬H(I) (Ca/Si = 0.67) and the C▬S▬H(II) (Ca/Si = 1.67)—are still far from the data in laboratory. That is why, the disordered structure is simulated and confirmed. For the monolithic structures of C▬S▬H(I) and C▬S▬H(II), using MD simulation, the averaged Young’s moduli are separately about 75.92 and 60.95 GPa. The last value of C▬S▬H(II) with the size about 5 nm enables to evaluate the value of the elastic modulus of LD C▬S▬H and HD C▬S▬H: 18.11 and 31.45 GPa, using the self-consistent scheme.

    The LD C▬S▬H phase can be modeled assuming a linear elastic and cohesive-frictional plastic material. The forced displacement method is used as mentioned above.

    From the load-depth curve of LD C▬S▬H phase nanoindentation experiment, the pressure is about 1.5 mN and the maximum indent depth is about 441.79 nm, the representative curve of LD C▬S▬H is in Ref. [44]. The yield stress can be averaged by the experimental data [21] and the simulation [38], which is 0.198 GPa. The LD C▬S▬H parameters were separately set as follows: elastic moduli E are separately 18.22 GPa (by the self-consistent scheme, μ = 0.25), 29.78 GPa (by the Mori-Tanaka scheme, μ = 0.25), and 17.35 GPa (nanoindentation experiment). The region range of Poisson ratio μ is 0.20–0.25 and the contact friction is about 0.6 [26].

    Stress distribution and deformed displacement zones of LD C▬S▬H are shown in Figure 14.

    Figure 14.

    Stress distribution and displacement zones of LD C▬S▬H under various frictions. (a) Stress distribution zone (f = 0.6). (b) Deformed displacement zone (f = 0.6).

    As is shown in Figure 14, the maximum residual stress after unloading on the ideal contact condition distributes just below the indenter, while the maximum residual stress after unloading on the friction condition distributes around the top of indenter and the residual stress below the indenter is not too large. Similarly, the comparison of experimental and simulated unloading P-h curves of LD C▬S▬H phase is shown in Figure 15.

    Figure 15.

    Comparison of experimental and simulated unloading P-h curves of LD C▬S▬H phase. (a) Comparison of experiment and MD simulation. (b) Unloading curves under various Poisson’s ratio.

    As is shown in Figure 15(a), we can see that the simulated unloading part is close to the representative experimental part. It can be seen that the smaller Poisson’s ratio is, the smaller the slope of unloading curve will be. Comparison of purely elastic unloading curve of LD C▬S▬H phase is shown in Figure 15(b); we can see that, the Mori-Tanaka scheme (29.78 GPa, μ = 0.25) and the self-consistent scheme (18.22 GPa, μ = 0.25) are very consistent with experimental results (17.35 GPa), which verifies the reliability of both simulation methods. Comparison of purely elastic unloading curve of LD C▬S▬H phase shows that the simulated P-h curve (18.22 GPa, by the self-consistent scheme of C▬S▬H(II) structure (C/S = 1.67) with the size about 5 nm simulated using MD method in Chapter 5) is close to the experimental curve (17.35 GPa), which verifies the reliability of both simulation methods.

    5.2.3 Simulation and discussion of elastic modulus of HD C▬S▬H

    From the representative load-depth curve of HD C▬S▬H phase nanoindentation experiment, the pressure is about 1.5 mN and the maximum indent depth of HD C▬S▬H phase is about 278.08 nm, the representative curve of HD C▬S▬H is in Ref. [44]. Similarly, the HD C▬S▬H parameters were separately set as follows: elastic moduli E are separately 31.63 GPa (by the self-consistent scheme, μ = 0.25), 37.71 GPa (by the Mori-Tanaka scheme, μ = 0.25), and 32.70 GPa (nanoindentation experiment). The region range of Poisson ratio μ is 0.20–0.25, and the yield stress in the Eq. (27) is 0.347 GPa [38]. By the Ref. [26], the contact friction is about 0.6.

    Stress distribution and deformed displacement zones of HD C▬S▬H are in Figure 16.

    Figure 16.

    Stress distribution and displacement zones of HD C▬S▬H under various frictions. (a) Stress distribution zone (f = 0.6). (b) Deformed displacement zone (f = 0.6).

    As is shown in Figure 16, the maximum residual stress after unloading on the ideal contact condition distributes just below the indenter, while the maximum residual stress after unloading on the friction condition distributes around the top of indenter and the residual stress below the indenter is not too large.

    Using ABAQUS software, under the loading force of 1.5 mN, the comparison of experimental and simulated unloading P-h curves of HD C▬S▬H phase is shown in Figure 17.

    Figure 17.

    Comparison of experimental and simulated unloading P-h curves of HD C▬S▬H phase. (a) Comparison of experiment and MD simulation. (b) Unloading curves under various Poisson’s ratio.

    As is shown in Figure 17(a), we can see that the simulated unloading part is close to the representative experimental part. It can be seen that the smaller Poisson’s ratio is, the smaller the slope of unloading curve will be. Comparison of purely elastic unloading curve of HD C▬S▬H phase is shown in Figure 17(b); we can see that the Mori-Tanaka scheme (37.71 GPa, μ = 0.25) and the self-consistent scheme (31.36 GPa, μ = 0.25) are very consistent with experimental results (32.70 GPa), which verifies the reliability of simulation method. This indicates that the selected simulation parameters are reliable to characterize its plastic deformation behaviors. It is found that the comparison results of the experimental and simulated unloading curves are basically coincided and have a quite good agreement. Thus, the comparison of the experimental results and simulated values shows the feasibility and rationality of the simulation methods mentioned above.

    In order to verify the simulation results of the structures, especially for the LD C▬S▬H, nanoindentation simulation is carried out. The simulated elastic modulus is close to the average by experiment, which has demonstrated the reasonableness of simulation. From the simulation results and experimental data [44], it indicates that the selected simulation parameters are reliable, which thus provides the basis and parameters for the multiscale simulation [35, 37].

    In conclusion, the properties of each phase contained in the cement paste are considered in such a way that hydrated and nonhydrated cement phases are assumed to be perfectly elastic. The values of the elastic modulus and of the Poisson coefficient of these different phases are obtained from nanoindentation measurements. The comparison of elastic modulus between references [45, 46] and present work in thesis is shown in Table 5.

    PhasesDensity (g/cm3)Mixing volume (%) [45]Paste volume (%)[45]E (GPa) by referenceE (GPa) by simulationE (GPa) by experiment
    C3S3.1523.401.17117.6 [47]
    C2S3.287.350.78117.6 [47]
    C3A3.034.420.00117.6 [47]
    C4AF3.732.871.39117.6 [47]
    Gypsum2.323.470.0045.7 [48]
    Portlandite (CH)2.240.0013.9646.6 [49]43.13 (AFEM); 45.46 (DFT)≈44.7
    C▬S▬H1.900.0029.0322.4 [50]Monolithic C▬S▬H: C▬S▬H(I) ≈ 75.92; C▬S▬H(II) ≈ 60.95≈22.05
    C▬S▬H pozz2.650.0049.9922.4 [50]By SC scheme: LD C▬S▬H = 18.22; HD C▬S▬H = 31.63≈33.52
    Afm2.020.0015.1242.3 [50]
    Ettringite (AFt)1.780.006.8722.4 [50]
    Porosity total058.4931.690

    Table 5.

    Young’s modulus of phases in Portland cement paste by references and our work.

    The nanoindentation technique has identified in various phases of the cement pates to determine their elastic moduli. The unloading P-h curves under various Poisson’s ratio have been finally obtained, which provides basis for theoretical guidance and practical analysis.

    6. Conclusions

    Nanoindentation simulation is done on ABQUS to get mechanical properties of CH and C▬S▬H, with the comparison of nanoindentation experiment. The effect of Poisson’s ratio on the nanoindentation simulation of CH and C▬S▬H has been quantified through FE simulations. A numerical model was developed to validate and assess the elastic parameters associated with each phase. The P-h curves of the simulation are compared to conclude on the feasibility of the methodology. The experimental elastic moduli of LD C▬S▬H and HD C▬S▬H deviate up to ±20% in our previous work. The present simulation helps to understand more advanced constitutive relations for C▬S▬H (time-dependent deformation, yield-criteria, hardening phenomena, etc.).

    1. The simulated Young’s modulus of CH phase is about 45.459 GPa (DFT methods, by LDA method) and 43.13 GPa (AFEM method) by nanoindentation simulation, which is close to the results of 40.3 ± 4 GPa by Constantinides.

    2. Based on the evaluated elastic modulus values of LD C▬S▬H 18.11 GPa and HD C▬S▬H 31.45 GPa using the self-consistent scheme of the monolithic C▬S▬H(II) structure 60.95 GPa with the size about 5 nm, the simulated results for LD C▬S▬H structure 18.22 GPa and HD C▬S▬H structures 31.63 GPa, which is close to the results of 18.2 ± 4.2 and 29.4 ± 2.4 GPa by Constantinides, thus has demonstrated the reasonableness of simulation parameters.

    3. Finally, on behalf of the purely elastic stage, the unloading curves of the three (CH, LD C▬S▬H, and HD C▬S▬H structures) models are simulated on ABAQUS software. It is found that the comparison results of the experimental and simulated unloading curves are basically coincided and have a quite good agreement. Thus, the analysis of unloading curves shows the feasibility and rationality of the simulation methods mentioned above.

    Overall, as the length scale of C▬S▬H structure with hundreds of nanometers can be called C▬S▬H gel, and these colloids have a number of C▬S▬H particles having an average particle radius about ∼5 nm, the C▬S▬H model composed of C▬S▬H particles with radius of 5 nm for nanoindentation simulation thus can also be precisely described by MD in future, which can be used to compare the results using FEM method, although simulation of these models at larger scale is time-consuming. In the future, the upscaling of the obtained properties at the atomic level will be concerned in order to connect this modeling to continuum models by using relevant input form at small level, carrying forward the critical information to represent the continuum model with the intrinsic nanoscale features incorporated. Other mechanical properties may also be investigated, not only fracture properties and their upscaling toward upper scales, but also creep behavior of the C▬S▬H phase which exhibits such a time-dependent behavior.

    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.

    © 2019 The Author(s). Licensee IntechOpen. This compact is distributed under the terms of the Creative Commons Attribution-NonCommercial 4.0 License, which permits use, distribution and reproduction for non-commercial purposes, provided the original is properly cited.

    How to cite and reference

    Link to this compact Copy to clipboard

    Cite this compact Copy to clipboard

    Jia Fu (May 10th 2019). Atomistic Simulation of Anistropic Crystal Structures at Nanoscale, Atomistic Simulation of Anistropic Crystal Structures at Nanoscale, Jia Fu, IntechOpen, DOI: 10.5772/intechopen.84597. Available from:

    compact statistics

    318total compact downloads

    More statistics for editors and authors

    Login to your personal dashboard for more detailed statistics on your publications.

    Access personal reporting

    Related Content

    This Book

    Atomistic Simulation of Anistropic Crystal Structures at Nanoscale

    Authored by Jia Fu

    Next compact

    Atomistic Simulation of Anistropic Crystal Structures at Nanoscale

    By Jia Fu

    Related Book

    First chapter

    Fourier Series and Fourier Transform with Applications in Nanomaterials Structure

    By Florica Matei and Nicolae Aldea

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

    More About Us