Computation of Elastic Properties of Portland Cement Using Molecular Dynamics

Wu, Weidong; Al-Ostaz, Ahmed; Cheng, Alexander H.-D.; Song, Chung R.

doi:10.1061/(ASCE)NM.2153-5477.0000026

Free access

TECHNICAL PAPERS

May 16, 2011

Computation of Elastic Properties of Portland Cement Using Molecular Dynamics

Authors: Weidong Wu [email protected], Ahmed Al-Ostaz, M.ASCE, Alexander H.-D. Cheng, M.ASCE, and Chung R. Song, M.ASCEAuthor Affiliations

Publication: Journal of Nanomechanics and Micromechanics

Volume 1, Issue 2

https://doi.org/10.1061/(ASCE)NM.2153-5477.0000026

PDF

Abstract

There is a growing interest in relating nanostructures to the macro properties of engineering materials such as composites and cement materials. Better understanding of structure and elastic properties of nanoparticles in concrete by modeling and experiment could lead to nanoengineered concrete with much better performance and energy efficiency. In this study, the molecular dynamics (MD) atomistic simulation technique was applied to study the elastic properties of major portland cement compounds (i.e., alite, belite, and aluminate). Applicability of three commonly used force fields: COMPASS, Universal force field (UFF), and Dreiding were evaluated in the MD simulation. The combination of different simulation cell sizes and force fields was investigated. MD simulation results of cement were comparable to the experimental data. The results could be used as nanoparticle properties for multiscale modeling of concrete, cementitious composites, and aggregate.

Introduction

There is a growing interest (Bernard et al. 2003; Feng and Christian 2007; Haecker et al. 2005; Ulm et al. 2004; Wu et al. 2009) in relating the nanostructure to the macro properties of engineering materials such as concrete and cement paste. Modeling of concrete gives a better understanding of its microstructure and is a time and cost effective way to improve its performance. It is essential in the modeling to have knowledge of the accurate properties of the different constituents (residual cement, hydrated cement phases, sand, aggregate, etc.) involved in the final micro- and nanostructure of concrete.

Efforts in determining lastic properties for cement both experimentally and numerically can be found in current literatures (Boumiz et al. 1997; Granju 1987; Manzano et al. 2009; Velez et al. 2001). Boumiz et al. (1997) showed that elastic modulus of unhydrated cement is about 117.6 GPa. Granju (1987) used polished blocks to measure the elastic moduli of cement clinkers, obtaining scattered results that are in between 60 and 300 GPa. Velez et al. (2001) applied nanoindentation technology to measure the elastic modulus and hardness of constituents of portland cement clinker. They found that at the microscopic scale, the elastic moduli of these phases range between 125 GPa for

C_{4} AF (4. CaO . {Al}_{2} O_{3} . {Fe}_{2} O_{3})

to 145 GPa for

C_{3} A (3. CaO . {Al}_{2} O_{3})

.

The work most closely related to the present study of the constituent parts of cement is Manzano et al. (2009). They used force-field atomistic methods to obtain the elastic properties of the main varieties, including alite, belite, portlandite, and C-S-H structurally related mineral crystals tobermorite and jennite. The Young’s modulus they obtained for alite and belite are 138.9 GPa and 137.9 GPa, respectively. The simulations were carried out by General Utility Lattice Program (GULP) code (Gale 1997) which is essentially a lattice dynamics (LD) rather than a molecular dynamics method. Their work did not address some of the important issues in MD and force-field simulation such as the choice of force-field and simulation cell sizes. Moreover, the properties of aluminate (

C_{3} A

) existing in cement was not reported in their paper. Other than cement properties computation, molecular dynamics was also applied to the study of some cement materials issues such as structure and properties of C-S-H (Manzano et al. 2007; Faucon et al. 1997) and water-tobmorite interface (Kalinichev et al. 2007). Pellenq et al. (2009) developed a realistic cement hydrated product C-S-H molecular model, which is important since current molecular dynamics study on C-S-H is mostly based on minerals Tobermorite or Jennite.

Force fields are a fundamental issue underlying all atomistic simulations and are the most important requirement for any molecular mechanics simulation (Cygan and Kubicki 2001). A proper choice of force field (Burkert and Allinger 1989; Lifson and Warshel 1968) that uses a set of empirical formulas to describe the interatomic interactions among various atoms in a molecule or a group of molecules is a crucial step in MD simulation. In the present paper, three different force fields: COMPASS (Condensed-phase Optimized Molecular Potentials for Atomic Simulation Studies); UFF (Universal force field) (Rappe et al. 1992); and Dreiding (Mayo et al. 1990) were applied to explore their validity in cement simulation. COMPASS is incorporated in both Discover and Forcite Plus, which are two atomistic simulation modules in the Materials Studio software used for the present study. UFF and Dreiding are only available in Forcite. Discover is a powerful atomistic simulation code which has a range of well validated force fields. It may be used to study relationships between structural and molecular behavior, provide insight into key molecular interactions, and predict important properties of solids, liquids, and gases. Forcite Plus is a classical molecular mechanics simulation tool which allows energy calculations, geometric optimizations, dynamics simulations, and mechanical properties calculation. Song et al. (2007) found that the size of the MD simulation cell may have an effect on the mechanical properties results when they used MD to study crystal quartz properties. They obtained a negative Poisson’s ratio for a single quartz unit cell and homogenized properties were captured by increasing the size of the simulating cell. Therefore, this study considered different sizes of the simulation supercell (using multiunit cells to create a supercell) to better represent the structure of the crystals.

This paper proposes a molecular dynamics approach to calculate elastic properties (specifically bulk, shear, Young’s moduli, and Poisson’s ratios) of cement constituents alite, belite, aluminate mineral crystal. In conjunction with the present study, the following sections first briefly describe the simulating crystal structures of

C_{3} S (3. CaO . {SiO}_{2})

,

C_{2} S (2. CaO . {SiO}_{2})

, and

C_{3} A

; the approach of MD simulation which includes introduction of MD. The comparison of the three force fields and the MD conditions adopted in this work were introduced next. Simulation results and discussions were given in the end.

Structures of Simulating Crystals

Cement has four main clinker minerals: alite, belite, aluminate, and brownmillerite (Taylor 1997).

C_{3} S

is formed above 1,250°C by a reaction of

C_{2} S

and

C

.

C_{3} S

undergoes several polymorphic phase transformations at different temperatures.

C_{2} S

also has several polymorphs:

α - C_{2} S

,

β - C_{2} S

and

γ - C_{2} S . C_{3} A

is a typical random composite made from mixing cement, aggregate, and sand.

C_{3} S

is the most important compound and may exist in cement for 55–65% of total mass. The crystal structure of

C_{3} S

(Golovastikov et al. 1975) is built up from

{SiO}_{4}^{4 -}

tetrahedral with

{Ca}^{2 +}

in the corners and oxygen in the center (Fig. 1). Dicalcium silicate (Mideley 1952) or

C_{2} S

is another important constituent of portland cement. The structure is made up of isolated

{SiO}_{4}

tetrahedral and Ca ions, which look like strings of alternating Ca ions and tetrahedral since four of the eight Ca ions are positioned alternately above and below

{SiO}_{4}

tetrahedral in the

y

direction (Fig. 2). Compared with

C_{3} S

and

C_{2} S

,

C_{3} A

compound (Mondal and Jeffery 1975) does not exhibit polymorphism. The latter is a cubic structure built from

{Ca}^{2 +}

ions and highly puckered rings of six

{AlO}_{4}

tetrahedral (Fig. 3). The crystallographic data for the four single crystals used in this research are given in Table 1.

Fig. 1. Unit cell crystal structure of $C_{3} S$ ; projection view along the (100) direction, displaying calcium (Ca), silicate (Si), aluminum (Al), and oxygen (O)

Fig. 2. Unit cell crystal structure of $C_{2} S$ : projection view along the (100) direction

Fig. 3. Unit cell crystal structure of $C_{3} A$ : projection view along the (001) direction

Table 1. Cell Parameters of

C_{3} S

,

C_{2} S

, and

C_{3} A

Crystals

Crystal system	$C_{3} S$	$C_{2} S$	$C_{3} A$
Space group	P-1	$P 21 / n$	Pa3
$a$ [ $Å$ ]	11.67	5.48	15.263
$b$ [ $Å$ ]	14.24	6.76	15.263
$c$ [ $Å$ ]	13.72	9.28	15.263
$α$ [ $o$ ]	105.5	90	90
$β$ [ $o$ ]	94.33	94.33	90
$γ$ [ $o$ ]	90	90	90

Methods

Introduction of MD Simulation

Two atomistic simulation techniques—lattice dynamics, used by Manzano et al. (2009), and MD, adopted in this research—may be applied to model crystal structure at atomic level. LD, which was originally developed by Born and Huang (1954), uses an analytical approach to solve Newton’s equation of motion, followed by a statistical mechanical treatment. Compared with LD, MD uses a numerical solution to the equation of motion. The positions and velocities of atoms in the simulating system are also not constrained to solids as they are in LD, and all of the particles are effectively involved in time-dependent motion. Both methods are capable of modeling solids, and although MD requires more computational time, it has the potential to be more accurate. For more details about the differences between these two atomic simulation techniques, refer to Parker et al. (2001).

Molecular dynamics explore the macroscopic properties through simulations at the microscopic level, including atomic positions and velocities by statistical mechanics. MD simulation begins with energy minimization, and after a structure is created, it is refined to produce a stable conformation. However, geometry optimization can be allowed along with the energy minimization in certain molecular dynamic simulations. Minimization is an iterative procedure in which the coordinates of the atoms and possibly the cell parameters, are adjusted so that the total energy of the structure is reduced to a local minimum. Although the geometry of a system can be optimized through energy minimization process, thermal dynamic variables temperature and pressure control are only incorporated in dynamics simulation. The kinetic energy introduced by temperature control can help the molecular motions to overcome potential barriers since there is usually more than one minimum (Fig. 4) under periodic boundary conditions (Cygan and Kubicki 2001). To find the global minimum, temperature is periodically increased from an initial temperature (room temperature in this study) to a midcycle temperature e.g., 500 k, and then decreased again to the initial temperature. This strategy is called anneal dynamics, and it can be used to capture a global minimum.

The energy-optimized structure was then used as the initial structure for dynamics simulation. An initial velocity was assigned to each atom in the simulating system. It was determined from the Maxwell-Boltzmann distribution, and then each of these was scaled to heat the system to the target temperature. The forces acting on each atom, their velocities, and their positions are derived for all subsequent time steps based on the given potential through solving Newton’s equation of motion. A large number of iterations are usually needed to bring the system to an equilibrium configuration. To ensure thermodynamic equilibrium, this study constantly monitored thermodynamics quantities, such as energy and temperature, versus time during simulations. When equilibration was achieved, these quantities fluctuated around their averages, which remained constant over time.

After the dynamics simulation had been performed, the resulting deformed molecular structure was analyzed to determine elastic constants. Elastic constants of the final atomic configuration are computed using the static approach suggested by Theodorou and Suter (1986). The generalized expressions of elastic constants relating the various components of stress and strain are defined as

C_{lmnk} = \frac{\partial σ_{l m}}{\partial σ_{n k}} |_{T, ε_{n k}} = \frac{1}{V_{0}} \frac{\partial^{2} A}{\partial ε_{l m} \partial ε_{n k}} |_{T, ε_{l m}, ε_{n k}}

(1)

where

C_{lmnk}

= a fourth-order stiffness tensor, caused by the symmetry of the stress and strains, and the conservation of energy, this fourth-order may be reduced to the second-order tensor;

σ_{l m}

= second-order stress tensor;

ε_{l m}

and

ε_{n k}

= second-order strain tensors;

A

= Helmholtz free energy; and

V_{o}

= volume of the simulation cell in the undeformed configuration. The corresponding isotropic polycrystalline elastic moduli can be calculated based on the corresponding single-crystal elastic constants by Voigt-Reuss-Hill (VRH) approximation (Chung and Buessem 1968; Hill 1952). According to Hill (1952), VRH approximation can be simply expressed as

K_{V R H} = (K_{V} + K_{R}) / 2

(2)

and

G_{V R H} = (G_{V} + G_{R}) / 2

(3)

where

V

and

R

refer to the Voigt and the Reuss averages. They can be derived from elastic stiffness constants

C_{i j}

and elastic compliance constants

S_{i j}

by the following Eqs. (4) - (7); also

C_{i j}

and

S_{i j}

in these expressions are no longer tensors as in Eq. (1).

Voigt average:

9 K_{V} = (C_{11} + C_{22} + C_{33}) + 2 (C_{12} + C_{23} + C_{31})

(4)

15 G_{V} = (C_{11} + C_{22} + C_{33}) - (C_{12} + C_{23} + C_{31}) + 3 (C_{44} + C_{55} + C_{66})

(5)

Reuss average:

1 / K_{R} = (S_{11} + S_{22} + S_{33}) + 2 (S_{12} + S_{23} + S_{31})

(6)

15 / G_{R} = 4 (S_{11} + S_{22} + S_{33}) - 2 (S_{12} + S_{23} + S_{31}) + 3 (S_{44} + S_{55} + S_{66})

(7)

Force Fields

Unfortunately, there is no universal guidance that may be used to help choose a right force field for a certain mineral crystal for MD simulation; therefore, it is desirable to try various possible force fields. To address differences and to better understand the three force fields used in this study, a brief description of their functional forms of potential energy is discussed below.

COMPASS Force Field

The COMPASS force field is a high-quality general ab initio force field that has been widely used in simulations of liquids, crystals, and polymers. COMPASS is the first ab initio force field that has been parameterized and validated using condensed-phase properties, as well as various ab initio and empirical data for molecules in isolation. Consequently, this force field enables accurate and simultaneous prediction of structural, conformational, vibrational, and thermophysical properties for a broad range of molecules, and under a wide range of conditions of temperature and pressure. The functional forms of the COMPASS (Sun 1998) force field are

E_{total} = \sum_{b} [K_{2} (b - b_{0})^{2} + K_{3} (b - b_{0})^{3} + K_{4} (b - b_{0})^{4}] + \sum_{θ} [K_{2 θ} (θ - θ_{0})^{2} + K_{3 θ} (θ - θ_{0})^{3} + K_{4 θ} (θ - θ_{0})^{4}] + \sum_{ϕ} [K_{1 ϕ} (1 - \cos ϕ) + K_{2 ϕ} (1 - \cos 2 ϕ) + K_{3 ϕ} (1 - \cos 3 ϕ)] + \sum_{χ} K_{2 χ} (χ - χ_{0})^{2} + \sum_{b, θ} K_{b θ} (b - b_{0}) (θ - θ_{0}) + \sum_{b, ϕ} (b - b_{0}) [K_{1 b} \cos ϕ + K_{2 b} \cos 2 ϕ + K_{3 b} \cos 3 ϕ] + \sum_{θ, ϕ} (θ - θ_{0}) [K_{1 θ ϕ} \cos ϕ + K_{2 θ ϕ} \cos 2 ϕ + K_{3 θ ϕ} \cos 3 ϕ] + \sum_{b, θ} (θ^{'} - θ_{0}^{'}) (θ - θ_{0}) + \sum_{θ, ϕ} K_{θ, ϕ} (θ^{'} - θ_{0}^{'}) (θ - θ_{0}) \cos ϕ + \sum_{i, j} \frac{q_{i} q_{j} e}{r_{i j}} + \sum_{i, j} ε_{i j} [2 (\frac{r_{i j}^{0}}{r_{i j}})^{9} - 3 (\frac{r_{i j}^{0}}{r_{i j}})^{6}]

(8)

where

K

,

H

,

F

,

A

and

B

= constants;

b

= bond distance;

θ

= bond angle;

ϕ

is related to torsion; and

V

= coefficient of cosine-Fourier expansion. Terms in lines five to 11 in the Eq. (8) are cross terms, which are bond or angle distortions caused by nearby atoms. These off-diagonal cross-coupling terms make COMPASS a very unique and attracting force field. Although the COMPASS force field has high accuracy because of its ab initio approach and cross-coupling terms, it is constantly being revised to extend its coverage because input parameters are not available yet for all elements.

UFF

UFF is a purely diagonal, harmonic force field. Bond stretching for UFF is described by a harmonic term, angle bending by a three-term Fourier cosine expansion, and torsions and inversions by cosine-Fourier expansion terms. The van der Waals interactions are described by the Lennard-Jones potential. Electrostatic interactions are described by atomic monopoles and a distance-dependent Coulombic term. The most distinguishing advantage of UFF lies in the fact that it covers a full periodic table, which means when COMPASS or other force fields fail to work on some crystals, switching to UFF is possible. The potential energy of UFF is given by

E_{total}^{U F F} = E_{R} + E_{θ} + E_{ϕ} + E_{ω} + E_{v d w} + E_{e l} = \frac{1}{2} K_{I J} (r - r_{I J})^{2} + K_{I J K} \sum_{n = 0}^{m} C_{n} \cos n θ + K_{I J K L} \sum_{n = 0}^{m} C_{n} \cos n ϕ_{I J K L} + K_{I J K L} (C_{0} + C_{1} \cos ω_{I J K L} + C_{2} \cos 2 ω_{I J K L}) + D_{I J} {- 2 [\frac{x_{I J}}{x}]^{6} + [\frac{x_{I J}}{x}]^{12}} + 332.0637 (\frac{Q_{i} Q_{j}}{ε R_{i j}})

(9)

where

K_{I J}

,

K_{I J K}

, and

K_{I J K L}

= force constants that relate to bond stretch, angel bend, and bond torsion, respectively. The first term represents the bond stretching where only one bond

I J

is involved; the second term indicates an angle bending where two bonds

I J

and

J K

are contained, and the third term implies bonds torsion where three bonds

I J

,

J K

, and

K L

are involved. The symbol

θ

= periodic angle in Fourier expansion;

C_{n}

= expansion coefficient defined by the natural bond angle (

θ_{0}

);

r

= bond length;

r_{0}

= natural bond length;

ω_{I J K L}

= angle between the IL axis and IJK plane;

D_{i j}

= well depth.

x_{I J}

= van der Waals bond length;

Q_{i}

and

Q_{j}

= charges in electron units;

R_{i j}

= distance in Angstroms, and

ε

= dielectric constant.

The angle bending, torsions and inversions are all described by cosine-Fourier expansion terms. The expression does not contain cross terms that appeared in the COMPASS potential energy.

Dreiding

The Dreiding force field is purely diagonal, with harmonic valence terms and a cosine-Fourier expansion torsion term. The van der Waals interactions are described by the Lennard-Jones potential. Electrostatic interactions are described by atomic monopoles and a screened (distance dependent) Coulombic term. Hydrogen bonding is described by an explicit Lennard-Jones 12–10 potential. The potential energy of Dreiding force field is expressed as

E_{total}^{Dreiding} = E_{v a l} + E_{n b} = (E_{B} + E_{A} + E_{T} + E_{I}) + (E_{v d w} + E_{Q} + E_{h b}) = \frac{1}{2} k_{e} (R - R_{e})^{2} + \frac{1}{2} C_{I J K} [\cos θ_{I J K} - \cos θ_{J}^{0}]^{2} + \frac{1}{2} V_{J K} {1 - \cos [n_{J K} (φ - φ_{J K}^{0})]} + \frac{1}{2} K_{i n v} (ψ - ψ_{0})^{2} + D_{0} [ρ^{- 12} - 2 ρ^{- 6}] + 332.0637 (\frac{Q_{i} Q_{j}}{ε R_{i j}}) + D_{h b} [5 (\frac{R_{h b}}{R_{D A}})^{12} - 6 (\frac{R_{h b}}{R_{D A}})^{10}] \cos (θ_{D H A})

(10)

where

E_{B}

,

E_{A}

,

E_{T}

,

E_{I}

, and

E_{h b}

= energy of bond stretch, bond angle bend, dihedral angle torsion, inversion, and hydrogen bond, respectively;

K

,

C

, and

V

= constants;

D_{0}

= van der Waals well depth;

ρ = R / R_{0}

= scaled distance;

R_{0}

= van der Waals bond length;

Q

= charge;

D_{h b}

,

R_{h b}

, and

R_{D A}

= hydrogen parameters.

The three force fields have similar terms for stretching-van der Waals (COMPASS uses 6–9 type, and Universal and Dreiding use 9–12 type Lennard-Jones type expressions) and electrostatic interactions. Although UFF covers the full periodic table, its energy expression is purely diagonal and harmonic. The COMPASS has more complicated coupling cross terms, unlike the other force fields. The COMPASS is expected to be the best candidate, despite the longer computational time.

MD Conditions Used in This Study

Supercells

1 a \times 1 b \times 1 c

,

2 a \times 2 b \times 2 c

(

a

,

b

,

c

are the unit cell dimensions, corresponding to 1.8 unit cells) that were used in this study’s simulations. (The

3 a \times 3 b \times 3 c

supercell was also used for

C_{2} S

.) The number of atoms in the simulation systems ranges from 28 (for

C_{2} S

unit cell) to 2,112 (for the

2 a \times 2 b \times 2 c

C_{3} A

supercell).

For systems with less than 200 atoms, this study used a smart minimization method which combines a steepest descent, a conjugate gradient, and a Newton-Raphson method that began with the steepest descent method, followed by the conjugate gradient method, and ended with a Newton-Raphson method. In studies that dealt with systems larger than 200 atoms, Newton-Raphson is replaced with the conjugate gradients method. The cutoff distance (Fig. 5) in all the simulations is 12.50Å. There is a high computational cost to fully consider the nonbonded terms, because an atom is bonded to only a few atoms of its neighbors, but it interacts with almost all the atoms in the molecule. Fortunately, the van der Waals, nonbond energy term, declines rapidly—it is typically modeled using a (6–12 Lennard-Jones potential), which implies that the attractive forces vanish with distance as

r^{- 6}

and the repulsive forces as

r^{- 12}

, where

r

= the distance between two atoms. In general, a cutoff radius is used to speed up the calculation so that the atoms separated with a distance greater than the cut off radius have no nonbonded interaction. Ewald summation method is used to determine nonbond (van der Waals interactions and coulomb electrostatic) energies. Energy minimized models were used as initial structures for molecular dynamics simulations that were performed as

N P T

canonical ensembles at

P = 0.001 GPa

(air pressure) and

T = 298 K

(room temperature). A statistical ensemble is a collection of all possible systems which have different microscopic states (atomic positions) but have an identical macroscopic state (pressure, temperature, etc.). The choice of

N P T

ensemble is based on the work done by Kalinichev et al. (2007) and Cygan et al. (2004). Proper temperature and pressure controlling methods should be chosen after the choice of ensemble to generate the correct statistical ensemble. The Andersen thermostat temperature controlling method and the Parrinello barostat pressure controlling scheme were applied in the Discover runs. The Nosé thermostat temperature controlling method and the Berendsen barostat pressure controlling scheme were employed in the Forcite simulations. Another important parameter is time step. The criterion for choosing a time step is that the increments between steps should be small when compared to the average time between molecule collisions. Usually the time step should be approximately one-tenth of the shortest period of motion of molecules (Leach 2001). This study set the time step to 1 fs (femtosecond,

10^{- 15} s

), which is an acceptable compromise between computational cost and accuracy (Parker et al. 2001). Depending on the number of atoms in the system, the dynamic time used in this research ranges between 100 ps (picosecond,

10^{- 12} s

) to 400 ps.

Results

Molecular dynamics with three common force fields, two MD tools, and variously sized supercells were used to simulate single crystals of cement constituents

C_{3} S

,

C_{2} S

, and

C_{3} A

. Eqs. (11) - (12) are sample elastic constants for unit cell

C_{3} S

with a combination of the COMPASS force field and the Forcite Plus MD tool.

Elastic stiffness constants Cij (GPa) of

C_{3} S

:

C_{i j} = [\begin{matrix} 132.0832 & 71.3912 & 49.3632 & 21.6456 & - 6.4990 & 1.5686 \\ 71.3912 & 211.6644 & 124.8335 & - 4.1212 & 0.5856 & - 0.5979 \\ 49.3632 & 124.8335 & 234.2412 & 2.4543 & 5.7832 & - 1.9870 \\ 21.6456 & - 4.1212 & 2.4543 & 57.4192 & 0.1954 & 0.0241 \\ - 6.4990 & 0.5856 & 5.7832 & 0.1954 & 56.9348 & - 0.6179 \\ 1.5686 & - 0.5979 & - 1.9870 & 0.0241 & - 0.6179 & 50.6488 \end{matrix}]

(11)

Elastic compliance constants

S_{i j} (1 / T P a)

of

C_{3} S

S_{i j} = [\begin{matrix} 10.2395 & - 3.3109 & - 0.3845 & - 4.0854 & 1.2521 & - 0.3541 \\ - 3.3109 & 7.9904 & - 3.5783 & 1.9749 & - 0.1028 & 0.0543 \\ - 0.3845 & - 3.5783 & 6.2787 & - 0.3782 & - 0.6413 & 0.2083 \\ - 4.0854 & 1.9749 & - 0.3782 & 19.1155 & - 0.5125 & 0.1197 \\ 1.2521 & - 0.1028 & - 0.6413 & - 0.5125 & 17.7765 & 0.1519 \\ - 0.3541 & 0.0543 & 0.2083 & 0.1197 & 0.1519 & 19.7654 \end{matrix}]

(12)

This study combines Eqs. (11) - (12) with Eqs. (2) - (7) to arrive at the isotropic elastic properties of

C_{3} S

shown in Table 2, column 3. Table 2 provides a summary of the MD computation results of

C_{3} S

obtained by the following combinations: (1) discover with COMPASS force field (D-C); (2) Forcite Plus together with COMPASS force field (F-C); and (3) Forcite Plus with Universal force field (F-U). Compared with the experimental results, this study concludes that the best choice for the simulation of

C_{3} S

is the F-C combination, although D-C also yields good results.

Table 2. MD Simulation Results of

C_{3} S

Properties	$1 a \times 1 b \times 1 c$			$2 a \times 2 b \times 2 c$			Ref. Values	Difference^e
Supercell
MD tools-force field	D-C	F-C	F-U	D-C	F-C	F-U
$E$ (GPa)	168	137	45.4	139	148	6.11	135^a; 117^b; $60 \sim 300$ ^c; 138.9^d	1.5%
$ν$	0.34	0.29	0.26	0.36	0.35	0.47	0.31^b	$- 6.5 %$
$K$ (GPa)	180	110	31.8	171	168	33.3	—
$G$ (GPa)	63	53	18	51	54.8	2.08	—

Note:

E

: elastic modulus;

ν

: Poisson’s ratio;

κ

: bulk modulus;

G

: shear modulus;

D

: Discover;

C

: COMPASS;

F

: Forcite Plus;

U

: UFF.

^a

^b

^c

^d

^e

The errors are computed based on the numbers in bold.

A combination of D-C and F-C was used in this study based on the results of

C_{3} S

, we mainly used combinations of D-C and F-C in the simulation of

C_{2} S

. The computation results for

C_{2} S

are given in Table 3, which again indicated F-C combination is the best choice. The F-U combination underestimates the modulus values in both

C_{3} S

and

C_{2} S

simulation results compared with the reference moduli of 135 GPa for

C_{3} S

and 130 GPa for

C_{2} S

. The MD simulation results for

C_{3} A

are shown in Table 4, and the D-C combination gave the Young’s modulus of 141 GPa, which is only less than 2.8% of the value determined by nanoindentation, while the MD results under different simulation conditions have large differences.

Table 3. Molecular Simulation Results of

C_{2} S

Properties	$1 a \times 1 b \times 1 c$			$2 a \times 2 b \times 2 c$		$3 a \times 3 b \times 3 c$		Ref. value	Difference
Supercell
MD tools-force field	D-C	F-C	F-U	D-C	F-C	D-C	F-C
$E$	285	122	58	265	121	262	121	130^a	$- 6.9 %$
$ν$	0.23	0.19	0.35	0.22	0.19	0.22	0.19	0.31^b	$- 39 %$
$κ$	177	65	63.2	161	64	161	64		—
$G$	116	51	21.4	108	51	108	51		—

^a

Velez et al. 2001.

^b

Boumiz et al. 1997.

Table 4. Molecular Simulation Results of

C_{3} A

Properties	$1 a \times 1 b \times 1 c$		$2 a \times 2 b \times 2 c$	Ref. value	Difference
Supercell
MD tools-force field	F-U	F-D	F-U
$E$	70	141	66.3	145^a	$- 2.8 %$
$ν$	0.31	0.22	0.32	0.31^b	3.2%
$κ$	63	84	59.8	—
$G$	26.6	58	25.2	—

^a

Velez et al. 2001.

^b

Boumiz et al. 1997.

Conclusions and Discussions

This molecular dynamics study on the mechanical properties of major cement constituents showed that the elastic moduli of

C_{3} S

,

C_{2} S

and

C_{3} A

are 137 GPa, 121 GPa, and 141 GPa, respectively. The corresponding Poisson’s ratios are: 0.29, 0.19, and 0.22. These are the best achieved results when compared to the published nanoindentation experimental values. According to these results, this study concludes the following:

1.

A COMPASS force field and Forcite Plus combination are the best choice for the simulation of

C_{3} S

,

C_{2} S

. Dreiding force field with Forcite Plus gave reasonable properties of

C_{3} A

.

2.

From the results obtained from simulating multiple supercell sizes with a limited number of atoms, it concluded that the size of supercells has a minimal effect on the mechanical properties of cement compounds, compared to the effect of chosen force field. However, for some other materials such as quartz (Song et al. 2007), and Multiwall Carbon Nanotube (MWCNT) (Alkhateb et al. 2008), the mechanical properties are dependent on the supercell size.

The different combinations of molecular dynamics simulation engines, even for the same force fields might yield inconsistent results. One of the possible reasons could be that a Parinello barostat is used in the Discover module calculations to control the pressure and stress. The Parinello method allows changes in both the shape and the volume of the cell, so that the internal stress of the system can match the externally applied stress, and it is thus very useful for studying the stress-strain relationship in materials, but a Parrinello barostat may cause large pressure fluctuations. However, a Berendsen barostat is applied in the Forcite Plus simulations (the Parinello is not available in Forcite Plus). A Berendsen pressure control method changes the pressure by altering the coordinates of the particles and the size of the unit cell in the simulation system, thus avoiding the fluctuations brought on in a Parrinello barostat.

MD is a very useful atomic simulation tool in the computation of elastic properties for cementitious materials as illustrated in this study. However, the computational time of MD simulation is tremendous for systems with large number of atoms. Fortunately, faster computer hardware is more affordable, and better MD coding and algorithms are currently available. A real universal force field for the simulation of cementitious materials does not exist; therefore, the simulation is largely dependent on force fields. Thus, the careful selection of force fields in MD simulation is first and foremost to consider; it is always a good practice to choose more than one force field.

Acknowledgments

This work was partially supported by funding received under a subcontract from the Department of Homeland Security-sponsored Southeast Region Research Initiative (SERRI) at the Department of Energy’s Oak Ridge National Laboratory. The authors also would like to thank Miss Hunain Alkhateb for her comments and constructive suggestions.

References

Alkhateb, H., Al-Ostaz, A., and Cheng, A. H.-D. (2008). “Geometric and simulation parameter effects on elastic moduli of multi-walled carbon nanotubes using molecular dynamics approach.” Proc. of the 23rd Technical Conf. of American Society for Composites (CD-ROM), Memphis, TN.

Google Scholar

Bernard, O., Ulm, F.-J., and Lemarchand, E. (2003). “Multiscale micromechanics-hydration model for the early-age elastic properties of cement-based materials.” Cem. Concr. Res., 33(9), 1293–1309.

Crossref

Google Scholar

Born, M., and Huang, K. (1954). Dynamical theory of crystal lattices, Oxford University Press, London.

Google Scholar

Boumiz, A., Sorrentino, D., Vernet, C., and Tenoudji, F. C. (1997). “Modelling the development of the elastic moduli as a function of the degree of hydration of cement pastes and mortars.” Proc. 13 of the 2nd RILEM Workshop on Hydration and Setting: Why does cement set? An interdisciplinary approach, RILEM, Dijon, France.

Google Scholar

Burkert, U., and Allinger, N. L. (1989). “Molecular mechanics.” ACS monograph 77, Oxford University Press, New York.

Google Scholar

Chung, D. H., and Buessem, W. R. (1968). “The Voigt-Reuss-Hill (VRH) approximation and the elastic moduli of polycrystalline ZnO,

Ti O_{2}

, and α-

A l_{2} O_{3}

.” J. Appl. Phys., 39(6), 2777–2782.

Crossref

Google Scholar

Cygan, R. T., and Kubicki, J. D., eds. (2001). “Molecular modeling in mineralogy and geochemistry.” Molecular modeling theory: Applications in the geosciences, Geochemical Society/Mineralogical Society of America, 1–35.

Crossref

Google Scholar

Cygan, R. T., Liang, J.-J., and Kalinichev, A. G. (2004). “Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field.” J. Phys. Chem. B, 108(4), 1255–1266.

Crossref

Google Scholar

Faucon, P., Delaye, J. M., Virlet, J., Jacquinot, J. F., and Adenot, F. (1997). “Study of the structural properties of the C-S-H(I) by molecular dynamics simulation.” Cem. Concr. Res., 27(10), 1581–1590.

Crossref

Google Scholar

Feng, L., and Christian, M. (2007). “Micromechanics model for the effective elastic properties of hardened cement pastes.” Acta Materiae Compositae Sinica, 24(2), 184–189.

Google Scholar

Gale, J. D. (1997). “GULP: A computer program for the symmetry adapted simulation of solids.” J. Chem. Soc., Faraday Trans., 93, 629.

Crossref

Google Scholar

Golovastikov, N. I., Matveeva, R. G., and Belov, N. V. (1975). “Crystal structure of the tricalcium silicate

3 CaO . Si O_{2}

.” Sov. Phys. Crystallogr., 20(4), 441.

Google Scholar

Granju, J. L. (1987). Modélisation des pates de ciments durcies: Caractérisation de l'état d'hydratation, lois d'évolution de la résistance en compression et du module de déformation longitudinale, Université Paul Sabatier de Toulouse (in French).

Google Scholar

Haecker, C. J., et al. (2005). “Modeling the linear elastic properties of portland cement paste.” Cem. Concr. Res., 35(10), 1948–1960.

Crossref

Google Scholar

Hill, R. (1952). “The elastic behavior of a crystalline aggregate.” Proc. Phys. Soc. London Sect. A, 65, 349–354.

Crossref

Google Scholar

Kalinichev, A. G., Wang, J., and Kirkpatrick, R. J. (2007). “Molecular dynamics modeling of the structure, dynamics and energetics of mineral-water interfaces: Application to cement materials.” Cem. Concr. Res., 37(3), 337–347.

Crossref

Google Scholar

Leach, A. R. (2001). Molecular modelling principles and applications, Prentice Hall.

Google Scholar

Lifson, S., and Warshel, A. (1968). “Consistent force field for calculations of conformations vibrational spectra and enthalpies of cycloalkane and n-alkane molecules.” J. Chem. Phys., 49, 5116–5129.

Crossref

Google Scholar

Manzano, H., Dolado, J., and Ayuela, A. (2009). “Elastic properties of the main species present in portland cement pastes.” Acta Mater., 57(5), 1666–1674.

Crossref

Google Scholar

Manzano, H., Dolado, J. S., Griebel, M., and Hamaekers, J. (2007). “A molecular dynamics study of the aluminosilicate chains structure in Al-rich calcium silicate hydrated (C-S-H) gels.” Phys. Status Solidi A, 25(6), 1324–1329.

Crossref

Google Scholar

Mayo, S. L., Olafson, B. D., and Goddard, W. A., III. (1990). “Dreiding: A generic force field for molecular simulations.” J. Phys. Chem., 94, 8897–8909.

Crossref

Google Scholar

Mideley, C. M. (1952). “The crystal structure of

β

dicalcium silicate.” Acta Crystallogr., 5, 307–312.

Crossref

Google Scholar

Mondal, P., and Jeffery, J. W. (1975). “The crystal structure of tricalcium aluminate,

C a_{3} A l_{2} O_{6}

.” Acta Crystallogr. Sect. A, B31, 689–697.

Crossref

Google Scholar

Parker, S. C., de Leeuw, N. H., Bourova, E., and Cooke, D. J. (2001). “Application of lattice dynamics and molecular dynamics techniques to minerals and their surfaces.” Molecular modeling theory: Applications in the geosciences, R. T. Cygan and J. D. Kubicki, eds., Geogchemical Society/Mineralogical Society of America, 63–82.

Crossref

Google Scholar

Pellenq, R. J.-M., et al. (2009). “A realistic model of cement hydrates.” Proc. Nat. Acad. Sci., 106(38), 16102–16107.

Crossref

Google Scholar

Rappe, A. K., Casewit, C. J., Colwell, K. S., Goddard, W. A., III, and Skiff, W. M. (1992). “UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations.” J. Am. Chem. Soc., 114(25), 10024–10035.

Crossref

Google Scholar

Song, C. R., Cho, H., Jung, Y.-H., Cheng, A. H. D., and Al-Ostaz, A. (2007). “Bridging molecular, particulate, and continuum mechanics for geomechanics application.” Proc., Geo-Denver 2007, ASCE, Denver, 18–21.

Google Scholar

Sun, H. (1998). “COMPASS: An ab initio force-field optimized for condensed-phase applications-overview with details on alkane and benzene compounds.” J. Phys. Chem. B, 102, 7338–7364.

Crossref

Google Scholar

Taylor, H. F. W. (1997). Cement chemistry, Thomas Telford, London.

Crossref

Google Scholar

Theodorou, D. N., and Suter, U. W. (1986). “Atomistic modeling of mechanical properties of polymeric glasses.” Macromolecules, 19, 139–154.

Crossref

Google Scholar

Ulm, F.-J., Constantinides, G., and Heukamp, F. H. (2004). “Is concrete a poromechanics material?—A multiscale investigation of poroelastic properties.” Mater. Struct., 37, 43–58.

Crossref

Google Scholar

Velez, K., Maximilien, S., Damidot, D., Fantozzi, G., and Sorrebtino, F. (2001). “Determination by nanoindentation of elastic modulus and hardness of pure constituents of portland cement clinker.” Cem. Concr. Res., 31(4), 555–561.

Crossref

Google Scholar

Wu, W., Al-Ostaz, A., Cheng, A. H.-D., and Song, C. R. (2009). Concrete as a hierarchical structural composite material, Univ. of Mississippi, University, MS.

Google Scholar

Information & Authors

Information

Published In

Journal of Nanomechanics and Micromechanics

Volume 1 • Issue 2 • June 2011

Pages: 84 - 90

Copyright

History

Received: Mar 21, 2010

Accepted: Oct 1, 2010

Published online: May 16, 2011

Published in print: Jun 1, 2011

Permissions

Request permissions for this article.

Request Permissions

Authors

Affiliations

Weidong Wu [email protected]

Research Associate, Dept. of Civil Engineering, Univ. of Mississippi, 203 Carrier Hall, University, MS 38677 (corresponding author). E-mail: [email protected]

View all articles by this author

Ahmed Al-Ostaz, M.ASCE

Associate Professor, Dept. of Civil Engineering, Univ. of Mississippi, 202 Carrier Hall, University, MS 38677.

View all articles by this author

Alexander H.-D. Cheng, M.ASCE

Professor and Dean, School of Engineering, Univ. of Mississippi, 101 Carrier Hall, University, MS 38677.

View all articles by this author

Chung R. Song, M.ASCE

Associate Professor, Dept. of Civil Engineering, Univ. of Mississippi, 218 Carrier Hall, University, MS 38677.

View all articles by this author

Metrics & Citations

Metrics

Citations

Download citation

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

Computation of Elastic Properties of Portland Cement Using Molecular Dynamics

Abstract

Introduction

Structures of Simulating Crystals

Methods

Introduction of MD Simulation

Force Fields

COMPASS Force Field

UFF

Dreiding

MD Conditions Used in This Study

Results

Conclusions and Discussions

Acknowledgments

References

Information & Authors

Information

Published In

Copyright

History

Permissions

Authors

Affiliations

Metrics & Citations

Metrics

Citations

Download citation

Cited by

View Options

Media

Figures

Other

Tables

PREVIOUS ARTICLE

Verify Phone

Congrats!

Abstract

Introduction

Structures of Simulating Crystals

Methods

Introduction of MD Simulation

Force Fields

COMPASS Force Field

UFF

Dreiding

MD Conditions Used in This Study

Results

Conclusions and Discussions

Acknowledgments

References

Information

Published In

Copyright

History

Permissions

Authors

Affiliations

Metrics

Citations

Download citation

Cited by

Figures

Other

Share

Copy the content Link

Share with email

Share

Request Username

Create a new account

Change Password

Password Changed Successfully

Verify Phone

Congrats!