Multiscale mechanical modeling of hydrated cement paste under tensile load using the combined DEM-MD method

Abstract

In this paper, a combined DEM-MD method is proposed to simulate the crack failure process of Hydrated Cement Paste (HCP) under a tensile force. A three-dimensional (3D) multiscale mechanical model is established using the combined Discrete Element Method (DEM)-Molecular Dynamics (MD) method in LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator). In the 3D model, HCP consists of discrete particles and atoms. Simulation results show that the combined DEM-MD model is computationally efficient with good accuracy in predicting tensile failures of HCP.

Keywords

hydrated cement paste / multiscale / MD simulation / DEM

Introduction

Hydrated cement concrete is widely used in buildings and pavements. Cement concrete is inherently heterogeneous and consisting of different phases. In service life, most mechanical failures of hydrated cement concrete originate from the mechanical failures of hydrated cement paste. At macroscale, mechanical failures are relatively well understood with numerous models using finite element method (FEM), discrete element method (DEM), and so on. DEM was first proposed by Cundall and Strack [1] for simulating movement and interaction of assemblies of rigid particles. In their original model, the movements and interactions of unbound assemblies of spherical particles subjected to external loadings are simulated. Later the DEM technique is expanded to material mechanical properties evaluation area. However, mechanical failures are usually originated from atomic bond breakages along specific planes under external loadings, where this conclusions are valid for most infrastructure materials [2,3]. The underlying macrostructure/microstructure and interactions between particles and molecules for most civil infrastructure materials [4–6] are yet to be well understood. Consequently, there is an urgent need to simulate the failure of hydrated cement paste at different scales for better mix design.

Hydrated cement paste is the chemical reaction products of cement hydration. Calcium silicate hydrate (C-S-H) is formed by chemical reactions between silica and calcium hydroxide interstitial solutions and by chemical reactions between calcium salt and alkaline silicate solutions [7]. Research results show that the volume fractions of C-S-H and calcium hydrate (CH) fluctuate around 0.45 and 0.15, respectively, regardless the water-to-cement (w/c) ratio; conversely, the saturated porosity increases with the increase of w/c ratio. The C-S-H phase contributes to 26%–45% of bulk modulus and 28% to 50% of shear modulus, and the fractional contribution to both bulk and shear moduli increases with the increase of w/c ratio from 0.25 to 0.6 [8]. How CH and C-S-H phases contribute to mechanical properties of cement paste at macroscale is still not clearly understood yet. Consequently, it is essential to investigate mechanical properties of C-S-H at fine scales to better understand mechanical performances of cement concrete at macroscale.

Mechanical behaviors of C-S-H haven been intensively investigated through both laboratory experiments and numerical simulations. Pellenq et al. [9]proposed a molecular dynamics model of C-S-H with an average chemical composition of (CaO)_1.65(SiO₂)(H₂O) _1.75, predicting essential physical and mechanical properties at nanoscale, with validations from nanoindentation experimental results. Based on actually chemical compositions of C-S-H, Hou et al. [10] established a molecular dynamics (MD) model of C-S-H with a central void to investigate the tensile strength and tensile failure at nanoscale; stress concentration around the central void and bending sheet of calcium silicate are observed, leading to complicated tensile failure of heterogeneous layered C-S-H gel.

Even though there are various models of C-S-H at nanoscale available [7,8,11–14], the fundamental microstructure and mechanical properties of C-S-H is still not well understood yet. However, there are some common features discovered [15]. The fundamental units of globules are tobermorite- or jennite-like structures; as globules pack together, C-S-H gels with different packing density form into layer-based structures [14]. Depending on the packing density, there are two types of C-S-H: high density calcite-silicate-hydrate (HD C-S-H) and low density calcite-silicate-hydrate (LD C-S-H), shown in Fig. 1 [13,16,17]. HD C-S-H has a denser structure with a greater density and less porosity, compared to LD C-S-H.

Research shows that the C-S-H phase of cement paste consists of colloidal particle [18], and C-S-H has cohesive characteristics at nanoscale due to attractive surface forces [15,19]. The corresponding colloidal models are proposed to analyze the colloidal nature of C-S-H gel. Diamond [20] proposed a colloidal model named Munich model to analyze the mechanical behavior of cement hydration products; in this model, discrete spheres represent C-S-H gel with disjoining pressure between particles and spaces between spheres represent capillary pores outside of the C-S-H gel. Jennings [13] proposed a model to consider colloidal particles – globules without long-range layered structure for the C-S-H in cement paste, in order to accurately predict physical properties such as specific area and density. Thomas and Jennings [14] proposed a colloidal model to analyze chemical aging of C-S-H with respect to surface area, drying shrinkage, creep, and permeability, and this model addresses the change of sol-gel transition in the reversible and irreversible condensation reactions due to creep and shrinkage. The aforementioned colloidal models only investigate the material structure and mechanical properties of C-S-H at nanoscale, without taking the heterogeneous nature of C-S-H at multiple scales into considerations.

Multiscale mechanical modeling techniques are developed for concrete materials, such as multiscale lattice model, multiscale finite element model, combined FEM-molecular dynamics (MD) method [6,16,21,22]. There have also been significant researches in the multiscale simulation areas [23–26]. The aforementioned multiscale mechanical models consider concrete material as either homogeneous or heterogeneous continuum and evaluate mechanical properties at macro- and micro- scales. However, such models cannot represent the discrete nature of colloidal C-S-H gel at nanoscale and minerals at atomic-scale.

In the multiscale simulation area, it mainly includes three kinds: hierarchical method [27], concurrent method [28] and semi-concurrent method [29]. Among all the approaches, there have been significant progresses in the multiscale methods for mechanical properties including fracture related properties [23–26]. In order to properly address the colloidal characteristics of C-S-H, this paper proposes an innovative multiscale mechanical modeling method, named the combined discrete element model-molecular dynamics (DEM-MD) model approach. In this DEM-MD model, discrete particles and atoms represent C-S-H gel and capillary pores, respectively. The DEM-MD model can establish a colloidal model for C-S-H gel, with the spatial scale ranging from atomic-scale to microscale. The failure of C-S-H gel under tensile load is analyzed using the DEM-MD model to investigate the mechanical properties of C-S-H gel under tensile loading at both atomic-scale and microscale. The computational scheme of the DEM-MD approach is introduced, followed by simulations. Final research findings are summarized at the end of this paper.

The combined discrete element model-molecular dynamics model approach

In the DEM-MD model, interactions between discrete particles and atoms are described by contact models. There are three contact models: LARGE_LARGE (R₁, R₂>0), SMALL_LARGE (R₁>0, R₂=0), and SMALL_SMALL (R₁=R₂=0). LARGE_LARGE describes the particle-particle interaction, and that is the interactions between C-S-H colloidal gel. SMALL_LARGE contact describes the atom-particle interaction, i.e., the interactions between C-S-H colloidal gel and gas molecules in the capillary pore. SMALL_SMALL contact describes the atom-atom interactions, i.e., the interaction between gas molecules in the capillary pore.

The computation is performed by revising the pair_colloid package in LAMMPS. In the original pair_colloid package in LAMMPS, there are two type of errors that may come up in calculation.

(i) The original pair_colloid program will stop running and present an error message, while r≤R₁+R₂ for LARGE_LARGE case: “Overlapping large/large in pair colloid”.

(ii) The original pair_colloid program will stop running and present an error message, while r≤R₁+R₂ for SMALL_LARGE case: “Overlapping small/large in pair colloid”.

Figure 2 shows the interactions between atoms and particles for three contact models. To avoid the aforementioned two errors in computations for the atom-particle interaction and the particle-particle interaction, the following revisions were made in computation algorithm. The van der Waals force between overlapping particles is considered as the summation of van der Waals force between two separated particles and van der Waals force between two parallel surfaces (shown in Fig. 2(c)). Instead of considering both adhesive and repulsive components of van der Waals interaction, van der Waals force between overlapping particles is now considered as the summation of the repulsive force between two displaced particles with Hooke contact model and the adhesive component of van der Waals force between two overlapping particles (shown in Fig. 2).

Figure 3 shows the calculation flowchart in the DEM-MD model. It is a concurrent multiscale model, updating the interactions among particles and atoms at every time step. The input parameters in the model include particle radius, Hamaker coefficients for the three contact models, separation distanceD, cutoff radiuses in the particle-particle contact. The calculation output include temperature, velocity, stress, strain, and so on. The three contact models are as follows:

a) Atom-atom interaction: Case SMALL_SMALL, Fig. 2(a)

(1)

where U_{\mathrm{total,ss}} = van der Waals potential energy between two atoms, energy unit; \epsilon = the depth of potential well, distance unit; \sigma = the finite distance at which the inter-particle van der Waals potential is zero potential energy constant, distance unit;r = distance between two atoms, distance unit; A_ss = Hamaker constant for these two atoms, related to \epsilon, and \sigma, energy unit.

b) Overlapping atom-particle interaction: CASE SMALL_LARGE, Fig. 1(b)

(2)

(3)r\le R_{\mathrm{1}}+R_{\mathrm{2}},U_{\mathrm{total,cs}}=U_{\mathrm{atom-particle}}+U_{\mathrm{atom-surface}}\times \pi a^{\mathrm{2}},

(4)U_{\mathrm{atom-particle}}=U_{\mathrm{A,atom-particle}}+U_{\mathrm{R,atom-particle}},

(5)U_{\mathrm{A,atom-particle}}={\displaystyle \frac{\pi \lambda q}{D+R_{\mathrm{1}}}}\left[{\displaystyle \frac{D^{\mathrm{2}}+\mathrm{2}D{R}_{\mathrm{1}}}{\mathrm{4}}}\cdot {\displaystyle \frac{D^{\mathrm{4}}-{(D+\mathrm{2}{R}_{\mathrm{1}})}^{\mathrm{4}}}{D^{\mathrm{4}}{(D+\mathrm{2}{R}_{\mathrm{1}})}^{\mathrm{4}}}}+{\displaystyle \frac{\mathrm{2}(D+R_{\mathrm{1}})}{\mathrm{3}}}\cdot {\displaystyle \frac{{(D+\mathrm{2}{R}_{\mathrm{1}})}^{\mathrm{3}}-D^{\mathrm{3}}}{D^{\mathrm{3}}{(D+\mathrm{2}{R}_{\mathrm{1}})}^{\mathrm{3}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{{(D+\mathrm{2}{R}_{\mathrm{1}})}^{\mathrm{2}}-D^{\mathrm{2}}}{D^{\mathrm{2}}{(D+\mathrm{2}{R}_{\mathrm{1}})}^{\mathrm{2}}}}\right],

(6)U_{\mathrm{R,atom-particle}}=-{\displaystyle \frac{\pi {\lambda }_{\mathrm{2}}q}{\mathrm{72}r}}\left[{\displaystyle \frac{r+\mathrm{15}{R}_{\mathrm{1}}}{{(r+R_{\mathrm{1}})}^{\mathrm{9}}}}-{\displaystyle \frac{r+\mathrm{3}{R}_{\mathrm{1}}}{{(r-R_{\mathrm{1}})}^{\mathrm{9}}}}\right],\hspace{0.17em}{\lambda }_{\mathrm{2}}=\epsilon {\sigma }^{\mathrm{12}},

(7)U_{\mathrm{atom-surface}}=U_{\mathrm{A,atom-surface}}+U_{\mathrm{R,atom-surface}}\hspace{0.17em},per\hspace{0.17em}unit\hspace{0.17em}area,

(8)U_{\mathrm{A,atom-surface}}={\displaystyle \frac{\pi {\lambda }_{\mathrm{1}}q}{\mathrm{6}{D}^{\mathrm{3}}}},\hspace{0.17em}{\lambda }_{\mathrm{1}}=\epsilon {\sigma }^{\mathrm{6}},

(9)U_{\mathrm{R,atom-surface}}=-{\displaystyle \frac{\pi {\lambda }_{\mathrm{2}}q}{\mathrm{45}{D}^{\mathrm{9}}}},

(10){\lambda }_{\mathrm{2}}=\epsilon {\sigma }^{\mathrm{12}},

where U_{\mathrm{total,cs}} = van der Waals potential between atom and particle, energy unit; U_{\mathrm{atom-particle}} = van der Waals potential between atom and particle while there is overlapping, energy unit;U_{\mathrm{atom-surface}} = van der Waals potential between atom and surface while there is overlapping, energy unit;U_{\mathrm{A,atom-particle}} = van der Waals potential between atom and particle due to attractive force while there is overlapping, energy unit;U_{\mathrm{R,atom-particle}} = van der Waals potential between atom and particle due to repulsive force while there is overlapping, energy unit;U_{\mathrm{A,atom-surface}} = van der Waals potential between atom and surface due to attractive force while there is overlapping, energy unit;U_{\mathrm{R,atom-surface}} = van der Waals potential between atom and surface due to repulsive force while there is overlapping, energy unit;A_cs = Hamaker constant for the atom and particle, energy unit; q= number density of the particle, distance unit ^‒3; r = distance between the atom and the origin of particle, distance unit.

c) Particle-particle interaction: CASE LARGE_LARGE, Fig. 2(c).

The following equations present the energy between any two discrete particles. When the two particles are away from each other, Eq. (11) represent the energy due to attractive van der Waal interaction, and Eq. (12) represent the repulsive van der Waal interaction. When the two particle overlap with each other in the calculation, the energy equations are represented in Eqs. (14)‒(20). At every time step, the energy is updated for any two particles; force applied on every particle is calculated at every time step; the corresponding acceleration of every particle is determined based on Newton’s second law at every time step; velocity and displacement are then determined from acceleration at every time step.

(11)

(12)U_{\mathrm{A,cc}}=-{\displaystyle \frac{A_{\mathrm{cc}}}{\mathrm{6}}}\left[{\displaystyle \frac{\mathrm{2}{R}_{\mathrm{1}}{R}_{\mathrm{2}}}{r^{\mathrm{2}}-{(R_{\mathrm{1}}+R_{\mathrm{2}})}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{2}{R}_{\mathrm{1}}{R}_{\mathrm{2}}}{r^{\mathrm{2}}-{(R_{\mathrm{1}}-R_{\mathrm{2}})}^{\mathrm{2}}}}+\ln \left({\displaystyle \frac{r^{\mathrm{2}}-{(R_{\mathrm{1}}+R_{\mathrm{2}})}^{\mathrm{2}}}{r^{\mathrm{2}}-{(R_{\mathrm{1}}-R_{\mathrm{2}})}^{\mathrm{2}}}}\right)\right],

(13)U_{\mathrm{R,cc}}=-{\displaystyle \frac{A_{\mathrm{cc}}}{\mathrm{37800}}}{\displaystyle \frac{{\sigma }^{\mathrm{6}}}{r}}\left[\begin{array}{l}{\displaystyle \frac{r^{\mathrm{2}}-\mathrm{7}r(R_{\mathrm{1}}+R_{\mathrm{2}})+\mathrm{6}(R_{\mathrm{1}}^{\mathrm{2}}+\mathrm{7}{R}_{\mathrm{1}}{R}_{\mathrm{2}}+R_{\mathrm{2}}^{\mathrm{2}})}{{(r-R_{\mathrm{1}}-R_{\mathrm{2}})}^{\mathrm{2}}}}+{\displaystyle \frac{r^{\mathrm{2}}-\mathrm{7}r(R_{\mathrm{1}}+R_{\mathrm{2}})+\mathrm{6}(R_{\mathrm{1}}^{\mathrm{2}}+\mathrm{7}{R}_{\mathrm{1}}{R}_{\mathrm{2}}+R_{\mathrm{2}}^{\mathrm{2}})}{{(r+R_{\mathrm{1}}+R_{\mathrm{2}})}^{\mathrm{2}}}}\\ -{\displaystyle \frac{r^{\mathrm{2}}+\mathrm{7}r(R_{\mathrm{1}}-R_{\mathrm{2}})+\mathrm{6}(R_{\mathrm{1}}^{\mathrm{2}}-\mathrm{7}{R}_{\mathrm{1}}{R}_{\mathrm{2}}+R_{\mathrm{2}}^{\mathrm{2}})}{{(r+R_{\mathrm{1}}-R_{\mathrm{2}})}^{\mathrm{2}}}}-{\displaystyle \frac{r^{\mathrm{2}}-\mathrm{7}r(R_{\mathrm{1}}-R_{\mathrm{2}})+\mathrm{6}(R_{\mathrm{1}}^{\mathrm{2}}-\mathrm{7}{R}_{\mathrm{1}}{R}_{\mathrm{2}}+R_{\mathrm{2}}^{\mathrm{2}})}{{(r-R_{\mathrm{1}}+R_{\mathrm{2}})}^{\mathrm{2}}}}\end{array}\right],

(14)

(15)U_{\mathrm{particle-particle}}=U_{\mathrm{A,particle-particle}}+U_{\mathrm{R,particle-particle}},

(16)U_{\mathrm{A,particle-particle}}={\displaystyle \frac{{\pi }^{\mathrm{2}}{\lambda }_{\mathrm{1}}{q}_{\mathrm{1}}{q}_{\mathrm{2}}}{\mathrm{6}D}}\left({\displaystyle \frac{R_{\mathrm{1}}{R}_{\mathrm{2}}}{R_{\mathrm{1}}+R_{\mathrm{2}}}}\right),\hspace{0.17em}{\lambda }_{\mathrm{1}}=\epsilon {\sigma }^{\mathrm{6}},

(17)U_{\mathrm{A,particle-particle}}=-{\displaystyle \frac{{\pi }^{\mathrm{2}}{\lambda }_2{q}_{\mathrm{1}}{q}_{\mathrm{2}}}{72}}{\int }_{r-R_2}^{r+R_2}\frac{[R_2^2-{(r-R)}^2]}{r}\cdot [\frac{R+15R_1}{{(R+R_1)}^9}-\frac{R+3R_1}{{(R-R_1)}^9}]\mathrm{d}R,{\lambda }_2={\epsilon \sigma }^{12},

(18)U_{\mathrm{surface-surface}}=U_{\mathrm{A,surface-surface}}+U_{\mathrm{R,surface-surface}},per\hspace{0.17em}unit\hspace{0.17em}area,

(19)U_{\mathrm{A,surface-surface}}={\displaystyle \frac{\pi {\lambda }_{\mathrm{1}}{q}_{\mathrm{1}}{q}_{\mathrm{2}}}{\mathrm{12}{D}^{\mathrm{2}}}},\hspace{0.17em}{\lambda }_{\mathrm{1}}=\epsilon {\sigma }^{\mathrm{6}},

(20)U_{\mathrm{R,surface-surface}}=-{\displaystyle \frac{\pi {\lambda }_{\mathrm{2}}{q}_{\mathrm{1}}{q}_{\mathrm{2}}}{\mathrm{360}{D}^{\mathrm{8}}}},\hspace{0.17em}{\lambda }_{\mathrm{2}}=\epsilon {\sigma }^{\mathrm{12}},

where U_{\mathrm{total,cc}} = van der Waals potential energy between particle and particle, energy unit; U_{\mathrm{A,cc}} = attractive part of van der Waals potential energy between two particles, energy unit;U_{\mathrm{R,cc}} = repulsive part of van der Waals potential energy between two particles, energy unit; U_{\mathrm{particle-particle}} = van der Waals potential energy between two overlapping particles, energy unit;U_{\mathrm{surface-surface}} = van der Waals potential energy between two cutoff surfaces for two overlapping particles, energy unit;U_{\mathrm{A,particle-particle}} = attractive part of van der Waals potential energy between two overlapping particles, energy unit;U_{\mathrm{R,particle-particle}} = repulsive part of van der Waals potential energy between two overlapping particles, energy unit;U_{\mathrm{A,surface-surface}} = attractive part of van der Waals potential energy between two cutoff surfaces for two overlapping particles, energy unit;U_{\mathrm{R}\mathrm{,surface-surface}} = repulsive part of van der Waals potential energy between two cutoff surfaces for two overlapping particles, energy unit; A_{\mathrm{cc}} = Hamaker constant for the two particles, energy unit;R_{\mathrm{1}}, R_{\mathrm{2}} = radius of the two particles, distance unit; q_{\mathrm{1}}, q_{\mathrm{2}} = number density of the two particles, distance unit^-3; r = distance between the origins of two particles, distance unit; r_{\mathrm{cutoff,cc}} = cutoff radius of the van der Waals interaction between atom and particle, distance unit.

Note that our current research does not include the cohesive zone model [30–32], where it can be input into the current multiscale DEM-MD model by revising the potential function in the “cohesive zone area” defined in the simulation.

Computational model

Based on the packing density, C-S-H gels can be discretized into two categories: high-density (HD) C-S-H and low density (LD) C-S-H. Both laboratory experiments (including atomic force microscope (AFM) and indentation tests) and theoretical modeling are used to analyze mechanical properties of cement hydration products. Table 1 shows a brief summary on physical properties of C-S-H gel [33,34]. To better understand the tensile failure in the C-S-H gel, a multiscale model using combined molecular dynamics method and discrete element method is established for C-S-H gel. In this paper, there are two C-S-H models established, including HD C-S-H and LD C-S-H with the porosity of 24% and 37%, respectively.

The two-dimensional computation domain is set as 300×300 s² (s is a fundamental distance units set for a unitless LJ simulation in LAMMPS; in this paper,s represents a constant for C-S-H gel). Boundary conditions are periodic on each side. Isothermal-isobaric (NPT) ensemble is used to reach a stable state of the computational domain. There are a total of 90000 discrete components, including 68381 particles and 21619 atoms. Temperature is set as 1.0 K. Tensile strain is applied in y-direction at a strain rate of 0.001t^‒1. To investigate the influence of particle radius on the mechanical analysis, calculations are performed on HD C-S-H with the particle radius ranging from 2s to 40s, and the cutoff radius is set as 2.5 times of the particle radius. Fig. 4 presents the relationship between the maximum stress and particle radius. It is found that there is an exponential relationship between particle radius and the maximum stress.

Tensile failure of C-S-H colloidal model

Table 2 presents the input parameters for both HD C-S-H and LD C-S-H. Isothermal-isobaric (NPT) ensemble is used to reach a stable state of the two-dimensional computation domain.

Tensile loading is applied at a loading rate of 0.01/t in y-direction. Figure 5 plots the microstructure of HD C-S-H gel after the tensile test in y-direction, with red dots representing C-S-H gel and blue dots representing air in voids, shown in Visual Molecular Dynamics (VMD). Figure 6 shows the stress-strain relationship of HD C-S-H and LD C-S-H with the tensile loading applied iny-direction. As expected, the tensile stress in HD C-S-H is greater than the tensile stress in LD C-S-H. Around the peak stress, the calculated elastic modulus is 13.1 GPa for HD C-S-H and 60.5 GPa for LD C-S-H, falling into the range from 15 GPa to 190 GPa [34,35], based on measurement results and theoretical modeling of the elastic modulus of C-S-H. The unexpected problem is that the elastic modulus of HD C-S-H is smaller than that of LD C-S-H, indicating that the key parameters need further optimizations, especially the particle radiusR₂ and the inter-particle distance D for overlapping particles.

It should be pointed out that our current research is focused on the mechanical properties especially the tensile performance of HCP. However, the effective thermal conductivity is another important factor that affects HCP mechanical behaviors, where different numerical tools [36–45] can be employed to study the mechanical properties of such civil engineering materials. It is suggested that possible future studies on HCP effective thermal conductivity need to consider the temperature effect where it is constant in current simulation, also the potential functions from Eq. (1) to Eq. (20) need to be revised to address the thermal conductivity at nanoscale. Also note size effects are very important in multiscale simulation, since they present constitutive models with strain softening. However, due to our limited computation power, size effects are not considered in current research. In future studies, this issue should be taken into consideration.

Summary

Hydrated cement paste (HCP) is the hydration products of Portland cement. Uniaxial tensile tests on hydrated cement paste samples show a linear elastic behavior before the peak load is reached. After the peak load, HCP exhibits quasi-brittle failures. There are a lot of simulations to analyze tensile failures of HCP using the following methods: finite element method (FEM), discrete element method (DEM), and molecular dynamics (MD) method, etc. Tensile failures of HCP originate from atomic-bond breakages around inherent defects and inner voids, and the atomic failures further develop into microscales, leading to final failures. Therefore, this paper proposes a combined DEM-MD method to simulate the crack failure process of HCP under a tensile force. A three-dimensional (3D) multiscale mechanical model is established using the combined DEM-MD method in LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator). In the 3D model, HCP consists of discrete particles and atoms. Simulation results show that the combined DEM-MD model is computationally efficient with good accuracy in predicting tensile failures of HCP.

The multiscale mechanical model is established using combined DEM method and MD method. Mechanical performances of cement hydration products are simulated using this multiscale mechanical model. C-S-H is selected as a representation of the major cement hydration products. By considering the porosities of HD C-S-H and LD C-S-H, tensile failure properties of HD C-S-H and LD C-S-H are investigated using this multiscale model. As expected, HD C-S-H with a small porosity has a greater maximum tensile stress than LD C-S-H with a large porosity. Further research can be conducted on the optimization of key parameters, especially the particle radius and the inter-particle distanceD. The proposed DEM-MD modeling technique provide an alternative means to investigate the mechanical properties of C-S-H at both atomic-scale and microscale. Note that our current research is still in the initial stage and needs to be developed in future studies.

References

[1] Cundall P A, Strack  O D L. A discrete numerical mode for granular assemblies. Geotechnique, 1979, 29(1): 47–65

[2] Hou Y, Wang  L, Yue P ,  Sun W. Fracture failure in crack interaction of asphalt binder by using a phase field approach. Materials and Structures, 2015, 48(9): 2997–3008

[3] Hou Y, Wang  L, Pauli T ,  Sun W. Investigation of the asphalt Self-Healing mechanism using a Phase-Field model. J Mater Civil Eng, 2015, 27(3): 4014118

[4] Hou Y, Sun  F, Sun W ,  Guo M, Xing  C, Wu J . Quasi-brittle fracture modeling of preflawed bitumen using a diffuse interface model. Advances in Materials Science and Engineering, 2016, 2016: 8751646

[5] Hou Y, Sun  W, Das P ,  Song X, Wang  L, Ge Z ,  Huang Y . Coupled navier–stokes Phase-Field model to evaluate the microscopic phase separation in asphalt binder under thermal loading. J Mater Civil Eng, 2016, 28(10): 4016100

[6] Hou Y, Wang  L, Yue P ,  Pauli T ,  Sun W. Modeling mode I cracking failure in asphalt binder by using nonconserved Phase-Field model. J Mater Civil Eng, 2014, 26(4): 684–691

[7] Allen A J, Thomas  J J, Jennings  H M. Composition and density of nanoscale calcium-silicate-hydrate in cement. Nature Materials, 2007, 6(4): 311–316

[8] Zhang Q. Creep Properties of Cementitious Materials: Effect of Water and Microstructure, An Approach by Microindentation. Dissertation for PhD degree. Paris Est, 2014

[9] Pellenq R J, Kushima  A, Shahsavari R ,  Van Vliet K J ,  Buehler M J ,  Yip S, Ulm  F J. A realistic molecular model of cement hydrates. Proceedings of the National Academy of Sciences of the United States of America, 2009, 106(38): 16102–16107

[10] Hou D, Zhao  T, Wang P ,  Li Z, Zhang  J. Molecular dynamics study on the mode I fracture of calcium silicate hydrate under tensile loading. Engineering Fracture Mechanics, 2014, 131: 557–569

[11] Fonseca P C, Jennings  H M, Andrade  J E. A nanoscale numerical model of calcium silicate hydrate. Mechanics of Materials, 2011, 43(8): 408–419

[12] Hou D, Li  Z. Large-scale simulation of calcium silicate hydrate by molecular dynamics. Adv Cem Res, 2015, 27(5): 278–288

[13] Jennings H M. A model for the microstructure of calcium silicate hydrate in cement paste. Cement and Concrete Research, 2000, 30(1): 101–116

[14] Thomas J J, Jennings  H M. A colloidal interpretation of chemical aging of the C-S-H gel and its effects on the properties of cement paste. Cement and Concrete Research, 2006, 36(1): 30–38

[15] Hou D, Zhu  Y, Lu Y ,  Li Z. Mechanical properties of calcium silicate hydrate (C–S–H) at nano-scale: a molecular dynamics study. Materials Chemistry and Physics, 2014, 146(3): 503–511

[16] Sun W, Wei  Y, Wang D ,  Wang L. Review of multiscale characterization techniques and multiscale modeling methods for cement concrete: from atomistic to continuum. In: Kringos N, Birgisson B, Frost D, Wang L, eds. Multi-Scale Modeling and Characterization of Infrastructure Materials. Springer Netherlands, 2013

[17] Richardson I G . Tobermorite/jennite- and tobermorite/calcium hydroxide-based models for the structure of C-S-H: applicability to hardened pastes of tricalcium silicate, β-dicalcium silicate, Portland cement, and blends of Portland cement with blast-furnace slag, metakaol. Cement and Concrete Research, 2004, 34(9): 1733–1777

[18] Scherer G W. Structure and properties of gels- I. Theory. Cement and Concrete Research, 1999, 29(8): 1149–1157

[19] Powers T C. Structure and physical properties of hardened Portland cement paste. Journal of the American Ceramic Society, 1958, 41(1): 1–6

[20] Diamond S. Hydraulic cement pastes: their structure and properties. Cement and Concrete Association, Slough, UK, 1976, 1–334

[21] Nguyen VP, Stroeven  M, Sluys LJ . Multiscale failure modeling of concrete: micromechanical modeling, discontinuous homogenization and parallel computations. Computer Methods in Applied Mechanics & Engineering, 2012, s 201–204:139–156

[22] Qian Z. Multiscale Modeling of Fracture Processes in Cementitious Materials. Dissertation for PhD degree. Haveka B V, 2012

[23] Talebi H, Silani  M, Rabczuk T . Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Advances in Engineering Software, 2015, 80: 82–92

[24] Silani M, Talebi  H, Hamouda A M ,  Rabczuk T . Nonlocal damage modelling in clay/epoxy nanocomposites using a multiscale approach. J Comput SCI-Neth, 2016, 15: 18–23

[25] Silani M, Ziaei-Rad  S, Talebi H ,  Rabczuk T . A semi-concurrent multiscale approach for modeling damage in nanocomposites. Theoretical and Applied Fracture Mechanics, 2014, 74: 30–38

[26] Talebi H, Silani  M, Bordas S P A ,  Kerfriden P ,  Rabczuk T . A computational library for multiscale modeling of material failure. Computational Mechanics, 2014, 53(5): 1047–1071

[27] Aarnes J E, Krogstad  S, Lie K A . A hierarchical multiscale method for two-phase flow based upon mixed finite elements and nonuniform coarse grids. Siam Journal on Multiscale Modeling & Simulation, 2006, 5(2): 337–363

[28] Gu Y T, Zhang  L C. A concurrent multiscale method based on the meshfree method and molecular dynamics analysis. Siam Journal on Multiscale Modeling & Simulation, 2006, 5(4): 1128–1155

[29] Zhu H, Wang  Q, Zhuang X . A nonlinear semi-concurrent multiscale method for dynamic fractures. International Journal of Impact Engineering, 2015, 87: 65–82

[30] Jin Z H, Sun  C T. A Yoffe crack/cohesive zone model for a steady state moving crack. Mechanics Research Communications, 2016, 71: 44–47

[31] Wu T, Wriggers  P. Multiscale diffusion–thermal–mechanical cohesive zone model for concrete. Computational Mechanics, 2015, 55(5): 999–1016

[32] Karedla R S, Reddy  J N. Modeling of crack tip high inertia zone in dynamic brittle fracture. Engineering Fracture Mechanics, 2007, 74(13): 2084–2098

[33] Ulm F J, Constantinides  G, Heukamp F H . Is concrete a poromechanics materials?—A multiscale investigation of poroelastic properties. Materials and Structures, 2004, 37(1): 43–58

[34] Mondal P. Nanomechanical properties of cementitious materials. Dissertations & Theses- Gradworks, 2008

[35] Manzano H, Masoero  E, Lopez-Arbeloa I ,  Jennings H M . Mechanical behaviour of ordered and disordered calcium silicate hydrates under shear strain studied by atomic scale simulations. International Conference on Creep, Shrinkage, and Durability Mechanics, 2013, 86–97

[36] Alzina A, Toussaint  E, Béakou A ,  Skoczen B . Multiscale modelling of thermal conductivity in composite materials for cryogenic structures. Composite Structures, 2006, 74(2): 175–185

[37] Mortazavi B, Hassouna  F, Laachachi A ,  Rajabpour A ,  Ahzi S, Chapron  D, Toniazzo V ,  Ruch D. Experimental and multiscale modeling of thermal conductivity and elastic properties of PLA/expanded graphite polymer nanocomposites. Thermochimica Acta, 2013, 552: 106–113

[38] Mortazavi B, Pötschke  M, Cuniberti G . Multiscale modeling of thermal conductivity of polycrystalline graphene sheets. Nanoscale, 2014, 6(6): 3344–3352

[39] Tonks M R, Liu  X Y, Andersson  D, Perez D ,  Chernatynskiy A ,  Pastore G ,  Stanek C R ,  Williamson R . Development of a multiscale thermal conductivity model for fission gas in UO 2. Journal of Nuclear Materials, 2016, 469: 89–98

[40] Tonks M R, Millett  P C, Nerikar  P, Du S ,  Andersson D ,  Stanek C R ,  Gaston D ,  Andrs D ,  Williamson R . Multiscale development of a fission gas thermal conductivity model: coupling atomic, meso and continuum level simulations. Journal of Nuclear Materials, 2013, 440(1‒3): 193–200

[41] HouY, HuangY, SunF, Guo M. Fractal analysis on asphalt mixture using a two-dimensional imaging technique. Advances in Materials Science and Engineering, 2016, 2016: 8931295

[42] GeZ, WangY, HouY, Sun W, SunR . Evaluation of fracture in mortar subject to tension loading using phase field model and three point bending test.Materials and Design, 2015, 86: 121–128

[43] HouY, SunW, HuangY, Ayatollahi M, WangL , ZhangJ. Diffuse interface model to investigate the asphalt concrete cracking subjected to shear loading at a low temperature. Journal of Cold Regions Engineering, 2017, 31(2): 04016009

[44] HouY, WangL, WangD, Liu P, GuoM , YuJ. Characterization of bitumen micro-mechanical behaviors using AFM, phase dynamics theory and MD simulation. Materials, 2017, 10(2): 208

[45] HouY, GuoM, GeZ, WangL, SunW. Mixed-Mode I-II cracking characterization of mortar using phase-field method.Journal of Engineering Mechanics, 2017, 143(7): 04017033

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg