Cement paste is a material with heterogeneous composite structure consisting of hydrated and unhydrated phases at all length scales that varies depending upon the degree of hydration. In this paper, a method to model cement paste as a multiphase system at molecular level for predicting constitutive properties and for understanding the constitutive mechanical behavior characteristics using molecular dynamics is presented. The proposed method creates a framework for molecular level models suitable for predicting constitutive properties of heterogeneous cement paste that could provide potential for comparisons with low length scale experimental characterization techniques. The molecular modeling method followed two approaches: one involving admixed molecular phases and the second involving clusters of the individual phases. In particular, in the present study, cement paste is represented as two-phase composite systems consisting of the calcium silicate hydrate (CSH) phase combined with unhydrated phases tricalcium silicate (C_{3}S) or dicalcium silicate (C_{2}S). Predicted elastic stiffness constants based on molecular model representations employed for the two phases showed that, although the individual phases have anisotropic characteristics, the composite system behaves as an isotropic material. The isotropic characteristics seen from two-phase molecular models mimic the isotropic material nature of heterogeneous cement paste at engineering scale. Further, predicted bulk modulus of the composite system based on molecular modeling is found to be high compared to the elastic modulus, which concurs with the high compression strength of cement paste seen at engineering length scales.

## Introduction

Concrete is a heterogeneous composite material with random structure from macro to molecular scale [1]. The porosity and random composition at each scale is responsible for resultant characteristics that are seen on the properties of concrete [1–3]. Cement paste is the binder phase in concrete and contributes significantly to concrete strength [4]. Cement paste is a hierarchical, multiscale material where molecular level features from hydration and evolving microstructure play a key role on the engineering scale characteristics. Therefore, understanding the interactions and properties of cement paste not only at engineering scale but also at micron and at the submicron scale is essential to future developments and advancements that would allow tailoring of the mechanical properties of concrete by altering the nanoscale and material chemistry configuration of cement paste.

Cement paste is a composite material formed due to hydration when cement clinker is mixed with water. Its structure contains unhydrated phases, voids, and hydrated phases with randomness that is generally on the order of micrometers. The random composite nature of the hydrated and unhydrated phases exists and extends even at the nanometer scale [1]. The principal components in the most common types of unhydrated Portland cement are tricalcium silicate (C_{3}S), dicalcium silicate (C_{2}S), tricalcium aluminate (C_{3}A), and tetracalcium alumino ferrite (C_{4}AF), with a chemical composition of about 37–71%, 6–36%, 0–14%, and 4–18%, respectively. Additional common components in Portland cement are magnesium oxide (MgO) and sulfur trioxide (SO_{3}) [5].

The hydration of cement clinker is a complex process, in which multiple physical–chemical reactions occur, while the cement paste is setting and hardening. At initial stages, needle-shape crystals of ettringite (calcium sulfoaluminate hydrate) appear from the combination of calcium, sulfate, aluminate, and hydroxyl ions. Then, crystals of calcium hydroxide (CH) and fibers of calcium silicate hydrate (CSH) fill the spaces occupied by water and the dissolving cement particles. Depending on the alumina-sulfate ratio of the cement, ettringite may become instable and transforms into monosulfate hydrate (C_{6}A_{3}H_{32}) [6]. The main reaction product of cement hydration that is primarily responsible for its strength is CSH gel, forming approximately 50–60% of the volume of solids in cement paste. The nucleation of CH, that is approximately 20–25% of the volume of solids, takes place generally as massive crystals (up to many micrometers in length); but, it can also be found admixed with CSH at micro- and nano-scale [7]. Unhydrated C_{3}S and C_{2}S remain mixed with the CSH and other hydrated phases in distinct proportions depending on the degree of hydration, the initial composition, and the water cement ratio [7–10]. Due to the nature of the complex hydration process, the structure and properties of cement paste continuously transform during the entire life of the cement components [3,11] However, for engineering purposes, the properties of cement after 28 days of hydration are generally considered for design [11].

Multiscale predictive modeling methods that link the microscopic distribution of hydrated and unhydrated phases to the mechanical behavior of cement paste at macroscale have been reported [12–20]. Most of these methods are based on the definition of a representative volume element (RVE) to describe the morphology of the system. The RVE is established by continuum micromechanics, representing microheterogeneous regions as macrohomogeneous RVEs. To apply continuum micromechanics theories, phases or regions are considered as geometrical domains that may not necessarily correspond to material chemistry phases, i.e., they do not necessarily follow the molecular structures. The effective homogeneous properties of the RVEs are derived from homogenization methods assuming individual isotropic phases. Isotropic properties of individual cement paste phases are generally obtained as isotropic averages from experimental methods as reported in the literature [18,21–23].

At a lower-length scale level, experimental techniques such as nanoindentation [24–26], environmental scanning electron microscopy [27], and ultrasonic pulse velocity [22] have been used in an attempt to understand the variation in the local mechanical properties and its correlation to the phase distribution and the interfacial zones in cement pastes. In these techniques, the elastic properties can be considered to be isotropic in the sense of a statistical average that requires multiple points to be characterized to calculate the averages. The experimental results show that the observed mechanical behavior depends on the localized chemical composition of the small material region used in characterization. However, the effect of distribution and interaction of multitude of phases at nanoscale level on the engineering scale characteristics of cement paste is not yet fully understood. The present work proposes and presents a novel method to study the constitutive mechanical properties and understand characteristics of cement paste modeled as a molecular system formed by a combination of hydrated and unhydrated material chemistry phases, using molecular dynamics simulations. Prior work has primarily reported only such constitutive material predictions of individual single molecular phases for this cement paste composite system [28–35].

## Computational Model of Cement Paste as Multiphase Molecular Systems

### Model Description.

Characteristics and properties of CSH phase are considered to be one of the most dominant factors influencing the behavior of cement paste. Present work thus focused on two-phase systems of cement paste, consisting of CSH as the main hydrated phase, in combination with unhydrated phases C_{3}S and C_{2}S. Modeling nanoscale CSH features presents two big challenges due to the complex nature of the molecular CSH structure, namely: (i) it is a poor crystalline material, and (ii) different molecular representations of CSH phases can be considered. However, in mature cement paste with sufficient hydration, most of the CSH phase has been considered to be analogous to the crystalline mineral jennite [36,37]. Therefore, the molecular structure of mineral jennite is also followed in this present study to represent the CSH phase. The C_{3}S structure is represented by its triclinic form *T1*, stable at room temperature, and the C_{2}S structure corresponds to the monoclinic beta structure. The atomic configuration of each of the phases was obtained from the American Mineralogist Crystal Structure Database (AMCSD) [38]. Table 1 lists the cell parameters of the different molecular structures. Two-phase composite systems of cement paste were created using different mass fractions of CSH with C_{3}S or C_{2}S. The amount of unhydrated components used in the composite molecular configurations was varied to be between 56% and 92% for C_{3}S and 18–66% for C_{2}S. Molecular configuration of each of these single material structures is shown in Fig. 1.

### Construction of Multiphase Molecular Representations of Cement Paste.

To take into account the material chemistry of the individual phase components and the phase boundary effects, cement paste is modeled as two-phase molecular structures consisting of CSH and any one of the unhydrated components. Two approaches were used to generate the two-phase molecular representations; they were labeled “approach A” and “approach B” and are illustrated in Fig. 2. In approach A, the individual molecular structures of two phases are placed inside a scaled box of CSH jennite molecule, which has a total volume equal to the sum of the volumes of individual phases. Then, the scaled box was replicated three times in each of the principal directions to create a supercell. In the second approach, approach B, each individual phase is replicated three times in each principal direction creating individual supercells with larger molecular configurations, and then both are placed into a scaled box, whose total volume is equal to the sum of the volumes of the replicated cells. The density of the composite system forming the two-phase molecular representation of the cement paste is the same in both approaches and follows the rule of mixtures for composite materials.

A result of the two methods employed for the present two-phase molecular representation is that the simulation box in both cases is the same, but the structural configurations of the molecular material systems are different. Approach A represents admixed molecular phases, and approach B involves clusters of the individual phases. To minimize the size and boundary effects, as well as to estimate the macroscale elastic properties, periodic boundary conditions were applied.

### Parametrization of the Molecular Dynamics (MD) Simulations.

Initial atomic positions of the molecular systems phases were taken from AMCSD. The discover module of materials studio [42] MD analyzer and the condensed phase optimized molecular potentials for atomistic simulation studies (COMPASS) force field [43–45] were used to model the interatomic interactions in the two-phase molecular systems of cement paste. It was observed that the long range nonbonded interactions play an important role in the total energy of the cement paste phases. To account for the effect of these interactions, the Ewald summation technique was used, and to assure the charge neutrality of the systems a homogeneous background [46] was applied when required in MD analysis.

Geometrical optimization to obtain a minimal energy configuration of the molecular system was performed using the steepest descent method for an initial coarse minimization followed by the conjugate gradient technique to achieve a fine minimization [47]. Convergence in energy is taken to be obtained when the total energy difference between steps was less than 1 × 10^{−3 }kcal/mol. The temperature and pressure in the systems were equilibrated at 298 K and 1 atm using three sets of simulations. Andersen thermostat and Parrinello–Rahman barostat [48] were used in all simulations. MD modeling analysis methodology involved:

- (1)
Thermal equilibrium at 298 K using a canonical ensemble and a dynamic analysis duration of 100 ps with a 1 fs time-step.

- (2)
Isothermal–isobaric equilibration for 40 ps with a time-step of 0.5 fs.

- (3)
Isothermal–isobaric equilibration for 120 ps with a time-step of 1 fs.

For infinitesimal deformation, the stress–strain relation is considered to be linear, and the stiffness constants can be calculated according to Hooke's law, Eq. (1) [50]. This analysis was performed for the last 20 fs (20 steps) of the isobaric–isothermal equilibration in MD analysis methodology. Those 20 molecular configurations correspond to equivalent systems at the same thermodynamic state. Therefore, statistical averages of the stiffness constants were calculated and used in obtaining predicted values of mechanical stiffness.

## Constitutive Mechanical Properties and Characteristics of Two-Phase Molecular Structures of Cement Paste

### Stiffness Constants of Two-Phase Structures of Cement Paste.

The components of the stiffness tensor following Voigt notation [51] for the individual phases, calculated according to the MD methodology and parametrization presented in Sec. 2, are shown in Table 2. Stiffness tensor constant values close to zero are not shown. A material is considered to have isotropic characteristics if the relations between the components of the stiffness tensor satisfy:

- (i)
*C*_{22}and*C*_{33}are approximately equal to*C*_{11}. - (ii)
*C*_{13}is approximately equal to*C*_{23}and*C*_{12}. - (iii)
*C*_{14},*C*_{15},*C*_{16},*C*_{24},*C*_{25},*C*_{26},*C*_{34},*C*_{35},*C*_{36},*C*_{45},*C*_{46}, and*C*_{56}are close to zero. - (iv)
*C*_{44}is approximately equal to $(C11\u2212C12)/2$.

The predicted elastic stiffness constants presented in Table 2 for the three individual phases indicate that the C_{3}S phase satisfies characteristics to be considered as an approximately isotropic material. The elastic stiffness constants for C_{2}S and CSH phases however do not fully meet the criteria for being considered as isotropic materials.

The constitutive stiffness tensor values of the two-phase systems of cement paste predicted from MD analysis following the two-phase molecular representation approaches (approach A and approach B) indicate that these systems fulfill the conditions to be considered as isotropic materials for all the different compositions of hydrated phase CSH with one of the unhydrated phases C_{3}S and C_{2}S. This behavior is presented in Figs. 3 and 4 for systems composed of CSH and C_{3}S or C_{2}S at different percentages. Stiffness tensor constants values for the molecular systems composed by (*a*) 65.6% C_{3}S and 34.4% CSH and (*b*) 56.1% C_{2}S and 43.9% CSH for both approximations are shown in Table 3.

Clearly, two-phase molecular structure representations of cement paste show an isotropic material behavior based on the predicted stiffness constants. It can be noted that the commonly used homogenization technique indicates that the isotropic characteristics of cement paste are due to the random dispersion of the polycrystalline phases [12,20]. Results from present MD-based constitutive stiffness tensor predictions show clearly isotropic material characteristics, even with two phases combined in different proportions. The isotropic characteristics of the cement paste from the two-phase molecular composite structures are primarily derived from the interaction of phases, and not merely due to isotropic characteristics of individual components. It is evident from the present results that molecular level multiphase (in particular two-phase) representations and modeling techniques noticeably capture the effective material mechanics characteristics of cement paste observed at engineering scale.

#### Isotropic Elastic Properties From Two-Phase Molecular Structures of Cement Paste.

The Hashin–Shtrikman approximation was used to determine upper and lower limits [52], and the Voigt–Reuss–Hill approximation [53] was used to calculate the average of the effective isotropic elastic properties from the predicted stiffness constants of the two-phase systems. The results indicate that there was little to no difference between the Hashin–Shtrikman upper and lower bounds for all properties, with only a 0.5% maximum difference occurring between bounds for bulk modulus. These limits also compare quite closely to the Voigt–Reuss–Hill average for all the two-phase molecular representation systems studied, with a maximum deviation of only 2.7% from the Hashin–Shtrikman bounds for the bulk modulus.

Additionally, theoretical upper and lower limits for the bulk compressive modulus *K* of the two-phase systems were calculated based on the properties of individual components using the rule of mixtures for composite materials following Eqs. (2) and (3) [54]. In these equations, *V* is the volume fractions, *G* is the shear modulus, and the subscripts *m* and *n* refer to each of the components of the composite system (CSH—C_{3}S/C_{2}S). The isotropic bulk and shear modulus of each component were calculated from the stiffness constants presented in Table 2 following the Voigt–Reuss–Hill approximation. These elastic properties values of each individual phase are in agreement with data reported from computational and experimental studies [29,30,33,34,55–57] and are shown in Table 4

The estimated bulk modulus following Voigt–Reuss–Hill approximation and the calculated upper and lower limits from the rule of mixtures for the two-phase systems at various mass fractions of (*a*) C_{3}S and CSH and (*b*) C_{2}S and CSH are shown in Fig. 5. The predictive bulk modulus obtained directly from MD analysis shows a reasonably good agreement with the rule of mixtures, with some values outside the upper and lower bounds.

It is to be noted that the rule of mixture is only based on composition and does not take into account any material phase distribution that plays a critical role in the associated interaction energies and resultant properties and characteristics at the molecular chemistry level, as well as their engineering scale behavior. The two-phase molecular representation of cement paste in the present work effectively accounts for the interaction energies between the two phases and the distribution of the material phase components at molecular level.

The predicted bulk modulus is found to be slightly larger for the approach A representing admixed phases than for approach B which corresponds to clusters of individual phases. This increment can be attributed to the difference in the intermolecular energies between different phases. Clearly, there is more interphase boundaries present in the molecular representation of the two-phase systems employed in approach A compared to the approach B.

## Conclusions

A novel methodology to represent cement paste, consisting of both hydrated and unhydrated material constituents, as a multiphase molecular system at nanoscale level has been proposed and demonstrated. Present study considered a two-phase system comprised of primary hydrated component CSH with one of the unhydrated phase C_{3}S or C_{2}S. Molecular model configuration based on mineral jennite representation for CSH has been followed. The associated stiffness matrix and mechanical properties were calculated using two different approaches employed to represent the molecular structure of two-phase composite systems of cement paste. Predicted stiffness constants obtained from molecular dynamics modeling for the two-phase cement paste systems were found to show isotropic material characteristics. These isotropic characteristics from the two-phase molecular model representation were found to be independent of the isotropy or anisotropy of the individual material phase components. All the systems predicted a high bulk modulus, indicating the cement resistance to compression. It is evident from the present results that molecular level multiphase (two-phase in the present work) representations and modeling techniques capture the effective material mechanics characteristics of cement paste observed at engineering scale; can be extended to model and predict engineering scale characteristics of heterogeneous materials.

Mechanical behavior obtained from low length scale experimental characterization techniques such as nanoindentation depends upon the localized chemical composition of the material and interfaces present in the small loading region. This paper presented a framework for two-phase molecular level models suitable for predicting constitutive properties of heterogeneous cement paste and could provide potential for better comparisons with low length scale experimental characterization techniques.

Material behavior and characteristics at engineering scale can be considered to be based on homogenized scale effects from underlying lower-length scale features; effects of morphology and material chemistry. There is a tremendous interest within the material and mechanics research community in predictive constitutive material models; understanding material behavior and characteristics starting from molecular model representations, and using MD-based modeling. This paper presented and demonstrated two-phase MD-based molecular model representations and modeling methodology that captures heterogeneous features of cement paste. Current results show that these multiphase molecular model representations are effective to predict and mimic the engineering scale characteristics of cement paste and are potentially extendable to other heterogeneous materials.

## Acknowledgment

The authors acknowledge the support by the U.S. Army Research Office under a cooperative agreement award Contract No. W911NF-11-2-0043 (Program Manager: Dr. Joseph Myers) and the U.S. Army Engineering Research and Development Center's Military Engineering Basic/Applied “MMFP” Research Program.

## Nomenclature

*a, b, c*=lattice constant lengths (Å)

*C*=_{ij}the

*ij*th individual stiffness constants in Voigt notation (GPa)*C*=_{ijkl}the

*ijkl*th component of the stiffness tensor- C
_{2}S =dicalcium silicate

- C
_{3}A =tricalcium aluminate

- C
_{3}S =tricalcium silicate

- C
_{4}AF =tetracalcium alumino ferrite

*G*=_{m}effective isotropic shear modulus of hydrated phase

*m**G*=_{n}effective isotropic shear modulus of unhydrated phase

*n**K*_{lower}=theoretical lower limit of the bulk modulus for a composite system

*K*=_{m}effective isotropic bulk modulus of hydrated phase

*m**K*=_{n}effective isotropic bulk modulus of unhydrated phase

*n**K*_{upper}=theoretical upper limit of the bulk modulus for a composite system

- MgO =
magnesium oxide

- SO
_{3}=sulfur trioxide

*V*=_{m}volume fraction of hydrated phase

*m**V*=_{n}volume fraction of unhydrated phase

*n**α, β, γ*=lattice constant angles (deg)

*ε*=_{kl}the

*kl*th component of the strain tensor*σ*=_{ij}the

*ij*th component of the stress tensor (GPa)