Thermal conductivity of porous jennite by molecular dynamics method

TL;DR

This study uses molecular dynamics simulations to analyze how porosity and temperature affect the thermal conductivity of porous jennite, a key component in cement paste, revealing significant decreases with porosity and complex temperature dependence.

Contribution

It provides a detailed molecular dynamics analysis of porous jennite's thermal conductivity, including porosity and temperature effects, using validated computational methods.

Findings

Thermal conductivity decreases significantly with porosity.

Temperature influences thermal conductivity, increasing then decreasing with temperature.

The empirical coherent potential model accurately predicts conductivity changes.

Abstract

The thermal conductivity of porous jennite, as the major component of cement paste, and its porosity and temperature dependences are simulated by molecular dynamics methods using ClayFF force field. The porous jennite models with different porosities are created by removing atoms within the sphere from bulk jennite model. The thermal conductivity elements of bulk jennite calculated by equilibrium Green-Kubo method are almost identical to those by non-equilibrium M\"{u}ller-Plathe (MP) method at 300 K. The volumetric thermal conductivity of porous jennite is found to decrease from 1.141 W/m$\cdot$K to 0.144 W/m$\cdot$K as increasing the porosity from 0\% to 72.22\% at 300 K, following the empirical coherent potential model when the pore is assumed to be filled with air. As increasing temperature, the thermal conductivity is observed to increase from 240 K to 560 K and gradually decrease until 1100 K for porous jennites with the porosities of 0\%, 15.99\% and 32.79\%.

Tables

Table 1. Table 1: Lattice parameters of unit cell of jennite crystal, equilibrated through N ​ P ​ T 𝑁 𝑃 𝑇 NPT run, in comparison with experiment.

	$a$ (Å)	$b$ (Å)	$c$ (Å)	$\alpha$ (∘)	$\beta$ (∘)	$\gamma$ (∘)
ClayFF	10.865	7.402	10.946	102.06	94.57	110.33
Exp [18]	10.576	7.265	10.931	101.30	96.98	109.65

Table 2. Table 2: Sphere radius ( r 𝑟 r ), total number of atoms ( N tot subscript 𝑁 tot N_{\text{tot}} ), number of deleted atoms ( N del subscript 𝑁 del N_{\text{del}} ), mass density ( ρ 𝜌 \rho ), and porosities estimated using number of atoms ( P N subscript 𝑃 𝑁 P_{N} ) and mass density ( P ρ subscript 𝑃 𝜌 P_{\rho} ) in perfect and porous jennite models.

model	$r$ (Å)	$N_{\text{tot}}$	$N_{\text{del}}$	$\rho$ (g/cm3)	$P_N$ (%)	$P_{\rho }$ (%)
1	0	2484	0	2.205	0	0
2	3	2346	138	2.099	5.56	4.79
3	5	2208	276	1.976	11.11	10.39
4	5.5	2070	414	1.852	16.67	15.99
5	6	1863	621	1.667	25.00	24.39
6	7	1656	828	1.482	33.33	32.79
7	9	1173	1311	1.050	52.78	52.39
8	11	690	1794	0.612	72.22	72.22

Table 3. Table 3: Principal ( K 1 subscript 𝐾 1 K_{1} , K 2 subscript 𝐾 2 K_{2} , K 3 subscript 𝐾 3 K_{3} ) and volumetric averaged ( K v subscript 𝐾 v K_{\text{v}} ) thermal conductivities of jennite crystal as increasing the system length, i.e. number of repeating cells N Re subscript 𝑁 Re N_{\text{Re}} .

$N_{\text{Re}}$	$K_1$ (W/m$\cdot$K)	$K_2$ (W/m$\cdot$K)	$K_3$ (W/m$\cdot$K)	$K_{\text{v}}$ (W/m$\cdot$K)
2	1.521	1.357	0.566	1.148
3	1.239	1.506	0.748	1.164
4	1.182	1.339	0.886	1.136
5	1.140	1.378	0.917	1.145
8	1.100	1.354	0.967	1.141
10	1.085	1.368	1.003	1.152

Table 4. Table 4: Principal ( K 1 subscript 𝐾 1 K_{1} , K 2 subscript 𝐾 2 K_{2} , K 3 subscript 𝐾 3 K_{3} ) and volumetric averaged ( K v subscript 𝐾 v K_{\text{v}} ) thermal conductivities of the porous jennite with various porosities ( P ρ subscript 𝑃 𝜌 P_{\rho} ) at 300 K.

$P_{\rho }$ (%)	$K_1$ (W/m$\cdot$K)	$K_2$ (W/m$\cdot$K)	$K_3$ (W/m$\cdot$K)	$K_{\text{v}}$ (W/m$\cdot$K)
0.00	1.100	1.354	0.967	1.141
4.79	0.996	1.253	0.867	1.039
10.39	0.932	1.159	0.806	0.966
15.99	0.828	0.986	0.732	0.849
24.39	0.660	0.770	0.671	0.701
32.79	0.525	0.684	0.545	0.585
52.39	0.200	0.384	0.338	0.307
72.22	0.013	0.273	0.146	0.144

Table 5. Table 5: Thermal conductivity ( K 𝐾 K ) of the cement hydration product with different water to cement ratios (w/c)

	$K$ (W/m$\cdot$K)	w/c
This work	1.030	–
Yoon et al. [31]	1.28	0.4
Xu et al. [32]	0.53	0.35
Bentz [33]	0.90 $\sim$ 1.05	0.3, 0.4
Hansen et al. [75]	0.88 (early)	0.5
	0.78 (late)
Mounanga et al. [76]	1.0 (fresh)	0.348
	1.07 (28 days)
Wu et al. [63]	$0.9568\pm 0.0257$	0.45

Table 6. Table 6: The unit cell parameters and density of jennite crystal at different temperatures, obtained by MD simulation.

	300 K	500 K	700 K	900 K	1100 K	1300 K
$a$ (Å)	10.865	11.600	12.691	13.666	12.647	14.238
$b$ (Å)	7.402	7.627	7.188	7.246	6.183	7.312
$c$ (Å)	10.946	10.474	11.014	10.832	15.031	25.811
$\alpha$ (∘)	102.06	103.66	110.64	98.35	89.83	90.78
$\beta$ (∘)	94.57	93.95	84.53	93.47	99.12	89.64
$\gamma$ (∘)	110.33	111.55	110.63	114.55	104.31	106.85
V (Å3)	796.1	824.9	879.4	956.6	1123.6	2571.8
$\rho$ (g/cm3)	2.218	2.141	2.008	1.846	1.633	0.687

Table 7. Table 7: Coefficients of the 3rd-order polynomial relation between the thermal conductivity of porous jennite and temperature.

$P$ (%)	$c_0$	$c_1(\times {10}^{-3})$	$c_2(\times {10}^{-7})$	$c_3(\times {10}^{-10})$
0	1.095	0.145	7.510	$-$11.220
15.99	0.445	2.217	$-$30.583	9.635
32.79	0.313	1.535	$-$22.750	8.029

Equations

$$E_{\text{total}}=E_{\text{Coul}}+E_{\text{vdW}}+E_{\text{bond}}+E_{\text{angle}}$$

$$E_{\text{Coul}}=\frac{e^2}{4\pi {\epsilon }_0}\sum _{i\ne j}\frac{q_i{q}_j}{r_{ij}}$$

$$E_{\text{vdW}}=\sum _{i\ne j}{D}_{ij}^0\left[{\left(\frac{R_{ij}^0}{r_{ij}}\right)}^{12}-2{\left(\frac{R_{ij}^0}{r_{ij}}\right)}^6\right]$$

$$E_{\text{bond}}=\sum _{i\ne j}{k}_1(r_{ij}-r_0{)}^2$$

$$E_{\text{angle}}=\sum _{i\ne j\ne k}{k}_2({\theta }_{ijk}-{\theta }_0{)}^2$$

$$P=\left(1-\frac{{\rho }_p}{{\rho }_0}\right)\times 100\%$$

$$K_{\alpha \beta }=\frac{V}{k_B{T}^2}{\int }_0^{\infty }⟨q_{\alpha }(t)\otimes q_{\beta }(0)⟩dt$$

$$q_{\alpha }(t)=\frac{1}{V}\frac{d}{dt}\sum _{i=1}^N{r}_{i,\alpha }{H}_i$$

$$K_{\alpha }=-\frac{q_{\alpha }}{\partial T/\partial r_{\alpha }}=-\frac{q_{\alpha }}{\Delta T/\Delta r_{\alpha }}$$

$$K^{\text{G}}=K_b\frac{2K_b+K_p-2p(K_b-K_p)}{2K_b+K_p+p(K_b-K_p)}$$

$$K^{\text{CP}}=\frac{1}{4}\bigg{\{}(3p-1)K_{p}+(2-3p)K_{b}+\\ \left[\Big{(}(3p-1)K_{p}+(2-3p)K_{b}\Big{)}^{2}+8K_{b}K_{p}\right]^{\frac{1}{2}}\bigg{\}}$$

$$K^{\text{R}}=K_b\frac{v{p}^{2/3}+1-p^{2/3}}{v(p^{2/3}-p)+1+p-p^{2/3}}$$

$$K^{\text{F}}=K_b\frac{v(1-p^{1/3}+p)+p^{1/3}-p}{v(1-p^{1/3})+p^{1/3}}$$

$$K^{\text{P}}=K_b(1-p)+K_{p}p$$

$$K^{\text{S}}=\frac{1}{(1-p)/K_b+p/K_p}$$

$$K_v^{\text{SC}}=\frac{{\sum }_{s=1}^{n_p}{f}_s{K}_v^s{B}_s^{\text{sph}}}{{\sum }_{s=1}^{n_p}{f}_s{B}_s^{\text{sph}}}$$

$$K(T)=c_0+c_{1}T+c_2{T}^2+c_3{T}^3$$

Full text

††journal: Cement and Concrete Research

1 Introduction

Building materials have always been indispensable for survival of human beings, and among them, concrete has been the most widely used since the invention of modern cement. While concrete consists of cement matrix and aggregates such as sand and gravel or industrial slag, its main binding phase is known to be calcium silicate hydrates (CaO-$${{\text{SiO}_{\vphantom{\text{}}\text{2}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}}}$$-$${{\text{H}_{\vphantom{\text{}}\text{2}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak\text{O}}}; C$-$S$-$H) gel produced by cement hydration [1, 2]. The C$-$S$-$H gel is responsible for various important properties of concrete, including strength, toughness, and durability. Therefore, considerably mass of experimental and theoretical research has been carried out in revealing and improving mechanical properties of concrete and C$-$S$-$H gels [3, 4, 5, 6, 7, 8].

Although atomic modeling is widely accepted as a powerful tool to provide invaluable insights into the fundamental behavior of materials, the modeling of real C$-$S$-$H gel is known to be exceedingly challenging. This is partly associated with the difficulties in identifying the atomistic structures of C$-$S$-$H gel by material characterization experiments, such as X-ray diffraction (XRD), transition or scanning electron microscopy (TEM or SEM) and Raman spectroscopy [9]. In fact, the C$-$S$-$H gel in concrete is imperfect crystalline, in which the calcium to silicon (C/S) ratio is changeable from 1.2 to 2.1 according to the hydration process condition [10]. In such circumstances, they are often plausibly modeled by tobermorite and jennite minerals based on their structural similarities observed by XRD [11], although more realistic model with a C/S ratio of 1.65 for C$-$S$-$H has been suggested by randomly removing charge-neutralized ${{\text{SiO}_{\vphantom{\text{}}\text{2}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}}}$ units from tobermorite 11 Å [12, 13, 14, 15, 16]. The mineral tobermorite has much smaller C/S ratio of $\sim$0.83 compared with experiment, whereas for jennite it is reasonable to be about 1.5 due to the corrugation of Ca$-$O sheets [17]. So far, jennite and tobermorite family have been intensively studied as substitutes for the C-S-H gels. Bonaccorsi et al. have done pioneering work on characterizing the crystalline structures of jennite [18], tobermorite 11 Å [19] and tobermorite 14 Å [20] through XRD experiments. The structural and mechanical properties of jennite and tobermorite have been investigated using first-principles calculations, [17, 21, 22, 23, 24, 25, 26] and molecular dynamics (MD) simulations [27, 28, 29, 30]. However, there is only a limited number of studies on thermal properties of jennite and tobermorite as well as C$-$S$-$H gel [31, 32, 33, 34, 35], in spite of its importance in design and construction of energy efficient buildings with concrete [36]. In particular, studies of revealing the thermal conductivity of jennite at molecular scale are truly scarce.

Recently, there is a growing interest in high thermal insulation materials in construction to save energy use in buildings [37, 38]. In general, thermal insulation materials are porous [39, 40]. The thermal conductivity of such porous materials decreases as their porosity increases due to an increasing amount of inside small-sized air pores. Porous concrete is in particular advantageous due to its high thermal insulation, low density, fire-proofing and acoustic insulation [41]. Qomi et al. [34] reported the thermal properties of C$-$S$-$H, such as specific-heat capacity, thermal expansion coefficient and thermal conductivity, by means of MD method using core-only CSH force field potential [42]. Jin et al. [43] proposed the fractal model for predicting the thermal conductivity of porous concrete based on their measurements of autoclaved aerated concrete with various porosities.

In this work, we aim to determine the thermal conductivity of jennite and its dependence on the porosity and temperature by using MD simulations. We organized the paper as follows. In Section 2, we described computational methods, including the employed force field potential, jennite models and MD methods for predicting the thermal conductivity. Section 3 shows the results and discussion, where comparisons with the previous MD works, experimental results, and predictions based on the continuum models for heat transfer in porous materials, were provided. Temperature dependence of thermal conductivity was also presented together with a discussion in comparison with the available literature data. Finally, we provided the conclusions of this work in Section 4.

2 Methods

2.1 Force field

The accuracy of classical MD simulations strongly depends on the selected interatomic force field potential. To date, many force fields have been developed for cementitious materials [44, 45]. They are for example ClayFF [46], IFF (interface force field) [47], CementFF [48], ReaxFF [49], and CSH-FF [42]. Among them, the ClayFF potential was selected in this work, since it has been completely tested through lots of studies on various properties of clays and cement materials [30, 50, 51, 52]. As Mutisya et al. [29] pointed out, moreover, the ClayFF is able to capture thermodynamic properties of both hydrous and anhydrous cement phases.

In the ClayFF formalism, the total energy of a molecular system is represented as follows [44, 46],

[TABLE]

where $E_{\text{Coul}}$ is the Coulomb interaction energy, $E_{\text{vdW}}$ is the short-range van der Waals (vdW) energy, $E_{\text{bond}}$ is the bond stretching term, and $E_{\text{angle}}$ is the angle bending term, respectively. The electrostatic interaction energy between the charged particles is described by the Coulomb’s law as,

[TABLE]

where $q_{i}$ is the atomic partial charge of $i$th particle, $r_{ij}$ is the distance between the $i$th and $j$th particles, $e$ is the elementary charge, and $\varepsilon_{0}$ is the dielectric permittivity of vacuum. The attractive vdW dispersion energy is written as [44, 46],

[TABLE]

where $D^{0}_{ij}=\sqrt{D^{0}_{i}D^{0}_{j}}$ and $R^{0}_{ij}=(R^{0}_{i}+R^{0}_{j})/2$. $D^{0}_{i}$, $R^{0}_{i}$ are the empirical parameters derived from the experimental data for physical property. The bonded interactions, which are used to describe the O$-$H bonding in water molecules and hydroxyl groups, and the covalent bonding in polyatomic species, are simplified as harmonic terms as follows [44, 46],

[TABLE]

where $k_{1}$ and $k_{2}$ are the harmonic force constant parameters, and $r_{0}$ and $\theta_{0}$ are the equilibrium values of bond length and bond angle respectively. The parameters $q_{i}$, $D^{0}_{i}$, $R^{0}_{i}$ and $k$ for Ca, Si, O and H atoms considered in this work were generated by using the $msi2lmp$ tool of Large Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [53], and given in Table S1$-$S3 (Supplementary data).

2.2 Models for porous jennite

The chemical formula of jennite is ${{\text{Ca}_{\vphantom{\text{}}\text{9}}^{\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{18}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak{\text{(}}\text{}\nolinebreak\text{OH}\nolinebreak{\text{)}}\text{}_{\vphantom{\text{}}\text{6}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak{\cdot}\text{}}}$8${{\text{H}_{\vphantom{\text{}}\text{2}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak\text{O}}}$, and its crystalline lattice is triclinic with $P$1̄ space group [17, 18]. Fig. 1(a) shows the unit cell of jennite (see jennite-unitcell.cif in Supplementary material). For MD simulations, we constructed a supercell for bulk jennite model by using $3\times 4\times 3$ unit cells, which contains 2484 atoms, and subsequently carried out constant pressure and temperature ($NPT$) simulation with a time step of 0.1 fs for 100 ps at 300 K and atmospheric pressure, then it’s size was about $32.6\times 29.6\times 32.8$ Å, as shown in Fig. 1(b). At this NPT simulation the cutoff radii for the LJ and Coulombic interactions were set to be 15 Å. The periodic boundary condition was applied in all three crystallographic directions, and the $pppm$ method was used for calculation of electrostatic term. Nose-Hoover thermostat and barostat were used to control temperature and pressure of the system, and the velocity Verlet time integration algorithm was used to update the position and velocity of atoms. Fig. 1(c) shows the equilibrated atomistic structure of the jennite model. Table 1 presents the determined lattice parameters, which are in good agreement with experiment with a maximal relative error of 2.73% [18].

To verify the validity of jennite model and the force field, we also evaluated the radial distribution function (RDF) of Ca$-$O in comparison with available experiment and previous MD simulations. According to the difference of local structures, calcium atoms are classified into sheet (Ca${}_{\text{s}}$) and interlayer (Ca${}_{\text{w}}$) atoms, while oxygen atoms are divided into surrounding (O${}_{\text{s}}$) and water (O${}_{\text{w}}$) oxygen atoms. Therefore, we considered four different RDFs for Ca$-$O bonds (see Fig. S1). For the Ca${}_{\text{s}}-$O bonds, the first peaks were found at 2.41 and 2.59 Å, which are slightly larger than the previous MD results of 2.37 and 2.57 Å [27] while in reasonable agreement with the experimental value of 2.45$\pm$0.06 Å [54]. The distances between Ca${}_{\text{w}}$ and O${}_{\text{w}}$ or O${}_{\text{h}}$ (hydroxyl groups) atoms were estimated to be 2.35 or 2.50 Å, being in good agreement with the XRD experimental data of 2.36 or 2.50 Å [18]. These indicate that our jennite model can be used safely and the ClayFF employed in this work is reasonably accurate and transferable for jennite.

Then, the porous jennite models were created by removing all the atoms within sphere with a certain radius located at the center of perfect jennite model. Note that if one atom of a polyatomic molecule (e.g., CaO, ${{\text{SiO}_{\vphantom{\text{}}\text{2}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}}}$, ${{\text{H}_{\vphantom{\text{}}\text{2}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak\text{O}}}$) belonged to the spherical region, the molecule was totally removed. We increased the radius gradually from 3 Å to 11 Å, and estimated the corresponding mass density and porosity of each porous model. The porosity can be estimated by using the mass densities of perfect model $\rho_{0}$ and porous model $\rho_{p}$ as follows,

[TABLE]

For reference, we also provided the porosity calculated using the number of atoms instead of density. Table 2 lists the characteristics of porous jennite models with various porosities from 0 to 72.22%. Fig. 1(d), (e) and (f) show the equilibrated atomistic structures of porous jennite models with the porosities of 10.39%, 32.79% and 52.39% respectively.

2.3 Molecular dynamics methods for predicting the thermal conductivity

The thermal conductivity K of insulating solid is defined as a second-order tensor coefficient relating the temperature gradient $\nabla T$ to the heat flux q through Fourier’s law, $\textbf{\em q}=\textbf{\em K}\cdot\nabla T$. Typically, two different methods have been developed to calculate K via MD simulation: (1) equilibrium Green-Kubo (GK) method [55, 56] and (2) non-equilibrium Müller-Plathe (MP) method [57]. The former method does not impose temperature gradients and yields the full thermal conductivity tensor via one time simulation, while the latter approach imposes a heat flux on the simulation box and gives the heat conductivity separately in each direction for an anisotropic material. In this work, we applied these two approaches in order to guarantee an accuracy of calculation, relying mostly on the MP method.

In the GK approach, once the molecular system is equilibrated through $NPT$ simulation, the net heat flux fluctuates around zero during $NVE$ run. The dissipation rate of heat flux fluctuations is evaluated, and the thermal conductivity tensor for an anisotropic materials is measured at equilibrium by using the following GK relation [34, 58],

[TABLE]

where $\otimes$ is the dyadic product in tensor notation, $k_{B}$ the Boltzmann constant, $V$ the volume of system, and $\alpha$, $\beta$ the components of coordinate system. In this equation, $\phi_{\alpha\beta}(t)=\langle q_{\alpha}(t)\otimes q_{\beta}(0)\rangle/\langle q_{\alpha}(0)\otimes q_{\beta}(0)\rangle$ is the second-order tensor known as the heat flux autocorrelation function (HFACF), and thus, the thermal conductivity can be calculated from integration of HFACF. The heat flux in a polyatomic system with a number of atom $N$ is written as follows [34],

[TABLE]

where $\textbf{\em r}_{i}$ is the position vector of the $i$th atom, and $H_{i}$ is the Hamiltonian, which is composed of kinetic and potential (Eq. (1)) energies. The GK approach has been applied to C$-$S$-$H [34] and amorphous silicon dioxide thin films [59]. In this work, the equilibrium GK approach was applied to the perfect jennite crystal with the $3\times 4\times 3$ supercell using LAMMPS [53]. After equilibration of the system through $NPT$ run, we subsequently performed $NVE$ simulation with a time step of 0.1 fs for 20 ns ($2\times 10^{5}$ correlation time steps).

In the non-equilibrium MP method, thermal conductivity can be calculated by imposing a heat flux on the simulation box and measuring the resultant temperature gradient [57, 60, 61, 62]. If the heat propagates along $\alpha$ direction in the simulation box, the thermal conductivity in this direction is calculated by

[TABLE]

As implemented in the LAMMPS [53], the simulation box is divided into an even number of parallel slices, and the temperature of each slice is evaluated at every time step and then averaged over the whole time steps [61]. The heat flux is imposed through velocity rescaling, i.e., by exchanging velocities between atoms with the largest kinetic energy of the first slice and atoms with the lowest kinetic energy of the middle slice at a certain simulation time step. After reaching steady state by imposing heat flux for sufficiently long time, a temperature profile $T(r_{\alpha})$ is evaluated and thus the thermal conductivity is ready for estimation by using Eq. (9). For an anisotropic material, this procedure has to be performed separately to obtain the thermal conductivity in each direction of coordinate system. The MP method is easy to understand and has moderate simulation time, so it has been used in many studies of porous materials [59, 61, 62].

In this work, we applied the MP method mostly to the porous jennite models with various porosities for estimating their thermal conductivities at different temperatures, using the LAMMPS program [53]. The porous jennite models, constructed by removing atoms within the spherical region as described above, were equilibrated by performing constant volume and constant temperature ($NVT$) simulations with a time step of 0.1 fs for 100 ps. Then, the simulation boxes were created by extending these equilibrated porous jennite models by from two to ten times in [100], [010] and [001] crystallographic directions, and also equilibrated by $NVT$ run. To calibrate the temperature profile and the thermal conductivity, the extended simulations boxes are divided into 20 plane-parallel slices, and micro canonical ensemble ($NVE$) simulations were performed with a time step of 0.1 fs during 200 ps. Like the application of GK approach, the Nose-Hoover thermostat and barostat, and velocity Verlet time integration algorithm were adopted.

3 Results and discussion

3.1 Thermal conductivity of jennite by Green-Kubo method

We first presented the thermal conductivity of perfect jennite calculated by GK approach using the $3\times 4\times 3$ supercell model. Figure 2 shows the HFACF along the [100], [010] and [001] directions, and thus the principal elements of thermal conductivity tensor at temperature of 300 K. After approximately 200 fs of $NVE$ simulation at equilibrium state, the HFACF was found to fluctuate around zero while the thermal conductivity to converge into a certain value. The principal values were determined to be $K_{1}=1.125$, $K_{2}=1.393$ and $K_{3}=0.921$ W/m$\cdot$K along the [100], [010] and [001] directions respectively, and the average volumetric value was estimated to be $K_{\text{v}}=1.146$ W/m$\cdot$K. As expected earlier, the thermal conductivity of jennite was confirmed to be anisotropic like $K_{2}>K_{1}>K_{3}$, which is different from $K_{1}>K_{2}>K_{3}$ for C$-$S$-$H predicted by MD [34], possibly due to differences in model and force field. The reason of such anisotropy is that fairly heat conductive Si$-$O bonds attribute to higher thermal conduction along the [010] direction (the so-called Si$-$O dreierketten chains run along the [010] direction [18]), whereas in the [001] direction water molecules and looser Ca$-$O bonds scatter phonons, resulting in lower thermal conduction in this direction [27, 34].

3.2 Influence of system length on thermal conductivity

In the previous works, it has been found that in the MP method the accuracy of thermal conductivity calculation severely depends on the system length, whereas the effect of its cross-section area is independent for the case of amorphous nanoporous silica. [61]. Therefore, we tested the influence of simulation box length on the thermal conductivity of bulk jennite to determine the proper length. In this test, the perfect jennite model was used. Due to the anisotropy of thermal conductivity, the jennite model was repeated from two to ten times in the [100], [010] and [001] directions, yielding the simulation boxes with a range of system length from $\sim$60 to $\sim$300 Å, cross-section area of about $30\times 30$ Å2, and number of atoms from 4968 to 24 840.

As depicted in Fig. 3(b), the simulation boxes were divided into 20 plane slices perpendicular to the repeating direction (heat flux propagation direction), of which the 11th slice plays the role of heat source and the first slice is for heat sink. It is worth noting that this configuration is due to the periodic boundary condition in all the three directions, while under the fixed boundary condition the first and last slices can be used for heat source and sink [62]. Temperature was fixed at 300 K through $NVT$ simulation, and subsequently $NVE$ runs were performed with a time step of 0.1 fs during a total time of 200 ps. Figures 3(c)-(e) show temperature profiles in the simulation boxes with increasing lengths along the heat flux directions of [100], [010] and [001]. It was observed that, when the system length was not enough long, the temperature reduction going from the heat source to the heat sink slices did not follow the straight line but rather meandering path, raising difficulty and ambiguity in measuring the temperature gradient. For the system lengths over $\sim$150 Å (extending supercell model five times), the temperature profile became almost the straight line.

In Table 3, we summarize the principal and volumetric averaged thermal conductivities of jennite crystal as increasing the system length. As can be seen in Fig. 4, the principal thermal conductivities converged into the reliable values over the system length of $\sim$150 Å, which is comparable with the previous value of 100 Å for amorphous silica [61]. For the sake of high reliability, we selected the appropriate system length as about 240 Å (repeating jennite model 8 times, 19 872 atoms) by considering the computational cost and statistical error.

When compared with our previous result by GK approach, the obtained thermal conductivities, $K_{1}=1.100\pm 0.025$ W/m$\cdot$K, $K_{2}=1.354\pm 0.034$ W/m$\cdot$K, $K_{3}=0.967\pm 0.020$ W/m$\cdot$K and $K_{\text{v}}=1.141\pm 0.026$ W/m$\cdot$K, are in good agreement each other. Due to a lack of experimental data for jennite, we compared our results with the previous MD works for similar materials. For C$-$S$-$H, Qomi et al. [34] obtained $0.98\pm 0.2$ W/m$\cdot$K, which is comparable with ours. Wu et al. [63] also reported the similar value of concrete hydration product as 1.015 W/m$\cdot$K through the volumetric weighted average of the alite and belite. These indicate that our models, adopted force field, and computational method are sufficiently accurate for study of porous jennite in the following subsections.

3.3 Porosity dependence of thermal conductivity

We further investigated the porous jennite with porosities from 0% to 72.22% at fixed temperature of 300 K, comparing with the different models for describing the effective thermal conductivity of porous material and with the experimental data. It should be noted that introducing pores into the perfect crystal or media reduces the thermal conductivity substantially, and thus, porous media acts as a thermal insulator or a barrier to heat removal. As mentioned above, the simulation boxes extended by 8 times in each direction were employed (see Fig. 3(a) for porous jennite with a porosity of 24.39%). However, for the case of 72.22% porosity, the number of atoms was too small (5 520 atoms) and thus the width and height of simulation box were again doubled, keeping the length unchanged, so that the number of atoms became enough large (22 080).

Before presenting the results of this work, let us consider the previously suggested models that relate the thermal conductivity of porous media to the porosity. The Garnett model was first derived to determine the electrical conductivity of a binary composite, in which spherical inclusions with random size were dispersed, as follows [64],

[TABLE]

where $K_{b}$ and $K_{p}$ are the thermal conductivities of bulk and pore, and $p$ is the porosity. Landauer [65] proposed the coherent potential (CP) model to estimate the conductivity of a binary composite made of matrix and spherical inclusions as follows,

[TABLE]

By using the CP model, Cahill and Allen [66] well predicted the thermal conductivity of porous material with a porosity of 30%. The Russel model was developed to derive the effective thermal conductivity of a dry porous material, in which cubical inclusions with uniform pores were evenly distributed [67],

[TABLE]

where $v=K_{p}/K_{b}$. The Russell model assumed a parallel heat flow, whereas the Frey model treated a series heat flow [68],

[TABLE]

The Parallel and Series models were developed to give the maximum and minimum of thermal conductivity of a porous material because planes of constituent phases are in parallel and in series with the heat flow direction, as follows,

[TABLE]

[TABLE]

Meanwhile, Qomi et al. [34] derived the self-consistent (SC) model to estimate the homogenized thermal conductivity of matrix-inclusion system, as follows,

[TABLE]

where $K_{v}^{s}$ is the volumetric thermal conductivity of spherical inclusion, $f_{s}$ is the volume fraction between the matrix and inclusion, and $B_{s}^{\text{sph}}=3K_{v}^{\text{SC}}/(2K_{v}^{\text{SC}}+K_{v}^{s})$ is the spherical localization factor of the $s$th phase.

In Table 4, we present the thermal conductivities of the porous jennite with different porosities from 0 to 72.22% calculated by applying the MP method at 300 K. Figure 5 shows normalized thermal conductivity of the porous jennite as a function of porosity, comparing with the effective thermal conductivities of porous media derived from the various models described above. As expected, the thermal conductivity of porous jennite was observed to decrease as increasing the porosity, since the thermal conductivity of air (0.026 W/m$\cdot$K [69]) is much smaller than that of jennite. As can be seen in Fig. 5, although all models exhibit decreasing tendency of effective thermal conductivity as increasing the porosity, the CP and SC models (they are almost identical each other) were found to be in quantitative agreement with our MD simulation results. It should be noted that an iterative procedure is required in the calculation with the SC model, and thus, in order to calculate the thermal conductivity of porous jennite by using the SC model we made use of the result from the Parallel model as the initial value and iterated the calculation for 10 times during the self-consistent cycle, reaching convergence after 7 or 8 times. For the case of amorphous ${{\text{SiO}_{\vphantom{\text{}}\text{2}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}}}$, the CP model has also been found to reasonably produce its thermal conductivity for a range of porosities such as 10 $\sim$ 35% [61] and 0.54 $\sim$ 24.15% [62]. Furthermore, Santos [70] revealed that the thermal conductivity of conventional aluminous refractory concrete exhibited the decreasing tendency as increasing the porosity from 0 to $\sim$36.37% at various temperatures, nearly according to that from the CP model. Meanwhile, Coquil et al. [61] found a slight overestimation of the CP model compared with the MD result, whereas Zhu et al. [62] found an underestimation when the porosity was 24.15%. But in our study the MD results were in good agreement with those from the CP model over the porosity range from 0 to 52.39% and not ruled out for the porosity of 72.22%. Based on this analysis, it can be concluded that the CP model is able to well predict the thermal conductivity of not only amorphous porous silica but also porous cement based materials.

Next, we made a careful analysis of our MD result in comparison with those from the CP model, in which the pores were assumed to be filled with air or water, and from experiment. Whether the pores of porous cement hydration product would be filled with air or water (or even vacuum) depends on its hardening process condition. If the pores are assumed to be filled with water, the porous jennite has the maximal thermal conductivity, whereas filling with air will give the minimal values, as obtained by the CP model shown in Fig. 6. It is worth noting that for the case of water filling the thermal conductivity of porous jennite should decrease with increasing the porosity as well, due to smaller conductivity of water (0.6062 W/m$\cdot$K [69]) than jennite. According to the work of Constantinides and Ulm [74], the porosities of high and low density C$-$S$-$H were found to be 24 and 36% respectively, and therefore, if hydrated cement paste would be composed of only jennite, its thermal conductivity should be ranged from 0.553 W/m$\cdot$K with air filling to 0.992 W/m$\cdot$K with water filling on the CP model lines.

In our previous work [77], we demonstrated that when applying the multi-functional composite coating material [78] made from silica/alumina-rich minerals to the concrete wall, high dense ceramic layer with a thickness of $\sim$1 mm and a porosity of 10.05% can be formed on the wall surface by two-step hydration of cement and mineral powders. By using the laser flash method, we measured the thermal conductivity of this ceramic layer to be 1.030 W/m$\cdot$K at 223 K (see Table S2), which is placed between the red and blue lines as shown in Fig. 6. In Table 5, we summarize the thermal conductivity of cement hydration products with different water to cement ratios (w/c), measured from experiment or calculated from theory. Wu et al. [63] calculated the thermal conductivity of hardened cement paste to be $0.96\pm 0.03$ W/m$\cdot$K when its hydration degree was 0.95 and w/c =0.45. When compared with the available experimental data of cement paste with different w/c ratios such as $0.90\sim 1.05$ W/m$\cdot$K with the hydration degree of about $0.3\sim 0.8$, 1.28 W/m$\cdot$K [31], 0.78 and 0.88 W/m$\cdot$K for early and late stage of hydration process [75], 0.53 W/m$\cdot$K [32], and 1.0 and 1.07 W/m$\cdot$K for fresh and after 28 days [76], our predicted values with MD simulation can be said to be well suited with them.

3.4 Temperature dependence of thermal conductivity

At the final stage, we investigated the thermal conductivity of porous jennite as increasing temperature from 240 K to 560 K with a step of 60 K for details and from 300 K to 1300 K with a step of 200 K. For the sake of simplicity without losing generality, we selected the three different typical porosities of 0%, 15.99% and 32.79%. We followed the above-mentioned procedure to calibrate the thermal conductivity: doing $NPT$ simulation of perfect $3\times 4\times 3$ jennite model for 100 ps, creating the pore with a radius corresponding to the porosity and running $NVT$ for 400 ps, constructing simulation box by repeating 8 times the model and running $NVT$ for 100 ps, and finally performing $NVE$ run for 200 ps with imposing the heat flux, all these with a time step of 0.1 fs. These processes were repeated at different temperatures.

Table 6 lists the unit cell parameters and density of perfect jennite crystal obtained by doing $NPT$ simulations as increasing temperature from 300 K to 1300 K. It was observed that as increasing temperature the unit cell volume increases monotonically while the mass density decreases, although the lattice constants varies irregularly. It should be noted that at 1300 K the density was too small to think of crystal; it was actually transformed to liquid state. Moreover, as revealed by Yu and Kirkpatrick [79], the jennite crystal is known to be transformed to wollastonite and larnite above 1073 K. That is why we did not perform further simulation at 1300 K. Another point is that since the volume of jennite model increases with increasing temperature, the radius of spherical region for creating the pore should be adjusted to meet the target porosity.

Figure 7 shows the temperature profiles in the simulation boxes of perfect jennite along the principal crystallographic directions. At moderate temperatures, one can see a good linearity of temperature change along the slices, which makes it easy to measure the temperature gradients, while an inappreciable dispersion of temperatures was found at 1100 K.

In Fig. 8, we show the temperature dependence of thermal conductivities of porous jennite with porosities of 0, 15.99 and 32.79%. For the three cases of porosities, the principal thermal conductivities, especially volumetric averaged conductivities, increase to 500 K and decrease after that. Also, it was found that the [100] element of thermal conductivity became the biggest after 900 K in all the three cases. Figure 8(d) shows the volumetric thermal conductivities calculated in the detailed temperature range from 240 K to 560 K as well as in the temperature range from 300 K to 1100 K in comparison with our experimental values measured for the ceramic coating layer on the concrete wall by using the laser flash method. We note that our multifunctional composite building material was developed for use of thermal insulation, moisture barrier, sterilization of buildings [78]. Our experimental values were found to mostly locate between the porosity of 15.99% and 32.79%. Since the porosity of this ceramic coating layer was measured to be about 10% [77], these are slightly deviated from the predicted values, possibly due to not a pure jennite. Based on our calculation result, we suggested the 3rd-order polynomial relation between the thermal conductivity of porous jennite and temperature by

[TABLE]

where $c_{i}$ is the coefficient. Table 7 shows them resulting from interpolation process and in Fig. 8(d) we see that the 3rd-order interpolation fits well with our calculation data.

We compared our calculation results with other experimental data for cement paste and concrete. Ukrainczyk et al. [80] conducted experiments with calcium aluminate cement with w/c = 0.3, 0.4, finding that its thermal conductivity increases with increasing temperature from 293 K to 353 K due to loss of water from jennite at $343\sim 363$ K. At this temperature, the lattice constant $c$ of jennite crystal was found to shrink and simultaneously the jennite transformed from {{\text{Ca}_{\vphantom{\text{}}\text{9}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak{\text{[}}\text{}\nolinebreak\text{Si}_{\vphantom{\text{}}\text{6}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak\text{O}_{\vphantom{\text{}}\text{18}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak{\text{(}}\text{}\nolinebreak\text{OH}\nolinebreak{\text{)}}\text{}_{\vphantom{\text{}}\text{6}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak{\text{]}}\text{}}}$$\cdot8${{\text{H}_{\vphantom{\text{}}\text{2}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak\text{O}}}$ to {{\text{Ca}_{\vphantom{\text{}}\text{9}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak{\text{[}}\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{\text{]}}\text{}\nolinebreak{\text{(}}\text{}\nolinebreak\text{OH}\nolinebreak{\text{)}}\text{}_{\vphantom{\text{}}\text{8}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}}}$$\cdot2${{\text{H}_{\vphantom{\text{}}\text{2}}^{\vphantom{\text{}}\vphantom{\text{\smash[t]{2+}}}\text{}}\nolinebreak\text{O}}}$, called metajennite [79]. Xianzhi et al. [81] demonstrated that the thermal conductivity of cement paste decreased like 0.82, 0.78, 0.60 W/m$\cdot$K at 333, 363 and 393 K. Wang et al. [82] reported that the fly ash concrete had decreasing thermal conductivity from 1.69 to 0.95 W/m$\cdot$K as increasing temperature from 293 K to 793 K. As they pointed out, the change of microstructure was found to start at 573 K from SEM image observation, and Yu [79] stated that jennite lost the interlayer water molecules and thus basic layer structure upon heating to 623 K. It can be said that our calculation results are similar to these experimental data. On the contrary, some irregular changing tendency was observed from experiment for different kinds of concrete. For example, Gencel et al. [73] found that the thermal conductivity of concrete, made using vermiculite as dominant component, decreased from 293 K to 1193 K and again increased to 1373 K. Khaliq and Kodur [83] studied the temperature effect on the thermal properties for self-consolidating concrete, finding that its thermal conductivity decreased up to 693 K, marginally increased to 793 K, and again decreased to 1093 K. For aluminous refractory concrete, Santos [70] reported the increasing tendency from 300 K to 333 K, significantly decreasing to 793 K, and again gradual increasing to 1273 K. For explanation of such complicated behavior, the moisture loss by evaporation of free and pore water with the increase of temperature was considered as the main cause, which is different from our calculation models with the pore filled with air.

4 Conclusions

In this work, we have investigated the thermal conductivity of porous jennite by means of molecular dynamics method, aiming at clarifying the thermal behavior of cement hydration product C$-$S$-$H gel. The porous jennite models with different porosities have been created by removing atoms within spherical region with a radius corresponding to the porosity, using $3\times 4\times 3$ supercell model for perfect jennite. We have employed two different MD methods, namely equilibrium Green-Kubo and non-equilibrium Müller-Plathe methods, to calculate the thermal conductivity of jennite along the crystallographic [100], [010] and [001] directions. We found that the thermal conductivity of perfect jennite calculated by GK method using the supercell model were almost identical to those by MP method, exhibiting the crystallographic anisotropy like $K_{2}>K_{1}>K_{3}$. The non-equilibrium MP method has been mostly applied in this work to investigate the porosity and temperature dependences of thermal conductivity of porous jennite, using the simulation boxes extended the porous jennite model 8 times along the heat flux direction and divided into 20 plane slices. At room temperature of 300 K, the principal and volumetric thermal conductivities of bulk jennite were determined to be $K_{1}=1.100\pm 0.025$, $K_{2}=1.354\pm 0.034$, $K_{3}=0.967\pm 0.020$ and $K_{\text{v}}=1.141\pm 0.026$ W/m$\cdot$K in reasonable agreement with the experimental values for C$-$S$-$H. We have found the decreasing tendency of effective thermal conductivity of porous jennite as increasing the porosity, which was confirmed to follow the empirical coherent potential or self-consistent model proposed for porous materials among other kinds of models. Finally we have investigated the temperature dependence of thermal conductivity of porous jennite with the porosities of 0, 15.99 and 32.79%, finding its gradual increase from 240 K to 560 K and decrease from 560 K to 1100 K. By comparing with our own measurements and available experimental data for cement paste and concrete, the calculated values were said to be high reasonable and reliable. Therefore, we believe this work will be helpful for a design of high thermal insulating materials with cement and concrete for energy efficient buildings.

Acknowledgments

This work belongs to the state research project “Design of Innovative Functional Materials for Energy and Environmental Application” (No. 2016-20), supported partially by the State Committee of Science and Technology, Democratic People’s Republic of Korea. The simulations have been performed on the HP Blade System C7000 (HP BL460c) that is managed by Faculty of Materials Science, Kim Il Sung University.

Notes

Declarations of interest: none.