Advances in Materials Science and Engineering
Volume 2017 (2017), Article ID 9529602, 12 pages
https://doi.org/10.1155/2017/9529602

Research Article

A Multiscale Simulation Method and Its Application to Determine the Mechanical Behavior of Heterogeneous Geomaterials

Shengwei Li, Heping Xie, Ru Zhang, Mingzhong Gao, Zetian Zhang, Guo Li, and Jing Xie

State Key Laboratory of Hydraulics and Mountain River Engineering, College of Water Resources and Hydropower, Sichuan University, Chengdu, Sichuan 610065, China

Correspondence should be addressed to Ru Zhang; zhangru@scu.edu.cn

Received 4 March 2017; Accepted 9 July 2017; Published 9 August 2017

Academic Editor: Ying Li

Copyright © 2017 Shengwei Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

To study the micro/mesomechanical behaviors of heterogeneous geomaterials, a multiscale simulation method that combines molecular simulation at the microscale, a mesoscale analysis of polished slices, and finite element numerical simulation is proposed. By processing the mesostructure images obtained from analyzing the polished slices of heterogeneous geomaterials and mapping them onto finite element meshes, a numerical model that more accurately reflects the mesostructures of heterogeneous geomaterials was established by combining the results with the microscale mechanical properties of geomaterials obtained from the molecular simulation. This model was then used to analyze the mechanical behaviors of heterogeneous materials. Because kernstone is a typical heterogeneous material that comprises many types of mineral crystals, it was used for the micro/mesoscale mechanical behavior analysis in this paper using the proposed method. The results suggest that the proposed method can be used to accurately and effectively study the mechanical behaviors of heterogeneous geomaterials at the micro/mesoscales.

1. Introduction

Geomaterials, such as soil, rock, asphalt concrete, and cement concrete, are natural materials or heterogeneous composites composed of natural materials; these materials all involve structural forms at scales that range from an atomic structure to a grain structure and a mineral particle structure and all the way to macroscopic material units. Under the influence of temperature, pressure, water, and other external conditions, the correlation and progression of structural forms at different scales result in the formation of complex geomaterials. Geomaterials failures span these material levels and are therefore a typical multiscale phenomenon. Due to the presence of grain boundaries, lattice defects, intergranular pores, and other microstructures inside rock minerals, microdamage, such as plane cleavage and nucleation and movement of dislocations, forms and interacts at the microstructure locations under an external load; subsequently, this microdamage causes the sliding and fracture of rock mineral grains along grain boundaries with weak bonds. In addition, the convergence and connection of intergranular pores and microcracks can also lead to the accumulation of rock damage and the emergence of rock failures at the mesoscale. The micro/mesostructures inside geomaterials determine the local stress and strain distributions in geomaterials and affect the ultimate failure mode. Kassner et al. [1] believed that only after a multiscale mechanical model has been established can the initiation mechanism of a microfracture be truly understood. The researchers suggested that brittle and ductile materials should be further investigated at the mesoscale. Therefore, the deformation and failure of geomaterials are multiscale evolution processes; that is, it encompasses the formation of induced mesofractures via the accumulation of microdamage in rock mineral crystals and multiscale development.

In actual engineering applications, geomaterials exhibit a multitude of complex characteristics, all of which are related to structural characteristics and mechanical properties at the micro/mesoscales. The microstructures (grain boundaries, lattice defects, and intergranular pores) and mesostructures (size, shape, distribution characteristics, and microcracks of mineral grains) inside rock minerals are the root cause of the various macromechanical properties in a rock body [2]. However, the existing strength and deformation theories have not effectively correlated the atomic scale, mesoscale, lab scale, and engineering scale for geomaterials.

Therefore, by using a specific study method, beginning from the micromechanical behaviors of the mineral grains that form the heterogeneous geomaterials and correlating them with the mesofracture of mineral grains, a multiscale mechanical model for heterogeneous geomaterials can be established. This multiscale model can help to understand the microfracture mechanism of materials and to thus form the foundation for an in-depth analysis of macromechanical properties and fracture behaviors. The results will be of significant importance to the theory development of rock mechanics, the stability of rock engineering, and the adaptability of engineering construction.

With the rapid development of computer technologies in modern materials science, molecular simulations, molecular dynamics (MD) simulations, and other simulation methods for materials at the microscale have gradually become the primary methods used to simulate microscopic systems and to predict macroscopic properties. Since the 1990s, the important achievements made in the studies of geomaterials using molecular simulations were primarily represented by microscale and multiscale research.

Regarding microscale research, Seo et al. [3] reported on nanoscale mechanical properties, including uniaxial compression and shear, for three types of rock minerals using a molecular simulation for the first time. They obtained the corresponding stress-strain response and strength values, and the stiffness matrix elements were in good agreement with the reported experimental results [4]. Sato et al. [5] performed a molecular simulation of the flexibility characteristics for a single lamella of montmorillonite, and an ultimate axial compression of 0.8 GPa was obtained for the veneer bending failure; however, no other mechanical properties were mentioned. Manevitch and Rutledge [6] performed another molecular simulation of the mechanical properties for a single lamella of montmorillonite and obtained the corresponding elastic stiffness matrix; the researchers also noted that determining the lamella thickness is extremely important for determining the elastic modulus of veneers and further obtained the corresponding bending stiffness of veneers according to the thin plate theory. Based on nanocomposite studies, Xu et al. [7] performed a molecular simulation for a single lamella of montmorillonite and obtained an elastic modulus of approximately 277 GPa, which agreed extremely well with the values reported by Manevitch and Chen. Ayoub et al. [8] performed ab initio calculation to examine the structural and mechanical properties of calcite under high pressure conditions.

Regarding multiscale research, there are two classic research methods: the concurrent multiscale method and the sequential multiscale method [9]. For the concurrent multiscale method, a combined model is established by linking the several computational methods together. In the combined model, different scales of material behavior are considered concurrently and communicate using some of handshaking procedure. Currently, the main concurrent multiscale methods [10] include Macroscopic/Atomistic/Ab Initio Dynamics method [11], Quasi Continuum method [12], and coarse-grained molecular dynamics method [13]. For the sequential multiscale method, a series of hierarchical computational methods are linked in some methods, such as homogenization method [14]. In sequential multiscale method, the calculated quantities from a computational simulation at one scale are used to define the parameters of the model operative on the adjoining larger scale. Compared with the concurrent multiscale method, the sequential multiscale method is more convenient. Katti et al. [15, 16] proposed a multiscale framework for interlayer swelling in Na-montmorillonite based on the following: a molecular simulation of the stretching process at the microscale, the fact that lamella can be further divided into discrete elements at the mesoscale and, principally, scanning electron microscopy (SEM) measurements at the macroscale [17]. Following this research idea, Anandarajah and Amarasinghe [18] achieved a mechanistic understanding of the soil-water characteristic curve for unsaturated fine-grained soil at the molecular level based on molecular simulation results. Furthermore, several international scholars have proposed multiscale research that combines molecular dynamics (microscopic), dissipative particle dynamics (DPD) (mesoscopic), and the finite element method (FEM) (macroscopic) [9]. Scocchi et al. [19] reported a multiscale modeling approach for nanocomposites in detail based on the fact that a nonbonded interaction can bridge MD, the input parameters of the mesoscopic DPD, and the macroscopic FEM. Pereira et al. [20] further systematically elucidated the details of each step of the multiscale modeling approach. Similarly, Toth et al. [21] systematically simulated the overall composite strength of several types of polymer-surfactant systems after they were mixed with montmorillonite. They obtained the macroscopic elastic modulus, thermal conductivity, and permeability of the composite materials and then further optimized the composite formula based on these results.

To date, both domestic and international researchers have performed a large number of studies using molecular simulations of clay minerals, polymers, and natural silicate nanocomposites and have proposed micro/mesoscale research ideas that combine atomic scale molecular simulation with mesoscopic SEM or FEM methods [7, 15–20]. This research idea provides an important reference for the developmental-stage studies of the multiscale mechanical behaviors of geomaterials. The present paper focuses on using the molecular simulation method to investigate heterogeneous materials at the microscale. By combining the relatively mature image processing methods to determine the mesomechanical behaviors of geomaterials [22–24] and integrating the results with FEM, the mesoscopic mechanical properties were investigated to achieve a multiscale simulation of the mechanical behaviors of heterogeneous geomaterials.

2. Multiscale Simulation Method for Heterogeneous Materials

2.1. Introduction to Molecular Simulation

Molecular simulations encompass theoretical methods and computational techniques that are used to model or mimic the microscopic behaviors of molecular movements and are widely used in the fields of computational chemistry, computational biology, and materials science. Molecular simulations primarily include quantum mechanics simulations and classical mechanics simulations and are a highly efficient research method for directly revealing the structure and properties of materials at the atomic and molecular levels. The principal ideas for quantum mechanics simulation involve obtaining the electron wave function by solving the Schrödinger equation to calculate the probability of an electron being present in a certain interval. This approach describes the nonlocalized behaviors of electrons outside the nucleus and can fundamentally and accurately compute the molecular structure and material properties.

The Cambridge Serial Total Energy Package (CASTEP) module in the Materials Studio software (Accelrys) is a quantum mechanics simulation code that uses the density functional theory, which addresses problems related to the wave function in a periodic crystal lattice based on the Bloch theorem. Through a series of analyses and calculations, the material properties, including the surface properties, electronic structures, elastic properties, and optical properties, can be obtained. The structure and mechanical parameters, such as the elastic modulus and Poisson’s ratio of heterogeneous materials at the microscale, can be obtained using the CASTEP module in the Materials Studio software. Therefore, this molecular simulation can be used to study the micromechanical behaviors of heterogeneous materials.

2.2. Image Processing Method

In recent decades, with the development of experimental techniques, methods related to image processing, computed tomography, SEM, atomic force microscopy, and other new techniques have been developed that have made measuring parameters and observing physical processes at the mesoscale feasible [25]. Combining these advanced experimental techniques, Yue et al. [22–24] proposed a digital image-based numerical analysis method for heterogeneous geomaterials. In their studies, a digital image processing algorithm was first used to obtain the actual fine structure of the material, which was then vectorized. Then, using mesh automatic generation technology, the finite element mesh was generated for the actual material structure. The coupling of digital image processing techniques and finite elements is an effective method for correlating the fine structures of geomaterials with the response of macroscopic mechanics.

However, heterogeneous geomaterials are typically discontinuous, heterogeneous, and anisotropic. These materials contain defects at different scales that are of different types (grain boundaries, mineral cleavage, microcracks, pores, etc.); these defects decrease the macroscopic mechanical properties significantly more than they decrease the properties at the micro/mesoscales. In this study, digital image processing-based FEM uses the macroscopic mechanical properties of the materials. Therefore, using this method will result in a certain degree of error in the calculations.

2.3. Multiscale Simulation Method

To address the limitations of the digital image processing-based FEM that is widely applied in the area of heterogeneous geomaterials, this paper introduces a molecular simulation to study the micromechanical properties in rocks. By combining an analysis of polished rock slices at the mesoscale and a finite element numerical simulation, a multiscale simulation method based on the molecular simulation of the mechanical properties of heterogeneous materials is proposed to study the micro/mesomechanical behaviors of heterogeneous materials.

In this multiscale simulation method, the image processing method was used to obtain the mesoscale finite element model of the actual fine structure of heterogeneous geomaterials, and then, micromechanical parameters calculated by molecular simulation were assigned to the different structures of the heterogeneous geomaterials for numerical simulation on mesoscale.

Geomaterials are heterogeneous materials consisting of a variety of mineral grains. Each mineral grain shows different mechanical properties, while the interior of the mineral grains can be considered as homogeneous material relatively. At present, mechanical parameters of rock are mainly obtained using macromechanical tests, such as rock triaxial stress tests that can obtain the Young’s modulus, Poisson’s ratio, compressive strength, cohesion, and frictional angle, and the Brazil split test that can obtain the rock tensile strength. The macromechanical parameters are the comprehensive performance of many kinds of mineral grains, including the influence of joints and cracks. Besides, the macroscopic experimental method such as loading rate has great influence on it. Thus, the parameters obtained using macromechanical tests are not representative for the mechanical parameters of mineral grains with less defects. Because the interior of the mineral grains can be considered as homogeneous material, the mechanical parameters of minerals can be obtained by microscopic molecular simulation, which will be more suitable for microsimulation. The method would be more convenient, reliable, and stability, due to the acquisition of parameters not affected by the macrotest method (such as the loading rate). The multiscale method mainly improved the reliability of the calculation but did not improve the computing time obviously. At present, the computation time is acceptable. In the future, we will continue to explore the reasonable ways of simplifying the model to improve the computation efficiency.

The detailed study method is as follows:

(1)An X-ray fluorescence (XRF) spectrometer (XRF-1800 CCED, Shimadzu) and an X-ray diffractometer (DMAX-3C, Rigaku) (as shown in Figure 1) were used in the comprehensively analysis of the composition of heterogeneous geomaterials. The lithology identification and analysis were performed under a microscope at 50x magnification for the rock samples that had been tested by XRF and X-ray diffraction (XRD). Additionally, the polished rock slices were analyzed to verify the results from microscope examination and to obtain the mineral composition of the heterogeneous material.(2)At the microscale, by using the CASTEP module in Materials Studio [26], a computer program for quantum mechanics simulations based on the density functional theory [27, 28], the micromechanical properties of the primary minerals in the heterogeneous material were investigated using a molecular simulation method to generate the mechanical parameters, such as the elastic modulus and Poisson’s ratio.(3)At the mesoscale, using the analysis results of the polished rock slices, CAD modeling was performed to obtain the images of the mineral grain distribution inside the heterogeneous material. These models were then imported into the finite element software ANSYS to establish a finite element model of the actual fine structure. Furthermore, using the micromechanical parameters obtained from a molecular simulation, uniaxial compression and tension simulations were performed to analyze the mechanical behavior and failure characteristics of the material at the micro/mesoscales; the technology roadmap for the study is shown in Figure 2.

(a) X-ray fluorescence

(b) X-ray diffractometer

Figure 1: Composition analysis instrument.

Figure 2: The technology roadmap for the study.

3. The Micro/Mesostructures of Heterogeneous Materials

Rock material is a typical heterogeneous medium comprising aggregates of mineral crystals, cementitious noncrystalline particles, and various defects. Due to the different mineralization mechanisms, its mechanical properties exhibit extremely large spatial heterogeneity; that is, discontinuities exist in the spatial distribution of a rock’s mechanical properties (e.g., elastic modulus and strength). Therefore, kernstone, a discontinuous material, was studied in this paper to determine its micro/mesofracture behaviors using the multiscale simulation method proposed here.

3.1. Mesostructure

3.1.1. Composition Analysis

The XRF results of kernstone showed that the sample contained primarily SiO₂, Al₂O₃, K₂O, and Fe₂O₃; these four oxides accounted for more than 98% of all the oxides in the tested rock core. In addition, the XRD results showed that the kernstone sample was primarily composed of quartz and flint debris along with small amounts of calcarenaceous materials and boulder clay.

3.1.2. Mesostructure of Rock Based on the Analysis of Polished Slices

To verify the test results, the rock sample that had been tested using XRF and XRD was analyzed under a microscope at 50x magnification to identify its lithology, and an analysis of the polished slices was also performed. The identification of the sample composition and the analysis results of the mesostructure under a microscope are shown in Table 1 and Figure 3.

Table 1: Analysis results of the mesostructure.

Sample
Color
Structure
The identification under microscope

Kernstone
Grayish-white
Boulder-containing coarse-grain structures
The mineral compositions included boulder clay, single-crystal quartz, metaquartzite, and flint

Figure 3: Analysis results of the polished slices.

The identification of the kernstone sample under a microscope revealed that the slice was grayish-white in color and that transparent vitreous granulates accounted for 35–40% of the entire sample. Boulder-containing coarse-grain structures were poorly sorted and moderately to excellently rounded with no obvious stratification; silicon-containing particles were reflective and bright, while the others were dim. The mineral compositions included boulder clay, single-crystal quartz, metaquartzite, and flint, among which metaquartzite and flint were the primary components; calcium metasomatism was part of the metaquartzite. The single-crystal quartz particles, flint particles, and metaquartzite particles all exhibited calcification; the severe calcium metasomatism mostly occurred inside the flint debris, and the evenly distributed particles in sandstone exhibited excellent roundness and a variety of sizes. The contact between the crystals of the metaquartzite grains exhibited either a dissemination structure or a straight edge.

3.2. Microstructure

The analysis of the kernstone mesostructure showed that its primary components included quartz, calcite, kaolinite, and feldspar. To simplify the analysis, quartz and calcite were the only components used. Quartz, the chemical composition of which is SiO₂, belongs to the trigonal crystal system and the 154P3221 space group; its lattice parameters are as follows: a = 4.901 Å, b = 4.91 Å, c = 5.40 Å, interfacial angle = = 90°, and = 120°. The quartz unit cell is shown in Figure 4(a) (where the Si atoms are yellow and the O atoms are red).

(a) Quartz

(b) Calcite

Figure 4: Unit cell diagram of the mineral crystals.

Calcite, the chemical composition of which is CaCO₃, belongs to the trigonal crystal system and the 167R3C space group; its lattice parameters are as follows: a = 4.99 Å, b = 4.99 Å, c = 17.06 Å, interfacial angle = = 90°, and = 120°. The calcite unit cell is shown in Figure 4(b) (where the Ca atoms are green, the C atoms are gray, and the O atoms are red).

Figure 4 and Table 2 show a large difference between the microstructures of the two different rock mineral grains in kernstone. The typical structural unit for quartz is a silicon-oxygen tetrahedron comprising Si and O: each O sits at the vertices of the tetrahedron, the Si sits in the center of the tetrahedron, and the tetrahedra are connected at the vertices to form a network, where the Si and O are connected through the Si-O bond. Because the bond energy of a Si-O bond is 460 kJ/mol, the resultant crystal is an atomic crystal. In the crystal structure of calcite, the C-O bond is dominated by a covalent bond, and the bond energy of a C-O bond is 326 kJ/mol. In contrast, the ionic bond dominates in the Ca-O bond. Because the two bonds (covalent versus ionic) in the structure are markedly separated, the resultant structure is a crystal that contains various bond types. The physical properties are mostly dependent on the nature of the bond between O²⁻ and Ca²⁺. Therefore, regarding the classification of crystal types, calcite should be designated as an ionic crystal. Mineral fracture primarily occurs due to the breaking of the Ca-O bond; the C-O bond is difficult to break.

Table 2: Comparison of two different mineral crystals.

Crystal
Chemical composition
Crystal type
Crystal systems
Interparticle interaction
Density/g⋅cm⁻³

Quartz
SiO₂
Atomic crystal
Trigonal
Covalent bond
2.65

Calcite
CaCO₃
Ionic crystal
Trigonal
Ionic bond, covalent bond
2.71

4. Molecular Simulation of the Micromechanical Properties of Heterogeneous Geomaterials

4.1. Molecular Simulation Process

The molecular simulation of the micromechanical properties of rock mineral crystals using CASTEP included modeling, structural optimization, and property calculation.

The data of the α-quartz crystal and calcite crystal were selected from the database in the Material Studio software package to construct unit cell models of quartz and calcite. The quartz unit cell belongs to the trigonal crystal system and the 154P3221 space group; its lattice parameters are as follows: a = 4.901 Å, b = 4.91 Å, c = 5.40 Å, interfacial angle = = 90°, and = 120°. Calcite belongs to the trigonal crystal system and the 167R3C space group; its lattice parameters are as follows: a = 4.99 Å, b = 4.99 Å, c = 17.06 Å, interfacial angle = = 90°, and = 120°.

Then, structural optimization was performed for the α-SiO₂ crystal and calcite crystal. In the calculation, the Perdew-Burke-Ernzerhof (PBE) version of the generalized gradient approximation (GGA) was used to treat the exchange-correlation energy between electrons [29], while the ultrasoft pseudopotential plane-wave method was used to describe the interaction between the ion core and valence electrons. To improve the computation speed as well as to ensure sufficient computational accuracy, a computed cut-off energy of the plane wave of 380 eV was chosen. The Broyden-Fletcher-Goldfarb-Shanno (BFGS) method was used for model optimization [30], where the convergence precision for the iterative process was 5.0 × 10⁻⁶ eV/atom; the convergence criterion for the interaction force between atoms was set as 0.1 eV/nm; the maximum atomic displacement should be less than 5.0 × 10⁻³ nm; and the convergence criterion for internal stress was set as 0.02 GPa. The Monkhorst-Pack method [31] was used to integrate the special points in the Brillouin zone. For α-quartz and CaCO₃, the -mesh points were chosen at 665 and 662, respectively, to ensure the convergence of the energy system and conformation at the basal level of the quasi-complete plane wave.

Finally, the calculation was performed using the abovementioned parameters, and the geometry-optimized models for α-quartz and calcite crystals were obtained. The elastic properties of the optimized models were then calculated using the Elastic Anisotropy Measures code [32].

4.2. Micromechanical Properties

The elastic stiffness matrix and flexibility matrix of the α-quartz and calcite crystals can be obtained from the elastic property calculation, as shown in Tables 3–6. The values of elastic constants are in accordance with earlier experimental works and theoretical values [8].

Table 3: Elastic stiffness matrix of quartz crystal.

(GPa)
1
2
3
4
5
6

1
186.93967
74.68129
103.70831
−0.31848
0
0

2
74.68129
186.93967
103.70831
0.31848
0
0

3
103.70831
103.70831
177.97202
0
0
0

4
−0.31848
0.31848
0
88.98979
0
0

5
0
0
0
0
88.98979
−0.31848

6
0
0
0
0
−0.31848
56.12919

Table 4: Flexible stiffness matrix of quartz crystal.

(GPa)
1
2
3
4
5
6

1
0.008006
−0.0009
−0.00414

0
0

2
−0.0009
0.008006
−0.00414

0
0

3
−0.00414
−0.00414
0.010444
0
0
0

4

0
0.011238
0
0

5
0
0
0
0
0.011238

6
0
0
0
0

0.017816

Table 5: Elastic stiffness matrix of calcite crystal.

(GPa)
1
2
3
4
5
6

1
161.69417
68.80452
63.811
19.68031
0
0

2
68.80452
161.69417
63.811
−19.68031
0
0

3
63.811
63.811
92.89464
0
0
0

4
19.68031
−19.68031
0
35.24155
0
0

5
0
0
0
0
35.24155
19.68031

6
0
0
0
0
19.68031
46.44482

Table 6: Flexible stiffness matrix of calcite crystal.

(GPa)
1
2
3
4
5
6

1
0.0105519
−0.0035507
−0.0048092
−0.0078754
0
0

2
−0.0035507
0.0105519
−0.0048092
0.0078754
0
0

3
−0.0048092
−0.0048092
0.017372
0
0
0

4
−0.0078754
0.0078754
0
0.0371715
0
0

5
0
0
0
0
0.0371715
−0.0157509

6
0
0
0
0
−0.0157509
0.0282051

Because the strain energy density function, , must be greater than zero (positive definiteness), if the structure is deemed to be stable, the stability criterion in (1) must be met:

Further analysis of the elastic constants of the quartz and calcite crystals (as shown in Tables 3–6) shows that all the elastic constants of the two crystals meet the mechanical stability criterion, indicating that the calculated structures of the two crystals are stable.

Based on the Voigt-Reuss-Hill (VRH) approximation (the approximation solutions of Voigt, Reuss, and Hill are the maximum, minimum, and average of the modulus, resp.), the volume elastic modulus and shear modulus of the two mineral crystals can be obtained, as shown in Table 7.

Table 7: Bulk modulus and shear modulus of the two mineral crystals.

Crystal
Mechanical parameters
Voigt
Reuss
Hill

Quartz
The volume elastic modulus /Gpa
124
123.56
123.78

The shear modulus /Gpa
64.81
56.94
60.87

Calcite
The volume elastic modulus /Gpa
89.90
82.39
86.15

The shear modulus /Gpa
38.04
29.17
33.61

The basic micromechanical parameters, such as the elastic modulus and Poisson’s ratio , can be obtained by substituting the calculation results of the volume modulus and shear modulus into (2); the results are shown in Table 8.

Table 8: Elastic modulus and Poisson’s ratio of the two mineral crystals.

Crystal
/GPa

Quartz
156.9
0.289

Calcite
89.23
0.327

Tables 7 and 8 show that the elastic modulus, volume modulus, and shear modulus of the calcite crystal were much smaller than those of the quartz crystal; however, Poisson’s ratio of quartz was smaller than that of calcite. The elastic modulus, volume modulus, and shear modulus are parameters that measure a material’s capability to resist elastic deformation, whereas Poisson’s ratio reflects the elastic constant of transverse deformation. From the microscale viewpoint, these micromechanical parameters comprehensively reflect the bond strength between atoms, ions, and molecules. A quartz crystal is an atomic crystal, while a calcite crystal is an ionic crystal; the covalent bond energy is greater than the ionic bond energy. As a result, the elastic modulus, volume modulus, and shear modulus of a quartz crystal are greater than those of a calcite crystal.

5. Mesoscale Simulation of Heterogeneous Geomaterials Based on the Micromechanical Parameters Obtained from the Molecular Simulation

5.1. Mesoscale Simulation of Heterogeneous Geomaterials

Figure 3 shows that the kernstone sample studied in this paper was an extremely complex, discontinuous medium. The minerals inside primarily included boulder clay, single-crystal quartz, metaquartzite, and flint, of which the metaquartzite and flint were the primary components. Part of the metaquartzite exhibited calcium metasomatism. The components of the kernstone exhibited vastly different mechanical properties under an external load; complex interactions occurred simultaneously among the components. As a result, the mechanical properties of the rock materials (e.g., stress transfer, failure mode, crack propagation, and loading capacity) significantly differed from those of homogeneous rock materials. Therefore, the rock mesostructure model was constructed based on the analysis results of the polished slices of kernstone.

Because the boundary of each composition was clear and easy to distinguish, AutoCAD was used to trace the boundary lines for each composition, as shown in Figure 3; these boundary lines are shown in Figure 5. For convenience of analysis and calculation, the compositionally close components, such as quartz, quartz single crystal, flint, and metaquartzite, were all simplified to quartz; similarly, calcareous minerals and calcium metasomatism were all simplified to calcite.

Figure 5: The boundary lines for each composition in the kernstone.

The models for the polished slices of kernstone in Figure 5 were imported to ANSYS, a software for finite element analysis. Using the automatic mesh generation function in ANSYS, meshes of different materials were generated; the mesoscale finite element model of kernstone was obtained, as shown in Figure 6.

Figure 6: Mesoscale finite element model of the kernstone.

5.2. Parameter Selection and Boundary Condition Setup

5.2.1. Parameter Selection

The elastic constitutive model was used in the numerical simulation. The material parameters were taken from the micromechanical parameters, such as the elastic modulus and Poisson’s ratio, of the two crystals (i.e., quartz and calcite) obtained from the molecular simulation (Table 9). Then, these values were assigned to the corresponding mesoscale unit cells in the finite element model.

Table 9: Mechanical parameters of the materials.

Materials
/GPa

Flint, metaquartzite
156.9
0.289

Calcareous
89.23
0.327

5.2.2. Boundary Condition Setup

The simulation of the mesofracture behaviors primarily included the numerical simulation of uniaxial compression and tension. The specific boundary conditions and load applications are shown in Figure 6: the transverse and vertical degrees of freedom for the node at the left lower corner of the model (point A) were all restrained, the bottom used the symmetry boundary condition (Symmetry B.C.), and displacement loading was applied to the top of the model. To simulate uniaxial compression, the vertical displacement applied to the top was referenced to the peak strain obtained from the uniaxial compression test in the lab, which was 1. Similarly, to simulate uniaxial tension, the vertical displacement applied to the top was referenced to the peak strain obtained from the indirect tensile test in the lab, which was 1%.

5.3. Results and Discussion

By using the above boundary conditions and loading methods, the uniaxial compression and tension simulations for the kernstone at the micro/mesoscales were performed; the distribution nephograms of the first principal stress inside the rock were obtained, as shown in Figure 7. This figure shows that, in uniaxial compression, the tensile stresses for most regions were 0; large tensile stress occurred only close to the boundary of quartz and calcite, and the maximum tensile stress was 30.2 MPa. In the uniaxial tension simulation, tensile stress appeared throughout the entire simulation region and exhibited a more rapid change near the mineral boundary lines; the maximum tensile stress was 196.0 MPa. The results indicate that, because of the differences in the properties of the microscale unit cells, the weak structure produced in mineralization easily led to a stress concentration at the material interface.

(a) Uniaxial compression

(b) Uniaxial tension

Figure 7: The distribution nephograms of the first principal stress inside the kernstone under uniaxial simulations.

To better analyze the micro/mesomechanical behaviors of the rocks, the node stresses were extracted at the material boundary. The vector graphs of the first principal stress distribution at the boundary in the uniaxial compression and tension simulations of kernstone were obtained, as shown in Figure 8.

(a)

(b)

Figure 8: The distribution nephograms of the first principal stress in the mineral boundary of the kernstone under uniaxial simulations: (a) uniaxial compression and (b) uniaxial tension.

Under compression and tension, tensile stress was generated at the crystal boundaries of the kernstone. The tensile strength of the rock is much less than the compressive strength. A preliminary conclusion is drawn that, under continuous loading, the differences of the microstructure and micromechanical properties of different mineral crystals probably lead to a tensile stress concentration and inhomogeneous deformation at the mineral grain boundaries, as well as the failure of the mineral grain interior in the rock, which will cause the occurrence of mesocracks of mineral grains along the boundaries. The focus of the paper is to propose a multiscale simulation framework to study the micro/mesofracture behaviors of heterogeneous geomaterials. Thus, many factors were simplified in the model. Further research especially the deformation modes of interface will be analyzed more deeply to provide a more authentic result.

6. Conclusions

In this paper, to study the meso/micromechanics of heterogeneous geomaterials, a multiscale simulation method that combines molecular simulation at the microscale, a mesoscale analysis of polished slices, and finite element numerical simulation is proposed. This method was preliminarily explored by using it to study the micro/mesomechanical behaviors of kernstone. First, XRF and XRD were used to characterize the composition of a kernstone sample, which primarily included quartz and flint debris along with small amounts of calcarenaceous materials and boulder clay. Then, the primary minerals, that is, quartz and calcite, were chosen for molecular simulation at the microscale to obtain their micromechanical parameters. These micromechanical properties, such as the elastic modulus and Poisson’s ratio, of rock mineral crystals, including quartz and calcite, are known to differ significantly. Finally, after processing the images of the polished rock slices at the mesoscale, these images were imported into finite element analysis software to construct a finite element model. By assigning the micromechanical parameters of the mineral crystals obtained from the molecular simulation to the different minerals in the mesostructure, rock mechanical experiments were performed at the mesoscale. Under continuous loading, the fracture of mineral crystals at the microscale led to stress concentrations and inhomogeneous deformation at the mineral grain boundaries and to the failure of mesostructures inside the rock at the mesoscale, thereby resulting in mesocracks in the mineral grains along the boundaries.

This method is a two-dimensional numerical analysis method that considers the actual mesostructures of geomaterials and is therefore applicable to heterogeneous geomaterials. This method compensates for the limitations of traditional numerical methods that employ statistical theory to characterize the heterogeneity of materials. In the proposed method, a mechanical analysis of heterogeneous materials is obtained that more accurately reflects the mechanical properties of heterogeneous materials. Furthermore, the mechanical parameters of the geomaterials used in the numerical simulations were obtained from macroscopic experiments. However, heterogeneous geomaterials typically contain discontinuous, heterogeneous, and anisotropic bodies; these materials contain defects at different scales that are of different types (grain boundaries, mineral cleavage, microcracks, pores, etc.); these defects decrease the macroscopic mechanical properties significantly more than decreasing the properties at the micro- and mesoscales. Unlike the digital image processing-based FEM that is widely applied to heterogeneous geomaterials, the proposed method uses the micromechanical parameters obtained from the molecular simulation for the material parameters; therefore, the simulation results are more realistic. The application of this method in the study of the micro/mesofracture behaviors of kernstone indicates that this multiscale simulation method is effective for studying the mechanical behaviors of heterogeneous geomaterials.

Considering the complexity of the mechanical behaviors of heterogeneous geomaterials, further improvement is required. For example, in the study of the micro/mesofracture behaviors of kernstone in this paper, two-dimensional finite element models were constructed using only ANSYS software. In addition, only the possible fractures at the mineral grain boundaries were considered, and only a qualitative explanation was provided. In future studies, a mesoscale model needs to be established that is closer to actual heterogeneous materials, and suitable boundaries and boundary contact conditions must be selected to perform more in-depth studies of the multiscale mechanical behaviors of heterogeneous materials. Besides, further studies will be carried out to establish the link between mesoscale and macroscale by means of acoustic emission and computerized tomography (CT) technique.

Abbreviations

SEM: | Scanning electron microscope
MD: | Molecular dynamics
DPD: | Dissipative particle dynamics
FEM: | Finite element method
CASTEP: | Cambridge serial total energy package
XRF: | X-ray fluorescence
XRD: | X-ray diffraction
PBE: | Perdew-Burke-Ernzerhof
GGA: | Generalized gradient approximation
BFGS: | Broyden-Fletcher-Goldfarb-Shanno
VRH: | Voigt-Reuss-Hill
CT: | Computerized tomography.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to express sincere gratitude to the Analytical and Testing Center of Sichuan University for their software support of the Cambridge Serial Total Energy Package (CASTEP) module in the Materials Studio software. This paper was financially supported by the State Key Research Development Program (2016YFC0600702).

References

M. E. Kassner, S. Nemat-Nasser, Z. Suo et al., “New directions in mechanics,” Mechanics of Materials, vol. 37, no. 2-3, pp. 231–259, 2005.

Z. Bazant and E. P. Chen, “Scaling of Structural Failure,” Applied Mechanics Reviews, vol. 50, pp. 593–627, 1999.

Y.-S. Seo, Y. Ichikawa, and K. Kawamura, “Stress-strain response of rock-forming minerals by molecular dynamics simulation,” Materials Science Research, vol. 48, pp. 13–20, 1999.

T. J. Ahrens, Mineral Physics and Crystallography: A Handbook of Physical Constants, American Geophysical Union, 1995.

H. Sato, A. Yamagishi, and K. Kawamura, “Molecular simulation for flexibility of a single clay layer,” Journal of Physical Chemistry B, vol. 105, no. 33, pp. 7990–7997, 2001.

O. L. Manevitch and G. C. Rutledge, “Elastic properties of a single lamella of montmorillonite by molecular dynamics simulation,” The Journal of Physical Chemistry B, vol. 108, pp. 1428–1435, 2004.

W. Xu, Q. Zeng, and A. Yu, “Young's modulus of effective clay clusters in polymer nanocomposites,” Polymer (United Kingdom), vol. 53, no. 17, pp. 3735–3740, 2012.

A. Ayoub, A. Zaoui, and A. Berghout, “High-pressure structural phase transitions and mechanical properties of calcite rock,” Computational Materials Science, vol. 50, no. 3, pp. 852–857, 2011.

Q. H. Zeng, A. B. Yu, and G. Q. Lu, “Multiscale modeling and simulation of polymer nanocomposites,” Progress in Polymer Science, vol. 33, no. 2, pp. 191–269, 2008.

Y. Sheng and X.-G. Zeng, “The deformation mechanisms in process of crack propagation for alpha titanium with compounding microdefects,” Advances in Materials Science and Engineering, vol. 2016, Article ID 2156936, 2016.

F. F. Abraham, J. Q. Broughton, N. Bernstein, and E. Kaxiras, “Spanning the continuum to quantum length scales in a dynamic simulation of brittle fracture,” Europhysics Letters, vol. 44, no. 6, pp. 783–787, 1998.

E. B. Tadmor, M. Ortiz, and R. Phillips, “Quasicontinuum analysis of defects in solids,” Philosophical Magazine A, vol. 73, no. 6, pp. 1529–1563, 1996.

R. E. Rudd and J. Q. Broughton, “Coarse-grained molecular dynamics and the atomic limit of finite elements,” Physical Review B, vol. 58, article R5893, 1998.

J. M. Guedes and N. Kikuchi, “Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods,” Computer Methods in Applied Mechanics and Engineering, vol. 83, no. 2, pp. 143–198, 1990.

D. R. Katti, M. I. Matar, K. S. Katti, and P. M. Amarasinghe, “Multiscale modeling of swelling clays: A computational and experimental approach,” KSCE Journal of Civil Engineering, vol. 13, no. 4, pp. 243–255, 2009.

D. R. Katti and K. S. Katti, “Molecular interactions influence barrier and mechanical properties in swelling clays: A multiscale modeling and experimental investigation,” in Proceedings of the Geo-Frontiers 2011: Advances in Geotechnical Engineering, pp. 4263–4272, Dallas, Tex, USA, March 2011.

D. R. Katti and V. Shanmugasundaram, “Influence of swelling on the microstructure of expansive clays,” Canadian Geotechnical Journal, vol. 38, no. 1, pp. 175–182, 2001.

A. Anandarajah and P. M. Amarasinghe, “Microstructural Investigation of Soil Suction and Hysteresis of Fine-Grained Soils,” Journal of Geotechnical and Geoenvironmental Engineering, vol. 138, no. 1, pp. 38–46, 2011.

G. Scocchi, P. Posocco, A. Danani, S. Pricl, and M. Fermeglia, “To the nanoscale, and beyond!. Multiscale molecular modeling of polymer-clay nanocomposites,” Fluid Phase Equilibria, vol. 261, no. 1-2, pp. 366–374, 2007.

S. P. Pereira, G. Scocchi, R. Toth et al., “Multiscale modelling of polymer/clay nanocomposites,” Journal of Multiscale Modelling, vol. 03, pp. 151–176, 2011.

R. Toth, F. Santese, S. P. Pereira et al., “Size and shape matter! A multiscale molecular simulation approach to polymer nanocomposites,” Journal of Materials Chemistry, vol. 22, no. 12, pp. 5398–5409, 2012.

Z. Q. Yue, S. Chen, and L. G. Tham, “Finite element modeling of geomaterials using digital image processing,” Computers and Geotechnics, vol. 30, pp. 375–397, 2003.

S. Chen, Z. Q. Yue, and L. G. Tham, “Digital image based approach for three-dimensional mechanical analysis of heterogeneous rocks,” Rock Mechanics and Rock Engineering, vol. 40, pp. 145–168, 2007.

W. J. Xu, Z. Q. Yue, and R. L. Hu, “Study on the mesostructure and mesomechanical characteristics of the soil-rock mixture using digital image processing based finite element method,” International Journal of Rock Mechanics and Mining Sciences, vol. 45, pp. 749–762, 2008.

R. de Borst, “Challenges in computational materials science: Multiple scales, multi-physics and evolving discontinuities,” Computational Materials Science, vol. 43, no. 1, pp. 1–15, 2008.

S. J. Clark, M. D. Segall, C. J. Pickard et al., “First principles methods using CASTEP,” Zeitschrift fur Kristallographie, vol. 220, no. 5-6, pp. 567–571, 2005.

P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Physical Review, vol. 136, no. 3B, pp. B864–B871, 1964.

W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Physical Review, vol. 140, no. 4A, pp. A1133–A1138, 1965.

J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple [Phys. Rev. Lett. 77, 3865 (1996)],” Physical Review Letters, vol. 78, no. 7, p. 1396, 1997.

B. G. Pfrommer, M. Côté, S. G. Louie, and M. L. Cohen, “Relaxation of crystals with the quasi-newton method,” Journal of Computational Physics, vol. 131, no. 1, pp. 233–240, 1997.

H. J. Monkhorst and J. D. Pack, “Special points for Brillouin-zone integrations,” Physical Review. B. Solid State, vol. 13, no. 12, pp. 5188–5192, 1976.

A. Marmier, Z. A. D. Lethbridge, R. I. Walton, C. W. Smith, S. C. Parker, and K. E. Evans, “ElAM: A computer program for the analysis and representation of anisotropic elastic properties,” Computer Physics Communications, vol. 181, no. 12, pp. 2102–2115, 2010.