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: 1068

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 Uby Rappe is as [20]:

U=Ur+Uθ+Uϕ+Uω+Uvdw+Ucoul(1)

where Uris bond stretching, Uθis bond angle bending, Uϕis dihedral angle torsion, Uωis improper torsion, Uvdwis nonbonded van der Waals interaction, and Ucoulis 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 Cijklrefers 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 relationsCijnumber 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 εijis 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, Uθ, Uφ, and Uωof Eq. (1) has been investigated in the reference [5], where the molecule's chemical energy is equivalent to the beam strain energy. The coupled field of van der Waals force between atoms or charged ions may be considered by means of the Lennard-Jones potential. Nonlinear spring elements can represent Uvdwpotential.

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 Cijklin general. One can analytically average this tensor, using the full Cijklform 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 Cijcalculated. 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.

Advertisement

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/Lratio, 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/Lratio, 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 Uris for a bond stretch interaction, Uθfor a bond angle bending, Uϕfor a dihedral angle torsion, Uωfor an improper (out-of-plane) torsion, Uvdwfor 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, kθ, and kτare the bond stretching force constant, bond angle bending force constant, and torsional resistance, respectively, and the symbols Δr, Δθ, and Δϕrepresent the bond stretching increment, the bond angle change, and the angle change of bond twisting, respectively.

Assuming that the sections of carbon-carbon bonds are uniformly circular, thus Ix=Iy=I. So parameters of EA, EI, and GJcan 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 Lis the beam length of the C▬C bond, EArepresents tensile stiffness, EIrepresents bending stiffness, GJrepresents torsional stiffness of equivalent beam, kris the stretching constant, kθis the bending constant, and kτis the torsional constant. Obviously, these formulas have established the equivalent relation between the macrobeam structure and microscopic molecular structure. As long as the force constants kr, kθ, and kτare known, the sectional stiffness parameters EA, EI, and GJcan 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, Mifor node iand Fxj, Fyj, Mjfor node j. There are three nodal displacements for each node: ui, vi, θi for node iand 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, kris the stretching stiffness, kyiand kyjare the bending stiffness of joints iand j, and Lis 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 Natoms, 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 ivector and F¯iis the external force exerted on the atom ivector. 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 niterations, xnis the vector of generalized coordinates of atoms after niterations. 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 Kand 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)

Kand Pcan 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 Aand atoms Lamong 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 iis B-J, and the sum of the order number of KAand FAis related to the interaction between atoms Aand 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 Uris for a bond stretch interaction, Uθfor a bond angle bending, Uϕfor dihedral angle torsion, Uωfor an improper (out-of-plane) torsion, Uvdwfor a nonbonded van der Waals interaction, and Ucoulfor Coulombic (electrostatic) interactions. Equivalent beam elements [8] can be used to represent covalent bond to calculate Ur, Uθ, Uϕ, and Uωof Eq. (36), where the molecules’ chemical energy is equivalent to the beam strain energy. The coupled field of van der Waals force and Coulomb force between atoms or charged ions may be obtained by using Lennard-Jones potential or other potentials.

Nonlinear spring elements can represent Uvdwand Ucoulof 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, riis 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 rat 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 qiand qjare partial charges, ε0is 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 Uvdwand Ucoulrepresented 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 r0and F0are 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 rijis the distance between two atoms iand j; εijand σijare parameters of Lennard-Jones. εijrepresents the minimum potential energy between atoms iand j, which corresponds to the most stable interaction (the depth of the potential well under the distance of σij). σijis 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 rbetween 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 rbetween 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 iand jseparated 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 σijand 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 Cijis 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 σ3in z-axis direction is introduced, it will not result in shear strains, thus: γ12 = γ23 = γ31 = 0. We assume σ33expressed by ε33 as: σ33 = E3ε33. Substituting this equation into the ternary equations of stress-strain relation Eq. (42), if ε13, ε23, ε33are 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 E3and elastic constants Cij, the formula of z-axial modulus E3is thus obtained:

E3=c332c132c11+c12(45)

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

E2=c11c121+c12c33c132c11c33c132(46)

Finally, we can see that formulas of E1and E2axial 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 cijcalculated 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 nbecomes ∞ in zigzag nanotube parameters λz1 and λz2expression, 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 Ais the cross-sectional area, A = πDt; Dis the carbon nanotube diameter, and tis 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 △xis 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 r0is 0.3415 nm, and the F0is −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 θ0is 120°, balanced distance constant σc-cis 0.3414 nm, and bond energy constant εc-cis 0.00239 eV. After we obtain the reaction force of y-direction, we substitute all parameters into the formula of ESLGSsand ECNTin 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 nbecomes ∞ (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 Eof 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, kθ, kτ}. It can be seen in Figure 11(b) that the numerical result of E.tof zigzag CNT with various thicknesses is plotted, we can see that E.tis mainly distributed in the range of 0.33–0.37 kJ/m2, and our result of E.tis 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), kθ = 0.876 nN.nm, and kτ = 0.278 nN.nm [8]. The size 3.41 × 3.46 × 0.680 nm of TLGSs is modeled with ABAQUS in Figure 12(b) to determine its elastic constants. van der Waals force of C▬C bond is shown in Figure 12(c). Parameters of Lennard-Jones potential is as [30]: 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 θ0is 120°, balanced distance constant σc-cis 0.3414 nm, and bond energy constant εc-cis 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, Cijis 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-Yplane and the thickness direction of Zis 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 Cijof 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 C33larger 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 inx-direction. (b) Strain loading iny-direction. (c) Strain loading inz-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 Cijshows 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 Cijis 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,zand cijof hexagonal crystal is derived by crystal elastic theory. Besides, a reliable equation to calculate these Cijis 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 E3of hexagonal crystal as a function of elastic constants Cijare 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.tof 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.tis close to be 0.3604 kJ/m2. Besides, the calculated Eof 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 cijare in an acceptable range with experimental values.

  • CH crystal has shown anisotropy. Based on the stress-strain data of AFEM modeling, elastic constant cijis 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.

    Advertisement

    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 nin 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; riand rjare the nucleus position vector; while mε, m, respectively, stand for the quality of nuclei and electrons; ρ(r)is the electron density. Eqs. (3) and (4) include not only the nucleus and electron kinetic energy term, but also interaction term of nucleus-nucleus and electron-electron [8].

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

    ρ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)

    Tis 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 TSis 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 Ccorresponds to iand 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, iand jare 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 aand care 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 Kpoint 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 XBrillouin 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 Kpoints 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 aand bulk modulus B0are 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 Gof 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 cijof CaO at zero temperature and different pressures. It is found that c11varies largely under the effect of pressure, compared with the variations in c12and c44. The change of elastic constant c11represents a produced change with a longitudinal strain in length elasticity. The elastic constants c12and c44are related to the elasticity in shape, which is a shear constant. Besides, c12and c44are less sensitive of pressure as compared with c11. As pressure increases, c11and c12increase monotonically, but c44slightly 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 Gof the MgO agree well with the experimental data [29] and other calculated results [19, 22, 24, 33, 35, 39, 40]. It is noted that c12and c44for 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 c11varies largely under the effect of pressure as compared with the variations in c12and 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: GGAof 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 Cijof 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.

    Advertisement

    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. Yparameter 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 Cand the compliance coefficient Sdo not satisfy the relation CS=I. In other words, the elastic constants of Hill model are not self-consistent. That is why the Yparameter 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 L3direction, 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 L3within 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 L3is as [22, 23]:

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

    Therefore, whole stress σϕψat the direction L3is 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 C44can be calculated as follows:

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

    For cubic crystal, Yelastic 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 σL3in L3direction 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 Yelastic 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