Molecular dynamics study of the effect of moisture and porosity on thermal conductivity of tobermorite 14 Å

The structural model of tobermorite 14 Å was based on the ordered layer structure of monoclinic phase with a space group of B11b and the chemical formula of ${{\text{Ca}_{\vphantom{\text{}}\text{5}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak\text{Si}_{\vphantom{\text{}}\text{6}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak\text{O}_{\vphantom{\text{}}\text{16}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak{\text{(}}\text{}\nolinebreak\text{OH}\nolinebreak{\text{)}}\text{}_{\vphantom{\text{}}\text{2}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak{\cdot}\text{}}}$7${{\text{H}_{\vphantom{\text{}}\text{2}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak\text{O}}}$ [17]. Although calcium atoms and water molecules in interlayer space can occupy the sites randomly in natural sample of tobermorite 14 Å, they alternately reside in the locations in the case of such ordered model. Figure 1 shows the unit cell of tobermorite 14 Å.

To make modeling of porous tobermorite 14 Å, the monoclinic $9\times 9\times 2$ supercell with dimensions 61.020 $\times$ 65.835 $\times$ 56.968 Å was created, and then, the spherical voids with diameters of 25 Å, 32.7 Å and 40 Å were generated by removing all the atoms in the spherical regions. In this process, if a certain atom of a molecule was included in the void region, all other atoms in the molecule were also removed. Then the dimensions of voids were slightly enlarged and their shapes became tangled as shown in Fig. 2. By counting the number of removed atoms, the porosities corresponding to these porous models were estimated to be 24.7, 35.8 and 50.0% respectively. According to the classification of pores [16], these can be thought as gel pores and even capillary pores. It is worth noting that according to the previous studies for porous materials the porosity and temperature are the major factors influencing on their thermal conductivities, whereas shape and arrangement of pores are not significant factors [34, 35, 36]. Therefore, we expect that the porous models with the tangled spherical voids can be used safely to study thermal conductivity of porous cement paste.

We went further into the porous and moist models, which were made by inserting a certain number of water molecules into the pores. In fact, when the pores inside cement paste and concrete are vacuum or filled with air at dry condition, the concrete wall as a whole can insulate the building effectively against heat conduction due to low thermal conductivity of air (0.026 W/m$\cdot$K [37]). In reality, however, the inside pores of cement paste and concrete should be filled with water, which can occur during hydration and age or operation of buildings, resulting in dramatic increase of thermal conductivity of concrete wall due to relatively high thermal conductivity of water (0.606 W/m$\cdot$K [37]). We gradually increased the number of confined water molecules inside the pores, and then, the maximum numbers, which were from equilibrated big water simulation box at 300 K, were determined to be 1574, 2282 and 3187 for the pores with porosities of 24.7, 35.8 and 50.0% respectively. Figure 3) shows the models with 500 and 3187 water molecules in the case of 50% porosity.

It is necessary to define the moisture content of moist samples. The moisture content can be defined by mass as [8],

where the $m_{w}$ and $m_{d}$ are the masses of moist and dry samples respectively. Also, according to EN ISO 10456 [38], we used the definition of moisture content by volume as [14, 15],

where $\rho_{\text{wat}}$ is the water density, being 1000 kg/m³, and $V$ is the volume of the sample.

In the classical MD simulations, it is critical for accuracy to choose the proper interatomic force field for the subject under study. There have been developed lots of force field (FF) potentials for cementitious materials, including ClayFF, IFF (interface force field), CementFF, ReaxFF and CSH-FF [39, 40, 41, 42, 43, 44, 45]. In this work, we choose the ClayFF potential, which has been the most widely used to study cementitious materials over the years, providing sufficiently reliable results in comparison with experiments [27, 29, 32, 46].

The ClayFF potential is based on the single point charge (SPC) water model to represent the water, hydroxyl and oxygen$-$oxygen interactions. In this model, the bond stretching energy is represented by the following harmonic potential,

where $r_{ij}$ is the distance between the $i$th and $j$th atoms, and $k_{1}$ and $r_{0}$ are the stretching force constant and equilibrium bond length respectively. The bending energy of bond angles is written as follows,

where $\theta_{ijk}$ is the bond angle between the three atoms, and $k_{2}$ and $\theta_{0}$ are the bending force constant and equilibrium bond angle respectively. The ClayFF potential also includes the Coulombic and van der Waals (vdW) dispersion interactions as follows,

where $q_{i}$ is the partial charge of the $i$th particle, $e$ the elementary charge, and $\varepsilon_{0}$ the dielectric permittivity of vacuum. And $D^{0}_{ij}=\sqrt{D^{0}_{i}D^{0}_{j}}$ and $R^{0}_{ij}=(R^{0}_{i}+R^{0}_{j})/2$, where $D^{0}_{i}$ and $R^{0}_{i}$ are the empirical parameters derived from the experimental data for structural and energetic properties. Then, the total energy of the system can be written as the sum of the above force field potentials [39, 44],

We generated the force field parameters for Ca, Si, O and H atoms by using msi2lmp tool of LAMMPS (Large Scale Atomic/Molecular Massively Parallel Simulator) [47].

The thermal conductivity K of material, defined as a second-order tensor coefficient relating the temperature gradient to the heat flux through the Fourier’s law, $\textbf{\em q}=\textbf{\em K}\cdot\nabla T$, can be calculated via MD simulations using two different methods: (1) equilibrium Green-Kubo (GK) method [48, 49] and (2) non-equilibrium MD (NEMD) or Müller-Plathe (MP) method [50]. In the GK method, which is based on the fluctuation-dissipation theory at an equilibrium state of system, the thermal conductivity can be calculated by integrating the heat flux autocorrelation function (HFACF) obtained by sufficiently long time run. Although the thermal conductivity tensor can be obtained through one time simulation, the GK approach has not been applied to the porous materials in the previous works. Accordingly we also used the GK method for the bulk tobermorite 14 Å without any void in this work.

Meanwhile, the NEMD method has been widely used to predict the thermal conductivity of materials with and without porosity, specially porous materials [34, 36, 50, 51]. In this work, we also apply the NEMD method mainly to the porous and moist tobermorite 14 Å, through at least three simulations to determine the anisotropic thermal conductivity in each direction. In the NEMD method, the simulation box is divided into an even number of parallel slabs along the measuring direction, and the energy exchange between the 1st and middle slabs is imposed, i.e., the velocity of the atom having the lowest kinetic energy of the middle slab is exchanged with the velocity of the atom having the largest kinetic energy of the 1st slab at every proper time step. At this time, the temperature of the middle slab increases, whereas the temperature of the 1st slab decreases. Here, the temperature of the $i$th slab is calculated as,

where $n_{i}$ is the number of atoms in the $i$th slab, $k_{B}$ the Boltzmann constant, and $m_{\alpha}$ and $v_{\alpha}$ the mass and velocity of the $\alpha$th atom. After sufficiently long simulation time, the temperatures of all the slabs fluctuate around certain values, and the temperature gradient along the slabs should becomes linear. The periodic boundary condition (PBC) is in general applied in the three crystallographic directions, and then the thermal conductivity is calculated through the Fourier’s law as follows,

where $J$ is the heat flux that can be obtained from simulation, $\Delta T/\Delta x$ the temperature gradient along the heat flux direction, $A$ the cross sectional area through which the heat flux passes, and the number 2 is due to PBC.

In the NEMD approach, the length of simulation box is a critical factor for accuracy; for the cases of jennite and amorphous silica in the previous works the lengths were set to be larger than 110 Å [36] respectively. In accordance with this, we repeated the $9\times 9\times 2$ supercell model three times (about 170 Å) along the [100], [010] and [001] directions respectively for the porous models. For the case of perfect tobermorite 14 Å the supercell was repeated by four times so that the length became 230 Å. We set the number of dividing slabs to be 20 in this work. The energy exchange was conducted at every 100 step during the simulation of $10^{6}$ steps and further $5\times 10^{5}$ steps for averaging the temperature difference.

All the simulations in this work were performed by using the LAMMPS code. The timestep was set to be 0.1 fs, and the cutoff radii for force calculation to be 15 Å and 10 Å in molecular mechanics (MM) or energy minimization and dynamics respectively. The $pppm$ solver was used to evaluate the long-range Coulombic interaction.

We assessed the transferability of the ClayFF interatomic potential to tobermorite 14 Å through calculation of its lattice parameters, radial distribution functions (RDF) and elastic properties in comparison with the previous experimental and simulation data. In order to determine the lattice parameters of bulk tobermorite 14 Å, we conducted the energy minimization of the unit cell by using the conjugate gradient algorithm and the $pppm$ solver, where the lattice constants $a$, $b$, $c$ were controlled independently by external pressure. As can be seen in Table 1, the lattice parameters were determined to be in good agreement with the experimental values [17] with the maximum relative error of 1.8%.

Figure 4 shows RDFs for Si$-$O and Ca$-$O bonds, calculated by using the $4\times 4\times 1$ supercell (1664 atoms). In this figure, one can see that the distance between Si and O atoms is 1.58$\pm$0.04 Å and that between the intralayer Ca and O atoms is 2.48$\pm$0.15 Å, being agreed well with the experimental data of 1.56$\sim$1.74 Å and 2.27$\sim$2.74 Å [20, 17]. Note that for the case of tobermorite 11 Å, the intralayer Ca$-$O length is estimated to be shorter as 2.43$\pm$0.10 Å [52]. On the other hand, the distance between the interlayer Ca${}_{\text{w}}$ and O${}_{\text{w}}$ atoms was calculated to be 2.39$\pm$0.06 Å in our simulation (Fig. 4(c)), which is also in good agreement with the XRD experimental data of 2.36 or 2.55 Å [17].

We also calculated the elastic properties of bulk tobermorite 14 Å such as bulk modulus $B$, shear modulus $G$, Young’s modulus $E$ and Poisson ratio $\nu$. Table 2 presents the calculated values in comparison with the previous classical MM and ab initio simulations. Our calculated values are found to be in reasonable agreement with other classical MM simulations using the COMPASS force field conducted by Hajilar et al. [26] and Al-Ostaz et al. [24] and the Buckingham potential by Manzano et al. [23] and ab initio simulation by Dharmawardhana et al. [53]. When compared with the available experimental value for the bulk modulus of 47$\pm$3 GPa [54], our simulated value was slightly overestimated. These indicate that the ClayFF potential employed in this work is reasonably accurate and transferable for tobermorite 14 Å.

As a preliminary work, we calculated the thermal conductivity of bulk tobermorite 14 Å by using both the GK equilibrium and MP NEMD methods. Comparison with the available previous data was provided. For applying the GK method, we created the $4\times 4\times 1$ supercell (1664 atoms) as this method is less size dependent, and performed $NVE$ simulation for 10 ns ($10^{8}$ steps) with a time step of 0.1 fs. Figure 5 shows HFACFs and the thermal conductivities determined from those along the [100], [010] and [001] directions at 300 K. After the initial transition period of $\sim$200 fs, the system was found to reach the equilibrium state as showing the fluctuation of HFACF around zero point, and the thermal conductivities converged to a certain value. The principal elements of thermal conductivity tensor were calculated to be $K_{1}=1.704\pm 0.281$ W/m$\cdot$K, $K_{2}=2.092\pm 0.325$ W/m$\cdot$K and $K_{3}=1.149\pm 0.245$ W/m$\cdot$K in the [100], [010] and [001] directions respectively, giving the averaged volumetric value of $K_{\text{v}}=1.648\pm 0.283$ W/m$\cdot$K. These values are clearly bigger than those of tobermorite 11 Å of 0.98$\pm$0.2 W/m$\cdot$K, obtained from CSH FF potential GK simulation by Qomi et al. [7]

As for the anisotropic aspect, the relation of $K_{2}>K_{1}>K_{3}$ was found like the jennite model. The biggest value of $K_{2}$ is probably due to fairly heat conductive Si$-$O bonds that are aligned along the [010] direction [18], while the smallest value of $K_{3}$ is due to the interlayer water molecules and looser Ca$-$O bonds, which scatter phonons in the [001] direction [7, 18, 28]. In contrast to this, Qomi et al. [7] found the different relation of $K_{1}>K_{2}$ due to their different model, which contain many defects in the silicate chains along the [010] direction [56].

For another check, we also performed the NEMD simulations to compute the thermal conductivity of bulk phase, using the simulations boxes obtained by repeating the $9\times 9\times 2$ supercell model four times in the [100], [010] and [001] directions (about 230 Å and 67 392 atoms). The principal values of thermal conductivity were determined to be $K_{1}=1.596\pm 0.184$, $K_{2}=2.159\pm 0.25$ and $K_{3}=0.921\pm 0.232$ W/m$\cdot$K and the averaged value of $K_{\text{v}}=1.559\pm 0.222$ W/m$\cdot$K. These are more or less agreed with those by the GK method, and moreover, the same anisotropic feature of $K_{2}>K_{1}>K_{3}$ was observed.

From the calculated thermal conductivity of bulk tobermorite 14 Å, the thermal conductivity of cement paste can be estimated by using the upscaling model. In fact, cement paste is known to be a multiphase porous composite, consisted of C$-$S$-$H ($50\sim 70\%$ of the total volume of cement paste), CH ($20\sim 25\%$), ettringite ($15\sim 20\%$) and hydrogarnet (about 2%). Therefore, the thermal conductivity of cement paste should be calculated by averaging the values of these phases. For instance, Qomi et al. [7] evaluated the mean free path of phonons in various cement clinkers and C$-$S$-$H by applying the kinetic theory of heat transport, finding that these lengths were on the order of Ca$-$O and Si$-$O bond lengths. Therefore, it can be concluded that the heat transport in C$-$S$-$H and CH is in the vibration-dominated regime and the mean-field homogenization theory can be used in upscaling model for cement paste.

Hence we applied the self-consistent (SC) model derived from the mean-field homogenization theory [7] to evaluate the thermal conductivity of porous tobermorite 14 Å and cement paste by using the following equation,

where $n_{p}$ is the number of phases, and $K_{v}^{s}$, $f_{s}$ and $B_{s}^{\text{sph}}=3K_{v}^{\text{SC}}/(2K_{v}^{\text{SC}}+K_{v}^{s})$ are the volumetric thermal conductivity, the volume fraction and the spherical localization factor of the $s$th phase. The volumetric thermal conductivity of bulk tobermorite 14 Å was treated by using the average of the values obtained by the two different approaches in this work. Considering the various porosities of low-density (porosity, 36%) and high-density (porosity, 24%) and various wet conditions of dry and fully saturation by water, the thermal conductivity of C$-$S$-$H paste was estimated to be 0.766$\sim$1.312 W/m$\cdot$K by using Eq. 10. Using the previously calculated thermal conductivity of CH as 1.32 W/m$\cdot$K [7] and assuming that the cement paste is composed of 70% C$-$S$-$H and 30% CH, the thermal conductivity of cement paste was determined to be 0.911$\sim$1.314 W/m$\cdot$K. Table 3 lists the experimentally determined thermal conductivity of cement paste, indicating that our calculated value is in the range of experimental values [4, 6, 11].

As described above, the porous samples with different porosities of 24.7, 35.8 and 50% (see Fig. 2) were created using the $9\times 9\times 2$ supercell, and gradually increasing numbers of water molecules were added into the pores (see Fig. 3).

These samples with different porosities and water contents were equilibrated through $NVT$ runs at 300 K for $10^{6}$ steps. Then the extended samples in each crystallographic direction were generated by repeating the equilibrated samples 3 times along the desired direction, and again equilibrated through $NVT$ runs with $10^{6}$ steps. Finally $NVE$ simulations were performed for $10^{6}$ and $5\times 10^{5}$ steps to measure their thermal conductivities with the MP approach. The above process was performed 3 times for each direction of [100], [010] and [001]. Figure 6 typically shows the simulation box extended in [010] direction, with 2000 water molecules and the porosity of 50%, and the temperature profile along the slabs at the equilibrium state by the kinetic energy exchange. In the figures, water molecules were observed to be accumulated in the vicinity of pore surface and thus remain voids with small sizes, and the temperature gradient was seen to be nearly linear.

The volumetric thermal conductivities of porous and moist samples were obtained by averaging the three different values in the three different directions. Figure 7 shows the volumetric thermal conductivity of porous and moist tobermorite 14 Å as increasing the porosity and water content. It was found that the thermal conductivity decreases as increasing the porosity from 24.7% to 35.8% and to 50% while increases as a linear function of moisture content by mass (Eq. 1). As shown in Fig. 7, the linear coefficients derived from regression lines were found to decrease as increasing the porosity in agreement with the experimental observation [8].

In Fig. 8, we show the normalized thermal conductivity obtained versus dry tobermorite 14 Å ($K_{0}$) as a function of volume moisture content as suggested by standards [38, 57]. Unlike the thermal conductivity, the normalized ones were found to fit into the exponential functions with correlation coefficients of over 0.996. The exponents in the lines were obtained to be 1.583, 1.733 and 1.791 for the porosities of 24.7, 35.8 and 50% respectively. When compared with the experimental result 4.24 by Gomes et al. [15] and the standard value 4.00 for mortar with mass densities of 250$\sim$2000 kg/m³ [38, 57], our simulated exponents are clearly smaller. This is probably due to the different density (or porosity) and different components of mortar; in our case only tobermorite 14 Å as C$-$S$-$H paste was considered but other components of mortar such as sand and CH were not included. The density and porosity are also different; in the experimental works the densities of cement concretes were 415$\sim$630 kg/m³ (porosities of 76$\sim$84%) [8] and 430$\sim$830 kg/m³ [15], which are lighter than our density values of 1104$\sim$1662 kg/m³ (porosities of 24.7$\sim$50 %). Here it should be emphasized that the thermal conductivity of porous and moist tobermorite 14 Å increases exponentially when the moisture content increases and the exponent increases with the increase of porosity.

We also analyzed the data as an explicit function of porosity. Figure 9 shows the thermal conductivities as a function of porosity, together with the analysis lines predicted by the SC model for dry and fully wet porous tobermorite 14 Å. At a given porosity the thermal conductivities should lay in the vertical line to the porosity axis between the lines. In this figure, one can see that the thermal conductivities of dry porous tobermorite 14 Å decreases properly according to the SC model line (i.e., the minimum values at the porosities of 24.7, 35.8 and 50% are seen nearby the dash line). When the pores were fully saturated by water molecules, however, the simulated thermal conductivities were found to be located over the SC model line (i.e., the maximum values are observed above the black solid line), probably due to the confined water molecules in the pores. Although the thermal conductivity of bulk water was reported to be 0.6062 W/m$\cdot$K [37], those of confined water in the nano-size pore should be remarkably changed as their atomistic structure and thermodynamics vary in the C$-$S$-$H gel [7, 30, 58, 59]. In fact, the trapped water within the C$-$S$-$H gel was reported to exhibit the features of supercooled water [7] and the diffusion coefficient of the pore water was found to significantly reduce due to that the stable hydrogen bonds connected with Ca$-$OH and Si$-$OH groups restrict the transport of the surface water molecules [30]. At the position of 10$\sim$15 Å away from the surface of pore, however, the confined water molecules can be expected to behave bulk water. According to the quasielastic neutron scattering experiment [58], the confined water in C$-$S$-$H can be distinguished into chemically bound water, physically bound water and unbound water. The chemically or physically bound water was regarded to interact with the surface atom of the gel pore, whereas the unbound water to be confined within the capillary pore [58]. The diffusion coefficients of gel pore water were measured to be $D\approx 10^{-10}$ m²/s and those of capillary pore water to be $D\approx 10^{-9}$ m²/s [58], which are significantly smaller than those of bulk water, due to the interaction of the water molecules with C$-$S$-$H. Based on this consideration, we can expect that the thermal conductivity of confined nano-pore water in C$-$S$-$H is larger than that of bulk water.

To prove this expectation, we calculated the thermal conductivity of water confined in the nano-pore with a diameter of 25 Å (porosity 24.7%) by applying the GK approach. The thermal conductivity of confined water converged to 0.9257 W/m$\cdot$K, and by using this value and applying the SC model, the analysis line was obtained as shown by orange line in Fig. 9. Then the maximum thermal conductivities of porous and moist tobermorite 14 Å lay on the orange line. As increasing the diameter of pore, the thermal conductivity of confined water will decrease, approaching to the value of bulk water.

We developed a new kind of multifunctional composite building materials that has functions of good thermal insulation and moisture proofing for concrete buildings [60]. As possible mechanism for such multifunctions, two-step hydration of cement paste and mineral powder (SiO₂ and Al₂O₃) with particle sizes of 0.5$\sim$50 $\mu$m has been proposed to occur, forming the ceramics layer (porosity 10%) on the surface of concrete wall [61]. This ceramic layer can prohibit the percolation of water molecules from outside environment, which inhibits the concrete wall effectively from increasing the thermal conductivity and facilitate the thermal performance of concrete wall.

We calculated the phonon densities of states (PDOS) of various tobermorite 14 Å models since PDOS is related to variety of thermodynamic properties such as thermal conductivity and specific heat capacity. In general, heat conduction in materials can occur by different heat carriers, which include electrons, phonons and photons. Among such heat carriers, the most contributer can be different according to the material and temperature. For dielectric solids, the major heat conduction is through the phonon, of which vibration is harmonic at low temperature. The PDOS can be calculated with two different methods: the eigenvalue decomposition of the dynamic matrix and the Fourier transform method of the velocity autocorrelation function (VACF) [7, 27]. In this work we used the second method via Fourier transform of VACF, in which the PDOS $g(\omega)$ can be computed by the following equation,

where $\omega$ is the frequency, $v_{\alpha}$(t) the velocity of the $\alpha$th atom, and $\langle v_{\alpha}(t)\cdot v_{\alpha}(0)\rangle/\langle v_{\alpha}(0)\cdot v_{\alpha}(0)\rangle$ is the VACF.

Figure 10 shows the PDOS of bulk tobermorite 14 Å, porous sample, and porous and moist sample with 500 and 1000 water molecules in the case of 50% porosity. In this figure, the VACFs were observed to fluctuate around zero after about 300 fs of simulation time. To clarify the individual contributions from Ca, Si, O atoms, OH group, interlayer water and confined water molecules, their partial PDOSs were also analyzed.

In all the partial PDOSs, the range of Ca$-$O bond vibration was found to be $0\sim 430$ cm^-1, which is in good agreement with the experimental data of $0\sim 400$ cm^-1 [62] and the simulation data of $0\sim 450$ cm^-1 [7, 27, 63]. Meanwhile, the range of Si$-$O bond vibration was found to be $0\sim 1200$ cm^-1 in good agreement with the experiment [62] and the previous simulations [7, 27]. Here the phonon vibrations with $0\sim 400$ cm^-1 were found to be associated with the coupled rotation of ${{\text{SiO}_{\vphantom{\text{}}\text{4}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}}}$ [7]. One can see the peaks around 420 cm^-1 and 630 cm^-1, which are close to the experimentally observed peaks around $400\sim 500$ cm^-1 and 660 cm^-1, being associated with the deformation of ${{\text{SiO}_{\vphantom{\text{}}\text{4}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}}}$ and bending of Si$-$O$-$Si bond angles [62]. The phonon vibration with $800\sim 1200$ cm^-1 was reported to be associated with symmetric and asymmetric stretching of Si$-$O bonds. For the OH partial PDOS, the phonon of OH species was found to be smaller than that of Si atoms, which is similar to the previous ab initio work [27] but slightly different from the MM work [7] due to different C$-$S$-$H models. The phonon vibrations with 1610 cm^-1 were found to be associated with bending of the H$-$O$-$H angle of water molecule, being agreed well with the previous works [7, 27, 62]. The phonon frequencies of confined water species were found around 1595 cm^-1. On the other hand, the phonon vibrations found over 3600 cm^-1 correspond to the symmetric and asymmetric stretching of O$-$H bonds of water molecules and the stretching of hydroxyl species in the intralayer. In the porous and moist samples, the peaks of interlayer water, confined water molecules and hydroxyl groups were found at 3660, 3730 and 3740 cm^-1 respectively.