# Multi-Scale Multi-Physics Multi-Technique modeling for establishing composition-(micro)structure-property correlations in cement-based materials

Tulio Honorio

► To cite this version:

Tulio Honorio. Multi-Scale Multi-Physics Multi-Technique modeling for establishing composition-(micro)structure-property correlations in cement-based materials. Mechanics of materials [physics.class-ph]. Université Paris - Saclay, 2023. tel-04051504

HAL Id: tel-04051504

https://hal.science/tel-04051504

Submitted on 30 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Multi-Scale Multi-Physics Multi-Technique modeling for establishing composition-(micro)structure-property correlations in cement-based materials

## Modélisation Multi-Échelle Multi-Physique Multi-Technique pour la corrélation composition-(micro)structure-propriétés des matériaux cimentaires

Habilitation à diriger des recherches
de l'Université Paris-Saclay

présentée et soutenue à Gif-sur-Yvette, le 28 mars 2023, par

Tulio HONORIO DE FARIA

## Composition du jury

**Laurie LACARRIERE**
Professeur des universités,
INSA Toulouse

Présidente

**Christian HELLMICH**
Univ.-Prof. Dipl.Ing. Dr.techn.,
Vienna University of Technology

Rapporteur

**Roland J.-M. PELLENQ**
Directeur de Recherche,
CNRS, George Washington University

Rapporteur

**Gilles PIJAUDIER-CABOT**
Professeur des universités,
Université de Pau et des Pays de l'Adour

Rapporteur

**Farid BENBOUDJEMA**
Professeur des universités,
École Normale Supérieure Paris-Saclay

Examinateur

**Ian C. BOURG**
Professeur,
Princeton University

Examinateur

**Jorge S. DOLADO**
Scientist,
Spanish Council for Scientific Research,
Donostia International Physics Center

Examinateur

**Andrey KALINICHEV**
Directeur de Recherche,
Institut Mines-Télécom Atlantique

Examinateur
**Titre** : Modélisation Multi-Échelle Multi-Physique Multi-Technique pour la corrélation composition-(micro)structure-propriétés des matériaux cimentaires

**Mots clés** : Simulation moléculaire ; Micromécanique ; Matériaux de construction ; Matériaux poreux ; Échelle nanoscopique.

**Résumé** : Le paradigme central de la science des matériaux et de l'ingénierie des matériaux consiste à établir des corrélations entre les procèdes de fabrication, la (micro)structure, les propriétés et la performance des matériaux. Ces corrélations, en particulier les relations structure-propriété, sont fondamentales dans l'étude des matériaux cimentaires. L'étude et la compréhension de ces corrélations constituent une stratégie cruciale pour faire face à la montée en complexité des compositions minéralogiques et chimiques des systèmes cimentaires. Les approches de modélisation et de simulation apparaissent comme des outils importants pour établir ces corrélations, en particulier lorsqu'elles sont basées sur des principes scientifiques fondamentaux. Dans ce travail, je présente des techniques de modélisation et de simulation utilisées dans un cadre multi-échelles et multi-physiques dans le but d'établir des corrélations composition-propriété dans les systèmes cimentaires.

Ce manuscrit suit une organisation par échelle, avec le premier chapitre consacré aux informations de l'échelle moléculaire, suivi d'un chapitre sur la façon de mettre à l'échelle les informations provenant de l'échelle nanométrique jusqu'à l'échelle des applications de génie civil.

Tout d'abord, il est montré comment les simulations moléculaires peuvent être utilisées pour calculer les constantes élastiques, les propriétés thermiques, les propriétés de transport de masse, les propriétés électromagnétiques et les propriétés de couplage des phases pertinentes pour les systèmes cimentaires. Il est également montré comment ces simulations permettent de dévoiler les processus physiques associés à la sorption et aux couplages hygro-thermo-mécaniques dans une approche ascendente. Ensuite, une discussion sur les mésoéchelle(s) d'intérêt - dans le passage de l'échelle atomistique au continuum, mais aussi les mésoéchelles dans le domaine de la mécanique des milieux continus jusqu'à l'échelle macroscopique (celle des applications industrielles) des matériaux cimentaires – est présentée. On montre comment la théorie de la micromécanique fournit des estimations de propriétés effectives dans un contexte multi-physique mais aussi des informations sur les champs locaux. Les recherches en cours et les perspectives étendant la stratégie multi-échelle multi-physique et multi-technique à de nouveaux phénomènes et applications sont enfin présentées.

**Title** : Multi-Scale Multi-Physics Multi-Technique modeling for establishing composition-(micro)structure-property correlations in cement-based materials

**Keywords** : Molecular simulations ; Micromechanics ; Construction materials ; Porous materials ; Nanoscale.

**Abstract** : The central paradigm of materials science and materials engineering involves establishing processing-(micro)structure-property-performance correlations for materials. These correlations, in particular structure-property relationships, are fundamental in the study of cement-based materials, being a critical strategy for coping with the rise in the complexity of the mineralogical and chemical compositions of cement systems. Modeling and simulation approaches arise as important tools in establishing these correlations, especially when they are based on fundamental scientific principles. In this work, I present modeling and simulation techniques used in a multiscale and multi-physics framework with the goal of establishing composition-property correlations in cement systems.

This manuscript follows a scale-wise organization, with the first chapter devoted to information from the molecular scale, followed by a chapter on how to upscale nanoscale information up to the scale of Civil Engineering applications.

First, it is shown how molecular simulations can be used to calculate the elastic constants, thermal properties, mass transport properties, electromagnetic properties, and coupling properties of phases relevant to cement systems. It is also shown how molecular simulations help unveil the physical processes associated with sorption and hygro-thermo-mechanical couplings bottom-up. Then, a discussion on the mesoscale(s) of interest in the passage atomistic-to-continuum and in the realm of continuum mechanics up to the macroscopic scale (of industrial applications) of cement-based materials is presented. It is shown how Micromechanics theory provides estimates of effective properties in a multi-physics context but also information about local fields. Ongoing research and perspectives extending the multiscale multi-physics and multi-technique strategy to new phenomena and applications are finally presented.
# Contents

Contents i

List of Figures iii

List of Tables xiii

List of Symbols xv

1 Introduction 1

1 Multi-Scale and Multi-Technique 2

2 Multi-Physics 3

3 Outline 5

2 Information from the molecular scale 7

1 Introduction 8

2 Properties from molecular simulations 10

2.1 Structural features 10

2.2 Elastic constants 13

2.3 Thermal properties 16

2.4 Mass transport properties 21

2.5 Electromagnetic properties 34

2.6 Coupling properties 41

2.7 Acoustic response 45

3 Sorption processes and behavior of adsorbed water 45

3.1 Disjoining pressure and effective interactions 46

3.2 Hysteresis, cavitation and capillary condensation 50

3.3 Structure of water in microporous hydrates 53

3.4 Properties of confined and surface water 55

4 Conclusion 60

3 Upscaling properties: from nano to macro 61

1 Introduction 62

2 Mesoscale(s) 64

3 Coarse-grained simulations 64
<table><tbody><tr><td>3.1</td><td>ASR gels</td><td>64</td></tr><tr><td>3.2</td><td>C-S-H gels</td><td>68</td></tr><tr><td>4</td><td>Micromechanics</td><td>75</td></tr><tr><td>4.1</td><td>Representation of the microstructure</td><td>75</td></tr><tr><td>4.2</td><td>Conductivity-like properties</td><td>78</td></tr><tr><td>4.3</td><td>Elastic properties</td><td>86</td></tr><tr><td>4.4</td><td>Aging linear viscoelasticity</td><td>89</td></tr><tr><td>4.5</td><td>Thermal expansion</td><td>92</td></tr><tr><td>4.6</td><td>Heat capacity</td><td>93</td></tr><tr><td>4.7</td><td>Local information in elasticity and strength estimates</td><td>95</td></tr><tr><td>5</td><td>Conclusion</td><td>100</td></tr><tr><td>4</td><td>Conclusion and perspectives</td><td>101</td></tr><tr><td>1</td><td>Conclusions</td><td>101</td></tr><tr><td>2</td><td>Perspectives</td><td>101</td></tr><tr><td>2.1</td><td>Ongoing work</td><td>102</td></tr><tr><td>2.2</td><td>Other perspectives</td><td>107</td></tr><tr><td>A</td><td>Curriculum Vitae</td><td>111</td></tr><tr><td></td><td>Bibliography</td><td>129</td></tr></tbody></table>
List of Figures

1.1 Chain model of the central paradigm of materials science and engineering (adapted from [Olson, 1997]). 2

1.2 Time and space scales in modeling materials and structures: modeling and experimental techniques and the appropriate scales. 4

2.1 Snapshots of atomic structures of the systems studied: C-S-H [Honorio et al., 2021b, Honorio et al., 2022b] (Ca/Si=1.67 following [Kunhi Mohamed et al., 2018]), C-A-S-H [Hawchar and Honorio, 2022] (Ca/Si=1.7 following [Kunhi Mohamed et al., 2020]), tobermorite [Honorio, 2019] (Hamid's 11 Å tobermorite [Hamid, 1981]), crystalline ASR products [Honorio et al., 2020d, Honorio et al., 2020f, Honorio, 2021] (based on shlykovite [Zubkova et al., 2010]), montmorillonite [Honorio et al., 2017b, Honorio et al., 2018d, Honorio and Brochard, 2022], monosulfoaluminate [Honorio et al., 2020e], Friedel's salt [Honorio et al., 2022a], ettringite and metaettringite [Honorio et al., 2020e, Honorio et al., 2021a, Honorio, 2022b], bulk pore solutions [Honorio et al., 2019a, Honorio et al., 2020a], and clinker minerals (triclinic C₃S [Barbosa and Honorio, 2022] and C₃A [Honorio et al., 2020e]). 9

2.2 (A) XRD patterns of ettringite and metaettringite according to the RH at 300 K: comparison of simulation using AFFF [Honorio et al., 2021a] and data from the literature [Baquerizo et al., 2016, Skoblinskaya and Krasilnikov, 1975]. (B) Selected snapshots of simulated systems showing how the Ca₆Al₂(OH)₁₂ columns (indicated by the yellow shadows) in ettringite are not well defined in metaettringite. 12

2.3 Components of stiffness tensor (in Voigt notation) of Friedel's salt as a function of the temperature [Honorio et al., 2022a]. The main components are shown in (A)-(C), whereas (D) shows the components with a value close to 0 with respect to the standard deviation of the simulations (on the order of a few GPa). Savitzky-Golay filter was deployed to smooth simulation data in (A)-(D). (E) Elastic moduli and (F) Poisson ratio obtained from VRH approximation. 15
2.9 Shear viscosity of water confined in C-S-H as a function of the effective pore size *H* [Honorio, 2022a]: (A) Values computed for C-S-H from fitting NEMD velocity profiles, using Einstein-Stokes and (parallel to the pore wall) self-diffusion coefficient [Honorio et al., 2022b], and direct auto-correlation functions (Eq. 2.27); (B) Shear viscosity normalized with respect to the shear viscosity of bulk water ($\eta_{bulk}$, not to be confused with the bulk viscosity) comparison of results in (A) with experimental data from the literature on other nanoporous materials. The results of Low [Low, 1976] refer to interlayer distance in clays. The experimental studies [Haneveld et al., 2008, Persson et al., 2007, Thamdrup et al., 2007, Jarlgaard et al., 2006] refers to silica nanochannels. In both graphics, the value of the shear viscosity of bulk water is shown for comparison. 29

2.10 (A) Scheme showing a (parabolic) velocity profile from a (Poiseuille) flow confined between two surfaces with arbitrary textures described by the boundary conditions $b^+$ (upper) and $b^-$ (lower) (adapted from [Bazant and Vinogradova, 2008]). Each texture can be defined in their principal directions (here noted, 1 and 2); the textures are aligned when each pair {$b_1^+, b_1^-$} and {$b_2^+, b_2^-$} is colinear. (B) Negative slip, no-slip, and positive slip boundary conditions associated with hydrophilic, midwetting, and hydrophobic surfaces (adapted from [Zhan et al., 2020]). 31

2.11 Profiles of (A) density and (B) velocity for mesopores in C-S-H at 300 K for a pressure gradient applied along the x-direction [Honorio, 2022a]. The vertical dashed and solid lines indicate the position for the (negative) slip boundary; the distance between two vertical solid lines is the effective pore size according to the density profile of water oxygens. Density profiles are shifted one unit upwards for readability. Fittings (solid lines following the symbols in velocity profiles) using Eq. 2.31 are provided for each velocity profile. 32

2.12 (A) Intrinsic permeability calculated using Eq. 2.34 as a function of the interlayer space in slit channels with stick and slip boundary conditions [Honorio, 2022a]. The bilinear (defining micro- and meso-pore domains) and exponential fits of the slip length are considered. (B) Corrected intrinsic permeability $\kappa^c = f_\eta f_H \kappa$ calculated in slit channels with stick and slip boundary conditions. Only the exponential fit of slip lengths is considered. The inset zooms in the small micropores range. 33
2.17 Components of strain tensor obtained from simulations using the data from Fig.2.16 as a function of (A) the RH (with (B) zoom in the small strain domain), and (C) the water content (with (D) zoom in the small strain domain) [Honorio et al., 2021a]. Linear fitting of the small strain domain is provided for both cases. 44

2.18 Changes in the basal spacing of C-S-H quantified in GCMC-NVT simulations [Masara et al., 2023] compared against experimental data from the literature for C-S-H with different Ca/Si ratios during desorption [Gutteridge and Parrott, 1976, Smith and Bayliss, 1972, Bayliss, 1973]. 44

2.19 Disjoining pressure isotherms (top) and free energy profiles (bottom) as a function of the interlayer distance for various RH at 300 K during adsorption and desorption [Masara et al., 2023]. 48

2.20 Crystalline ASR products (Na- and K-shlykovite) controlled on the RH [Honorio et al., 2020d]: (a)-(b) water content, (c)-(d) disjoining pressure Π and (e)-(f) free energy Λ profiles as a function of the basal spacing for Na- and K-shlykovite (left and right, respectively) at 296 K and under an imposed water pressure $P_W(μ_W)$ of 1 atm. The insets in (c)-(f) zoom in on the oscillatory portions of the pressure and energy profiles. 49

2.21 (A) Water desorption isotherms for C-S-H basal distances between 13.7 and 25.3 Å and bulk water versus RH. The 13.7 Å basal distance pore is the equilibrium interlayer pore. The water content at a given RH ($N_w$) is normalized to the number of water molecules at 100 % RH (N). (B) Adsorption isotherms for C-S-H interlayer distances between 13.7 and 25.3 Å. (C) Sorption isotherms for the same range of C-S-H interlayer spacings and bulk water. The graphs are moved upward for clarity [Masara et al., 2023]. 52

2.22 Sorption isotherms of ettringite obtained from GCMC-NσT simulations [Honorio et al., 2021a] compared against experimental data [Skoblinskaya and Krasilnikov, 1975] ($T = 291$ K): (A) Desorption branch compared to the experimental sets "A"-"E" as labeled in [Skoblinskaya and Krasilnikov, 1975]. (B) and (C) Sorption cycles with close loop reversible hysteresis from 30 and 1 % RH (dark yellow labels), respectively. (D) Sorption cycle with open-loop hysteresis: simulation cannot capture the reversibility under adsorption. The simulated desorption and adsorption branches are fitted with Langmuir isotherms. The standard deviations of water contents obtained from simulations are on the order of the symbol size and are not shown for conciseness. 53
2.23 Radial distribution function of (A) Ow-Ow and (B) Hw-Hw pairs of water confined in C-S-H compared to bulk water [Honorio et al., 2021b]. The peaks indicating the nearest, second (first peak), second, third, etc., nearest correlations are labeled. The bifurcation of the second peak in Ow-Ow RDF for a c-length below 4 nm is also indicated. The positions of (C) the first peak in Ow-Ow RDF and (B) second peak in Hw-Hw RDF (the first peak is imposed by the equilibrium bond and angle in the respective harmonic potentials used in flexible SPC/E water). 54

2.24 Apparent density $\rho^*$; bulk modulus $K$, heat capacities $C_V$ and $C_P$; coefficient of thermal expansion $\alpha_{ii}$ and $\alpha_v = \text{Tr} \alpha_{ii}/3$; isosteric heat of adsorption $q_{st}$ as a function of the interlayer distance in C-S-H [Honorio et al., 2021b] compared to bulk SPC/E water and experimental data on bulk water. 57

2.25 (A) Sketch of the local-$\mu$VT conditions used to compute the properties of surface water using fluctuation formula in $\mu$VT ensemble [Honorio et al., 2022d]. (B) Initial layering (at time $t=0$) in C-S-H following a color code per layer and how the system evolves to a mixed state in a timescale of 100 ps. Layers were defined according to the density profiles of water oxygens Ow. 58

2.26 Properties per water layer according to two drained poromechanical conditions [Honorio et al., 2022d]: potential energy $U$, isochoric specific heat capacity $c_v$, molar incremental enthalpy $h_a$, bulk modulus $K$, and (volumetric) coefficient of thermal expansion $\alpha$. 59

3.1 Upscaling strategies in a multiscale multi-technique framework. 63

3.2 Snapshots of the mesostructures of Na-ASR and K-ASR gels as a function of the packing density $\eta_g$: networks (at top) and isosurfaces (at the bottom). 65

3.3 Eigenstresses evolution with densification of Na-and K-ASR gels for various packing densities [Ait-Hamadouche and Honorio, 2022b]. The PDFs are depicted in a log-scale for better identification of the behavior of the distribution tails. 66

3.4 Evolution of Young $E$ and shear $G$ moduli of Na-ASR (top) and K-ASR (bottom) as a function of the packing density at zero stress $\eta_0$ [Ait-Hamadouche and Honorio, 2022b]. The self-consistent estimate, using the Voigt-Reuss-Hill (VRH) bulk and shear moduli obtained from molecular simulations on Na and K-shlykovite as input, is provided for comparison. 67
3.5 Shear (left) and bulk (right) viscosities as a function of the packing density $\eta_g$ for Na- and K-ASR gels [Ait-Hamadouche and Honorio, 2022b]. The experimental data of [Yang et al., 2008] on static shear viscosity of sodium silicate gels (the gray curves) are shown for comparison. 67

3.6 (A) Separation of scales in nanolayers and stack of layers scales: conditions under which the assumptions of homogeneous plate theory holds [Honorio et al., 2018d]. (B) Representation of the "extreme" nanotextures of nanolayers materials (e.g., clays, C-S-H): (left) stacked or "brick-like" and (right) entangled or "spaghetti-like" types [Honorio et al., 2018d, Zabat et al., 1997]. (C) At left, sketch of a simple stack of clay layers undergoing a transition of hydration state from (n + 1)W to nW state. At right, a sketch of the transition between two hydration states and the superposition of solid layers: $l_{trans}$ characterizes the distance over which a layer bends in-between stack due to the interplay of flexibility and hydration forces [Honorio et al., 2018d]. 71

3.7 Layers aggregation (a) and isosurface (b) of C-S-H gel according to the flexibility of the layers for $\eta_g = 0.25$ [Ait-Hamadouche and Honorio, 2022a]. 74

3.8 Effective Young modulus of C-S-H a function of the packing density at zero stress $\eta_0$ according to the flexibility of the layers. The estimate of a self-consistent scheme for ellipsoidal inclusion with an aspect ratio of $a_r = 12.85/154.2 = 0.0833$ is provided for comparison. [Ait-Hamadouche and Honorio, 2022a] 74

3.9 Representation of the microstructure of concrete adopted in [Honorio et al., 2016b] to upscale aging linear viscoelastic properties: an approach combining analytical and numerical homogenization. 75

3.10 Cement paste microstructure considered for micromechanics upscaling schemes in [Honorio et al., 2022c]: Input volume fractions are obtained from Powers, Königsberger-Hellmich-Pichler (KHP) or Tennis and Jennings (TJ) models. SC or MT are deployed to upscale cement paste elastic properties. For 2-scale representations, C-S-H gel effective properties are upscaled using the SC scheme. All the other phases, including pores, are considered as spherical inclusions. SC scheme is deployed to upscale cement paste elastic properties. For the representations with ellipsoidal inclusions: with the Powers model, hydrates are considered as prolate particles $a_r = 10$; with the TJ model, C-S-H LD and HD are considered as prolate particles $a_r = 10$, and CH and AFm as oblate particles $a_r = 0.2$. 76
3.18 Micromechanics-based estimate performances for elastic prediction of the training set for some of the microstructure representations in Fig. 3.10 (only the representations with the best performance are shown here) [Honorio et al., 2022c]: experimental elastic constants from the training dataset (Young modulus *E*, Poisson ratio *ν*, bulk *K* and shear *G* moduli) plotted against the estimated values using various homogenization methods. Influence of the degree of hydration on the performances: the empty blue circles correspond to DOH <0.7, and the solid black circles correspond to DOH >0.7, i.e., predictions at late ages. 90

3.19 Basic creep at concrete scale [Honorio et al., 2016b]: comparison of homogenization estimates with experiments from [Craeye et al., 2009]. Mesostructures with spherical (full lines) and Voronoï (dashed lines) inclusions for scenarios with (right) and without (left) space-filling are shown. 91

3.20 Linear CTE of cement paste and mortar at early age [Honorio et al., 2018a]: estimates (full and dashed lines) compared to experimental data (dots) from [Wyrzykowski and Lura, 2013]. Full lines are "sealed" estimations, and dashed lines are "saturated" estimations. The corresponding estimations of the linear CTE of ITZ are shown in the upper inset. The lower inset shows the experimental heat flux as provided by [Wyrzykowski and Lura, 2013] in the context of COST TU 1404 Action. 93

3.21 Heat capacity for cement paste [Honorio et al., 2018a]: estimations (full and dashed lines) compared to experimental data (symbols) from [Mounanga, 2003] according to four different *w/c* ratios. Full lines are estimations without accounting for the coupling with thermal expansion; dashed lines are estimations accounting for the coupling with thermal expansion. 94

3.22 Effects of inclusion shape [Honorio et al., 2018b]: At the top, comparison of analytical estimates of stress distributions in matrix and inclusions with results from 3D finite element simulations on numerical samples with three different inclusion shapes. The Young moduli of the matrix and inclusion are $E_{mat} = 20$ GPa and $E_{inc} = 65$ GPa, respectively; the Poisson ratio of all phases is 0.2. A uniaxial load is applied in the x-direction. The analytical estimations are shown, assuming a normal distribution of stresses. At the bottom, numerical samples considered in 3D finite elements simulations: (a) spherical, (b) convex Voronoi polyhedral, and (c) convex Voronoi polyhedral flattened inclusions (with an aspect ratio of 3). 97
3.23 Comparison of analytical estimates [Honorio et al., 2018b] with experimental compressive strength of cement pastes obtained from the literature [Wyrzykowski et al., 2017, Honorio et al., 2016a, Pichler et al., 2013, Termkhajornkit et al., 2014, Lam et al., 2000]. At left, the results of Wyrzykowski et al. [Wyrzykowski et al., 2017] refer to a $w/c$ of 0.3; at the center, the results of Honorio et al. [Honorio et al., 2016a] refer to a $w/c$ of 0.5. In both cases, the evolution of Young moduli with age is presented in the respective articles. Analytical estimations are provided as a function of the critical volume fraction $f_{cr}$ and homogenization scheme (MT ou SC) used to get the effective properties of the hydration product layer scale. 99

4.1 Chemical potential of solid NaCl aqueous solution as a function of the molality [Hawchar, 2022]. The blue dashed lines correspond to Aragones et al. [Aragones et al., 2012] results. The solid magenta curves represent experimental results, and the red curves correspond to Hawchar Master thesis [Hawchar, 2022] simulations. Improved solubility estimations were obtained with new simulations using Yagasaki et al. [Yagasaki et al., 2020] potential. 105

4.2 Simulated aggregate volume fraction profile according to the granulometric curves of three concrete mix designs [Sayari et al., 2023]: F1 VeR-CoRs concrete with $D_{max}$ = 25 mm, F2 CETU concrete with $D_{max}$ = 25 mm and a different aggregates distribution, and F3 the B11 concrete with $D_{max}$ = 12.5 mm (the total volume fraction is the same and is equal to 40%). 106
<figure><img src="image_1.png" /></figure>
List of Symbols

$a, b, c$ three edge vectors

$\mathbf{A}_{pf}$ slip-driven plug flow factor

$\mathbf{A}_p$ (second-rank) localization tensor of phase $p$ for conductivity-like properties

$\mathbf{A}_p$ (fourth-rank) strain localization tensor of phase $p$ for elasticity

$\mathbf{A}_i^d$ strain localization tensor in the dilute case

$a_r$ aspect ratio

$b^+, b^-$ slip length of upper and lower surfaces, respectively

$\mathbb{B}_p$ (fourth-rank) stress localization tensor of phase $p$ for elasticity

$b_i$ apparent Biot coefficient

$\mathbf{B}_{sf}$ linear shear flow factor

$\bar{C}$ effective fourth-rank stiffness tensor

$c_j$ molar concentration of species $i$

$C_{ijkl} = \mathbb{C}$ (fourth-order) stiffness tensor

$C_{ij}$ Stiffness tensor in Voigt notation

$c_V$ specific heat capacity at constant volume

$c_P$ specific heat capacity at constant pressure

$C_P$ heat capacity at constant pressure

$C_V$ heat capacity at constant volume

$C_V^{class}$ heat capacity computed with the classical approach for the entire system

$C_{V-min}^{class}$ heat capacity of the dry mineral computed with the classical approach

$C_{V-min}^{semi-class}$ heat capacity of the dry mineral computed with the semi-classical approach

$d$ interlayer distance

$d_0$ characteristic size under which continuum mechanics is no more valid

$d_i$ characteristic length of the inhomogeneities

$D_i^f$ bending modulus

$d_{HB}$ distance of donor-acceptor oxygens

$\mathbf{D} = D_{ij}$ self-diffusion coefficient tensor

$D^{MSD}$ self-diffusion computed from the MSD

$D^{PBC}$ self-diffusion computed with periodic boundary conditions

$D^{VAC}$ self-diffusion computed from the velocity autocorrelation

$e$ elementary charge

$E$ Young modulus

$\mathbf{E}$ electric field
$E_a$ activation energy

$E_i$ internal energy

$E_{vib}$ vibrational energy

$E_{zz}$ out-of-pore-plane elastic modulus

$f$ frequency

$F$ Helmholtz free energy

$f_i(Q)$ weighting factor in $g_w(r)$ assumed to be equal to the effective number of electrons $Z$

$f_g$ gel porosity

$f_r$ volume fraction of phase $r$

$f_w$ mass fractions of water

$G$ shear modulus

$g(\omega)$ vibrational density of the states (VDOS)

$g_{i-j}(r)$ partial radial distribution function (RDF)

$g_w(r)$ weighted average of the partial pair distribution function

$h = [a, b, c]$ matrix grouping the 3 edge vectors

$h$ layer thickness

$H$ effective pore thickness

$h_a$ molar incremental enthalpy

$h_{bulk}$ molar enthalpy of the bulk fluid

$H_{cc}$ critical pore size for reversibility/continuity of capillary condensation transition

$H_{ch}$ critical pore size for the disappearance of hysteresis

$H_p$ enthalpy

$\hbar$ reduced Planck constant

$i$ imaginary unit

$I$ ionic strength

$I_i$ inertia moment

$\nabla i$ hydraulic gradient

$\mathbf{I}$ identity matrix

$\mathbf{I}$ identity fourth-order tensor

$\mathbf{J}$ vector flux

$J_i$ heat flux

$\mathbf{J}_I$ ionic current vector

$K$ bulk modulus

$\mathbf{K}$ hydraulic conductivity tensor

$k_{2b}, k_{3b}$ two- and three-body rigidities, respectively

$k_B$ Boltzmann constant

$k_g$ saddle-splay modulus

$K_g$ Gaussian curvature

$L$ length

$l_m$ mean free path of phonons

$l_p$ persistence length

$l_{REV}$ size of the Representative Elementary Volume (REV)

$l_{str}$ characteristic size of the structure in which the behavior of the material is to be considered as homogeneous
<table>
  <tr>
    <td>l<sub>trans</sub></td>
    <td>characteristic length associated with the interplay between the flexibility and the effective interactions</td>
  </tr>
  <tr>
    <td>m<sub>j</sub></td>
    <td>mass of particle j</td>
  </tr>
  <tr>
    <td>n(&omega;, T)</td>
    <td>Bose-Einstein distribution</td>
  </tr>
  <tr>
    <td>N</td>
    <td>number of particles</td>
  </tr>
  <tr>
    <td>N<sub>i-j</sub></td>
    <td>pair distribution function</td>
  </tr>
  <tr>
    <td>N<sub>t</sub></td>
    <td>total number of particles</td>
  </tr>
  <tr>
    <td><b>P</b></td>
    <td>polarization vector</td>
  </tr>
  <tr>
    <td>P<sub>ij</sub></td>
    <td>pressure tensor</td>
  </tr>
  <tr>
    <td>P<sub>co</sub></td>
    <td>cohesive pressure</td>
  </tr>
  <tr>
    <td>P<sub>w</sub></td>
    <td>water pressure</td>
  </tr>
  <tr>
    <td>q</td>
    <td>scattering vector</td>
  </tr>
  <tr>
    <td><b>q</b></td>
    <td>heat flux</td>
  </tr>
  <tr>
    <td>q<sub>i</sub></td>
    <td>charge pf particle i</td>
  </tr>
  <tr>
    <td>||</td>
    <td>stress fluctuation tensor</td>
  </tr>
  <tr>
    <td>q<sub>st</sub></td>
    <td>isosteric heat of adsorption</td>
  </tr>
  <tr>
    <td>r</td>
    <td>interparticle distance</td>
  </tr>
  <tr>
    <td>&#x2112;</td>
    <td>relaxation tensors</td>
  </tr>
  <tr>
    <td><b>r</b><sub>k</sub></td>
    <td>position vector of particle k of species i</td>
  </tr>
  <tr>
    <td>&lang;r<sub>i</sub>(t)<sup>2</sup>&rang;</td>
    <td>mean squared displacement (MSD)</td>
  </tr>
  <tr>
    <td>R</td>
    <td>gas constant</td>
  </tr>
  <tr>
    <td>S<sub>ijkl</sub> = &#x2112;</td>
    <td>(fourth-order) compliance tensor</td>
  </tr>
  <tr>
    <td>S<sub>2</sub></td>
    <td>pair (first-order) excess entropy</td>
  </tr>
  <tr>
    <td><b>S</b><sub>i</sub><sup>E</sup></td>
    <td>(second-rank) Eshelby-like tensor for conductivity-like properties</td>
  </tr>
  <tr>
    <td><b>S</b><sub>i</sub><sup>E</sup></td>
    <td>(fourth-rank) Eshelby tensor</td>
  </tr>
  <tr>
    <td>S<sub>ij</sub></td>
    <td>structure factor</td>
  </tr>
  <tr>
    <td>S<sub>l</sub></td>
    <td>longitudinal sound velocity</td>
  </tr>
  <tr>
    <td>S<sub>s</sub></td>
    <td>surface area</td>
  </tr>
  <tr>
    <td>S<sub>t</sub></td>
    <td>transverse sound velocity</td>
  </tr>
  <tr>
    <td>S<sub>v</sub></td>
    <td>volumetric sound velocity</td>
  </tr>
  <tr>
    <td>t</td>
    <td>time</td>
  </tr>
  <tr>
    <td>T</td>
    <td>temperature</td>
  </tr>
  <tr>
    <td>T<sub>cc</sub></td>
    <td>critical temperature for reversibility/continuity of capillary condensation transition</td>
  </tr>
  <tr>
    <td>T<sub>ch</sub></td>
    <td>critical temperature for the disappearance of hysteresis</td>
  </tr>
  <tr>
    <td><b>v</b><sub>j</sub></td>
    <td>velocity vector of particle j</td>
  </tr>
  <tr>
    <td>u<sub>i</sub></td>
    <td>displacement in direction i</td>
  </tr>
  <tr>
    <td>U<sub>2b</sub>, U<sub>3b</sub></td>
    <td>potential energy for two- and three-body interactions, respectively</td>
  </tr>
  <tr>
    <td>U<sub>LH</sub></td>
    <td>Landau-Helfrich bending free energy</td>
  </tr>
  <tr>
    <td>V</td>
    <td>volume</td>
  </tr>
  <tr>
    <td>V<sub>e</sub></td>
    <td>voltage</td>
  </tr>
  <tr>
    <td>V<sub>c</sub>(z)</td>
    <td>confinement potential</td>
  </tr>
  <tr>
    <td>V<sub>cav</sub></td>
    <td>volume of the cavity occupied by water</td>
  </tr>
  <tr>
    <td>x<sub>i</sub> = N<sub>i</sub>/N<sub>t</sub></td>
    <td>fraction of species i</td>
  </tr>
  <tr>
    <td>w/c</td>
    <td>water-to-cement ratio</td>
  </tr>
  <tr>
    <td>Z</td>
    <td>effective number of electrons</td>
  </tr>
  <tr>
    <td>z<sub>i</sub></td>
    <td>charge number of species i</td>
  </tr>
</table>
$\alpha_{ij}$ thermal expansion tensor

$\alpha_v$ volumetric thermal expansion

$\alpha_{CC}$ Cole-Cole exponent

$\gamma_s$ surface energy

$\gamma$ external tension

$\delta(x, y)$ displacement

$\epsilon_{ij}$ strain tensor

$\epsilon^{dry}$ drying shrinkage tensor

$\epsilon_{ij}$ dielectric permittivity tensor

$\epsilon_0$ vacuum permittivity

$\epsilon_{CC}$ Cole-Cole dielectric constant

$\eta$ viscosity

$\eta_0$ packing density realized at zero stress

$\eta_{bulk}$ bulk viscosity

$\eta_g$ packing density

$\phi_W, \phi_{IW}, \phi_I$ microscopic relaxation correlations function associated with water-water, ion-water and ion-ion contributions, respectively

$\kappa_0$ spontaneous curvature

$\bar{\kappa} = (\kappa_1 + \kappa_2)$ mean curvature in terms of principal curvatures $\kappa_1$ and $\kappa_2$

$\kappa^*$ dimensionless permeability

$\kappa^c$ intrinsic permeability

$\Lambda$ thermodynamic potential associated with drained experimental conditions (grand potential)

$\lambda_{ij}$ thermal conductivity tensor

$\lambda_{ij}^0$ static thermal conductivity

$\lambda_v$ volumetric thermal conductivity

$\mu$ chemical potential

$\nu_{ij}$ Poisson ratio

$\sigma_{ij}$ stress tensor

$\sigma_e$ electrical conductivity

$\sigma_{PS}$ electrical conductivity of the pore solution

$\Phi_{ij}$ heat flux auto-correlation function (HFACF)

$\Pi$ disjoining pressure

$\rho_0$ density number

$\rho$ density

$\tau_{CC}$ Cole-Cole relaxation time

$\tau_m$ relaxation time of phonons

$\theta_{HB}$ angle between O-H and O...O vectors in two water molecules

$\xi$ self-term for a cubic lattice at the room temperature

$\chi$ susceptibility

$\omega$ angular frequency
<table><tr><td>⟨★⟩<sub>ens</sub></td><td>ensemble average according to the ensemble "ens"</td></tr><tr><td>∇</td><td>gradient operator</td></tr><tr><td>⊗</td><td>dyadic product</td></tr><tr><td>[★]<sup>T</sup></td><td>matrix transpose operation</td></tr><tr><td>Tr(★)</td><td>trace operator</td></tr><tr><td>★<sub>∥</sub></td><td>parallel component</td></tr><tr><td>★<sub>⊥</sub></td><td>perpendicular component</td></tr><tr><td>=</td><td>complex conjugate</td></tr></table>
# Chapter 1

## Introduction

My research focuses on analytical and numerical methods to model the behavior of civil engineering materials and structures, mainly cement-based materials and concrete structures, in a multi-scale and multi-physics framework.

The central paradigm of materials science and materials engineering involves establishing materials Processing-(Micro)Structure-Property-Performance correlations [Olson, 1997]. These four primary elements can be understood as forming a chain (Fig. 1.1), allowing to

* *Understand and predict materials behavior* (material science relationships involving causes and effects): models based on the main mechanisms of the problem can be proposed for the understanding and prediction of behaviors, for example, by defining a certain property of the material as a function of its processing, composition and (micro-)structure.

* *Optimize and design materials with tailored properties or performance* (material engineering relationships involving objectives and means): from a given performance specification, inverse models are proposed for the optimization and design of materials. For example, a goal would be to meet a specific performance specification (e.g., mechanical strength), and one of the means would be to change the composition of the material.

These correlations, particularly structure-property relationships, are also fundamental in the study of cement-based materials [Bullard et al., 2019]. They are important for understanding the behavior of these materials and, from a more ambitious perspective, optimizing these materials for tailored performances. The "systematic development of structure-property relationships, through curation of fundamental material component data and validated modeling based on fundamental scientific principles" is recognized as a critical strategy for coping with the rise in the complexity of the mineralogical and chemical compositions of cement systems [Bullard et al., 2019]. Although modeling strategies are less frequently adopted than
<figure><img src="image_2.png" /><figcaption>**Figure 1.1:** Chain model of the central paradigm of materials science and engineering (adapted from [Olson, 1997]).</figcaption></figure>

experimental "trial-and-error exploration of the design space and mixture qualification process", they have the potential of "revolutionizing" concrete materials design, including design for service life [Bullard et al., 2019].

Modeling approaches serve to (i) predict the behavior of the materials and structures; (ii) understand its physical origin; (iii) extrapolate the predictions to conditions in which experiments can be prohibitive, and (iv) optimize the design of civil engineering materials and structures according to performance specifications or tailored properties. A **multi-scale-multi-technique** framework is needed to develop models based on fundamental scientific principles (e.g., [Vassaux et al., 2020]). Understanding the behavior of cement-based materials, or engineering them for tailored applications (e.g., for durability, for enhanced or new properties), requires a **multi-physics** approach. The relevance of Multi-Scale, Multi-Physics, and Multi-Technique modeling for establishing composition-(micro)structure-property correlations in cement-based materials is briefly discussed in the sequel.

# 1 Multi-Scale and Multi-Technique

In materials modeling, various scales can be critical for establishing Process-Structure-Property-Performance correlations. Figure 1.2 shows the relevant scales in materials science and engineering and how time and space scales are generally correlated in the context of multi-scale analysis. In a multi-scale approach, it is possible to associate relevant physical phenomena to the appropriate scale and then model the macroscopic behavior of the materials, properly accounting for their multi-scale nature. Following the example of what is observed in the experimental evaluation and characterization techniques, specific modeling and simulation techniques are more appropriate depending on the scale considered. A multi-technique framework is therefore required. The techniques used in my research are:
* Theoretical developments, including analytical homogenization and micromechanics [Honorio et al., 2016b, Honorio, 2017, Honorio et al., 2018a, Honorio et al., 2017a, Wyrzykowski et al., 2017, Honorio et al., 2018b, Honorio et al., 2020a, Honorio et al., 2020b, Königsberger et al., 2021, Honorio et al., 2022b];

* Numerical methods for continuum mechanics and structural mechanics [Honorio et al., 2014, Honorio et al., 2016a, Honorio et al., 2016b, Hervé-Secourgeon et al., 2021, Azenha et al., 2021];

* Classical atomistic (classical Monte Carlo and molecular dynamics) [Honorio et al., 2017b, Brochard et al., 2017, Honorio et al., 2018d, Honorio, 2019, Honorio et al., 2019b, Honório et al., 2019, Honorio et al., 2019a, Honorio et al., 2020e, Honorio et al., 2020a, Honorio et al., 2020f, Honorio et al., 2020d, Honorio et al., 2021b, Honorio, 2021, Honorio et al., 2021a, Honorio et al., 2022b, Honorio et al., 2022a, Barbosa and Honorio, 2022, Honorio, 2022b, Honorio and Brochard, 2022, Hawchar and Honorio, 2022] and coarse-grained simulations [Honorio, 2018, Ait-Hamadouche and Honorio, 2022a]; and more recently,

* Data science tools [Honorio et al., 2022c].

This manuscript will focus on Molecular Simulations and Micromechanics theory. Micromechanics is a powerful theory enabling for estimating the macroscopic properties of cement-based materials in a multi-physics framework using as input the system's composition, assumptions on the microstructure, and component property (i.e., the properties of each phase present in the material). Component property data is, however, lacking for some key phases in cement systems. Molecular simulations arise as a tool to provide such missing property data. These simulations provide a fundamental description of nanoscale processes at the origin of the various properties of materials. With this multi-scale and multi-technique approach, it is possible to tackle problems at the intersection of civil, mechanical, and materials engineering as well as physics and chemistry.

## 2 Multi-Physics

In my works, composition-property correlations have been investigated in a multi-physics framework, including

* **Mechanical properties** (elastic [Wyrzykowski et al., 2017, Honorio et al., 2020e] and viscoelastic properties [Honorio et al., 2016b, Honorio, 2017, Königsberger et al., 2021], mechanical strength [Honorio et al., 2018b, Honorio, 2021]). The mechanical behavior is crucial for most of the applications in which mechanical strengths and information on elastic response are key design parameters;
<figure><img src="image_3.png" /><figcaption>**Figure 1.2:** Time and space scales in modeling materials and structures: modeling and experimental techniques and the appropriate scales.</figcaption></figure>

* **Thermal properties** (heat capacity [Honorio et al., 2018a, Honorio et al., 2020d, Honorio et al., 2020f, Barbosa and Honorio, 2022, Hawchar and Honorio, 2022, Honorio, 2022b, Honorio et al., 2022a], thermal conductivity [Honorio et al., 2018a, Barbosa and Honorio, 2022, Honorio, 2022b]). The thermal behavior is critical for building envelopes aiming at ensuring thermal insulation (cement-based rendering mortars are often used for this purpose) and in massive structures at early-age (which are subjected to differential temperature rise due to the heat released by cement hydration that cannot be evacuated from the structure surfaces [Honorio, 2015]).

* **Mass transport properties** (diffusion [Honorio et al., 2019a, Honorio et al., 2022b] and permeability [Honorio, 2022a]). The mass transport behavior is crucial for understanding durability, which is directly affected by humidity and aggressive agent transport into the porous structure of cement-based materials.

* **Electromagnetic properties** (electrical conductivity [Honorio et al., 2020b, Honorio et al., 2020c] and dielectric permittivity [Honorio et al., 2020a]). The electromagnetic response is key to non-destructive evaluation and testing and structural health monitoring in civil engineering.

* **Coupling properties in thermo-poro-mechanics** [Honorio et al., 2018a,
Honorio, 2019]. The coupled phenomena (thermal expansion, drying shrink-
age, and poromechanics-related problems such as crystallization pressure de-
velopment) are also important to understanding the behavior and durability
of cement-based materials.

* **Sorption processes** [Honorio, 2019, Honorio et al., 2020d, Honorio et al., 2021a, Abdolhosseini Qomi et al., 2021, Honorio et al., 2021b], which are at the origin of various critical phenomena in cement-based materials, including drying shrinkage, reactivity, creep, etc.

Thus, a holistic view of cement-based materials properties accounting for these vari-
ous physics is critical to understand property development and aging, and durability.

The studies above focus mainly on cement-based materials, but I have worked
on other porous inorganic materials such as clays [Honorio et al., 2017b, Brochard
et al., 2017, Honorio et al., 2018c, Honorio et al., 2018d, Brochard and Honório,
2021, Bore et al., 2021, Honorio and Brochard, 2022] and hydroxyapatite [Honorio
et al., 2019b, Honório et al., 2019].

**3 Outline**

This manuscript follows a scale-wise organization, with the first chapter devoted to information from the molecular scale, followed by a chapter on how to upscale such nanoscale information up to the scale of Civil Engineering applications.

Chapter 1 shows how molecular simulations can be used to calculate the elastic constants, thermal properties, mass transport properties, electromagnetic properties, and coupling properties of phases relevant to cement systems. It is also shown how molecular simulations help unveil the physical processes associated with sorp- tion and hygro-thermo-mechanical couplings bottom-up.

Chapter 2 discusses the mesoscale(s) of interest in the passage atomistic-to-
continuum and the realm of continuum mechanics up to the macroscopic scale (of
industrial applications) of cement-based materials. It is presented how Microme-
chanics theory provides estimates of effective properties in a multi-physics context
but also information about local fields. It is also shown how coarse-grained simula-
tions give information on the mesocale between the atomic level and the micrometer
scale.

The manuscript finishes with the conclusions and perspectives.
# Chapter 2

## Information from the molecular scale

### Contents

<table><tr><td>1</td><td>Introduction</td><td>8</td></tr><tr><td>2</td><td>Properties from molecular simulations</td><td>10</td></tr><tr><td>2.1</td><td>Structural features</td><td>10</td></tr><tr><td>2.2</td><td>Elastic constants</td><td>13</td></tr><tr><td>2.3</td><td>Thermal properties</td><td>16</td></tr><tr><td>2.4</td><td>Mass transport properties</td><td>21</td></tr><tr><td>2.5</td><td>Electromagnetic properties</td><td>34</td></tr><tr><td>2.6</td><td>Coupling properties</td><td>41</td></tr><tr><td>2.7</td><td>Acoustic response</td><td>45</td></tr><tr><td>3</td><td>Sorption processes and behavior of adsorbed water</td><td>45</td></tr><tr><td>3.1</td><td>Disjoining pressure and effective interactions</td><td>46</td></tr><tr><td>3.2</td><td>Hysteresis, cavitation and capillary condensation</td><td>50</td></tr><tr><td>3.3</td><td>Structure of water in microporous hydrates</td><td>53</td></tr><tr><td>3.4</td><td>Properties of confined and surface water</td><td>55</td></tr><tr><td>4</td><td>Conclusion</td><td>60</td></tr></table>
# 1 Introduction

Constituent (or "phase" property) is one of the main inputs in multi-scale modeling. In cement-based materials (or more generally in complex multi-phase geomaterials), experimentally obtaining phases properties might be challenging because of

* The difficulty of isolating the phase in heterogeneous materials, in particular when these phases appear mixed with other with similar composition, as is often the case in cement systems. (Resorting to the synthesis of the phase is a solution to this problem but that brings new challenges, including knowing if the synthesized phase is equivalent to the phase found in cement-based materials.)

* Key phases exhibit compositional and structural variability, including notably C-(A)-S-H, AF-phases, and pore solutions [Vollpracht et al., 2016].

* The difficulty of probing the appropriate scale, which is the case of assessing the behavior of solid phases that occur as small crystallites or in small quantities (e.g., ASR products), quantifying the behavior of water in hydrates (including water confined in micro- and meso-pores or surface water), and quantifying ion sorption and dynamics.

The increasing complexity of chemical compositions of new cement systems make the presence of some *a priori* less common hydrate phases (e.g., M-S-H, other AF-phases beyond ettringite and monosulfoaluminate, C-S-H with alkali substitutions) more likely; and in some cases, these phases are even produced in significant quantities. This aspect encourages the effort of completing property data databases [Bullard et al., 2019] for the various potential phases that can be found in cement systems.

Furthermore, some basic properties are not yet experimentally measured even for some well-known phases in ordinary Portland cement (OPC) systems. This is the case of the thermal conductivity of ettringite (which is quite surprising knowing that the thermal problem is central in the durability problems involving ettringite [Taylor et al., 2001] and in the application of this mineral for thermal energy storage [Struble and Brown, 1984, Chen et al., 2019]). Also, surface-absorbed and confined water behavior, which is fundamental to understanding C-S-H behavior and durability processes in cement-based materials, remains elusive and difficult to assess experimentally [Abdolhosseini Qomi et al., 2021].

This chapter shows how molecular simulations have been used in my works to provide fundamental property data (elastic constants, thermal properties, mass transport properties, electromagnetic properties, and coupling properties) for various phases relevant to cement systems. Figure 2.1 shows the snapshots of the atomic structures of the systems studied. In some cases, molecular simulations have provided for the first time property data for the phase studied (e.g., some properties
<figure><img src="image_4.png" /><figcaption>**Figure 2.1:** Snapshots of atomic structures of the systems studied: C-S-H [Honorio et al., 2021b, Honorio et al., 2022b] (Ca/Si=1.67 following [Kunhi Mohamed et al., 2018]), C-A-S-H [Hawchar and Honorio, 2022] (Ca/Si=1.7 following [Kunhi Mohamed et al., 2020]), tobermorite [Honorio, 2019] (Hamid's 11 Å tobermorite [Hamid, 1981]), crystalline ASR products [Honorio et al., 2020d, Honorio et al., 2020f, Honorio, 2021] (based on shlykovite [Zubkova et al., 2010]), montmorillonite [Honorio et al., 2017b, Honorio et al., 2018d, Honorio and Brochard, 2022], monosulfoaluminate [Honorio et al., 2020e], Friedel's salt [Honorio et al., 2022a], ettringite and metaettringite [Honorio et al., 2020e, Honorio et al., 2021a, Honorio, 2022b], bulk pore solutions [Honorio et al., 2019a, Honorio et al., 2020a], and clinker minerals (triclinic C₃S [Barbosa and Honorio, 2022] and C₃A [Honorio et al., 2020e]).</figcaption></figure>
of crystalline ASR products, (meta)ettringite, Friedel's salt, C-A-S-H). The second part of this chapter will be devoted to the molecular modeling of sorption processes, which are critical for the drying and durability of cement-based materials. It will be shown how molecular simulations are a powerful tool to analyze complex aspects of sorption in cement systems, such as disjoining pressure and effective interactions, hysteresis (under sorption), the role of (heterogeneous) cavitation, and capillary condensation.

# 2 Properties from molecular simulations

Classical molecular simulations can be considered as a toolbox to compute properties once the atomic structure of the phase is known and relevant force fields are available.

Various empirical force fields have been developed for (or are appropriate for molecular modeling of) the various constituents of cement systems (see [Mishra et al., 2017] for a review). In my works, I have used ClayFF [Cygan et al., 2004, Cygan et al., 2021] and CSHFF [Shahsavari et al., 2010] to model calcium silicate phases. In the work on tricalcium silicate, we have tested various force fields [Barbosa and Honorio, 2022]. To better capture the structural features and elastic properties of AF-phases, we have proposed AFFF: a force field for AF-phases [Honorio et al., 2020e]. This force field is a reparametrization of ClayFF combined with parameters for sulfates of Cannon et al. [Cannon et al., 1994] adapted for AF-phases. An updated version of AFFF introduces harmonic covalent bonds and angles for sulfates [Honorio et al., 2021a]. AFFF parameters are available in LAMMPS format in https://github.com/tuliohf/Cement_Molecular/.

The primary strategies to get properties using molecular simulations are presented in the sequel. The most practical choices in simulation strategy are highlighted.

## 2.1 Structural features

The comparison of the structures obtained from molecular simulation with experimentally resolved structural descriptors provides a way of validating the choices of force fields. The structural features of the systems shown in Figure 2.1 were first validated against the available experimental data on various structural descriptors, including X-ray diffraction patterns, radial distribution functions, structural factors, lattice parameters, and atomic distances. Other interesting descriptors can be used to measure order/disorder (e.g., bond-orientational order [Steinhardt et al., 1983], excess entropy [Baranyai and Evans, 1989]) and assess water structure (H-bonds statistics). These structural descriptors are briefly presented in the following.

**Lattice parameters and interatomic distances** are basic parameters to verify the structural consistency of molecular models (atomic structure combined with force fields).
**Radial and pair distribution functions** provide information on the structural features of ordered and disordered materials that can be used to validate atomistic models and simulations. The partial radial distribution function $g_{i-j}(r)$ (RDF) of species i and j reads (e.g. [Thorneywork et al., 2014]):

$$g_{i-j}(r) = \frac{\rho_{i-j}(r)}{\rho_0} = \frac{1}{N_t \rho_0 x_i x_j} \sum_{k=1}^{N_i} \sum_{l \neq k}^{N_j} \delta\langle\mathbf{r} + \mathbf{r}_k - \mathbf{r}_l\rangle \quad (2.1)$$

where $N_t$ is the total number of particles, $\rho_0$ is the number density, $x_i = N_i/N_t$ is the fraction of species i (respectively, $x_j = N_j/N_t$), and $\mathbf{r}_k$ is the position vector of particle k of species i (respectively, $\mathbf{r}_l$ is the position vector of particle l of species j). The weighted average of the partial pair distribution function reads [Soyer-Uzun et al., 2012]:

$$g_w(r) = \sum_{ij} \frac{x_i x_j f_i(Q) f_j(Q) g_{i-j}(r)}{\langle f(Q) \rangle^2} \quad (2.2)$$

where $\langle f(Q) \rangle = \sum_i x_i f_i(Q)$, $f_i(Q)$ is assumed to be equal to the effective number of electrons Z for species i. The pair distribution function $N_{i-j}$ can be computed from $g_w$ by (e.g. [Fischer et al., 2006, Meral et al., 2011]):

$$N_{i-j} = 4\pi r^2 \rho_0 g_w(r) \quad (2.3)$$

In the case of C-S-H, radial and pair distribution functions obtained from simulations and experiments are compared in ref. [Honorio et al., 2021b]. The simulated structure exhibits similar features to the experimental peaks observed in previous studies in C-S-H with Ca/Si of 1.6 and 1.75 [Soyer-Uzun et al., 2012].

**Structure factor** quantifies material scattering and is obtained from the Fourier transform of radial distribution function:

$$S_{ij}(q) = 1 + 4\rho_0\phi \int_0^\infty (g_{i-j}(r) - 1) r^2 \frac{\sin(qr)}{qr} dr \quad (2.4)$$

where is $\rho_0$ is the average number density of particles, q is the scattering vector. This parameter has been used to evaluate medium-range order in C-S-H [Bauchy et al., 2014]. The structure factor was also computed for the C-S-H model studied in [Honorio et al., 2021b] and was compared with the experimental data provided by Soyer-Uzun et al. [Soyer-Uzun et al., 2012].

**X-ray diffraction (XRD) patterns** are related to the structure factor, the average number density of particles, and the orientation averaged form factor. XRD patterns have been computed for various systems in Figure 2.1 to validate the molecular modeling approach in yielding structurally sound systems. For example, Figure 2.2 shows the X-ray diffraction patterns obtained from simulations with AFFF at various RH at ambient temperatures compared with experimental patterns. Experimental and simulated patterns exhibit the main peaks at the same orientations.
The loss of well-defined peaks in metaettringite was fairly captured by the simulations. This example illustrates how simulation can provide structures consistent with the experimental data. This result builds confidence in using molecules simulation to explore the structural changes obtained under certain conditions that may be difficult to explore experimentally.

<figure><img src="image_5.png" /><figcaption>**Figure 2.2:** (A) XRD patterns of ettringite and metaettringite according to the RH at 300 K: comparison of simulation using AFFF [Honorio et al., 2021a] and data from the literature [Baquerizo et al., 2016, Skoblinskaya and Krasilnikov, 1975]. (B) Selected snapshots of simulated systems showing how the Ca₆Al₂(OH)₁₂ columns (indicated by the yellow shadows) in ettringite are not well defined in metaettringite.</figcaption></figure>

**Excess entropy** provides a quantification of the relative order in a molecular system, which is related to the structure, thermodynamics, and transport properties of liquid and disordered systems [Pacaud and Micoulaut, 2015]. The concept has been applied to measure order in C-S-H [Bauchy et al., 2014, Honorio et al., 2021b] and (meta)ettringite [Honorio et al., 2021a]. The pair (first-order) excess entropy S₂ can be computed from the pair radial function gᵢⱼ using [Baranyai and Evans, 1989, Pacaud and Micoulaut, 2015]:
$$ \frac{S_2}{k_B N} = -\frac{2\pi}{V} \sum_{i,j} \int_0^\infty [g_{i-j}(r) \ln g_{i-j}(r) - g_{i-j} + 1] r^2 dr \quad (2.5) $$

where N is the total number of particles, $k_B$ is the Boltzmann constant, V is the volume, and $x_i$ is the molecular fraction of particle of type i in the system. This formula allows calculating the pair-wise contribution $S_{2i-j}$ to the total entropy according to the choice of pairs i, j. In the case of C-S-H, the excess entropy is a non-monotonous function of the interlayer spacing d [Honorio et al., 2021b]: for $d \approx 5$ nm, the total, water, and electrolyte excess entropies are an increasing function, indicating that adding water to the system only contributes to the increase in its entropy. A (less intuitive) decrease in system excess entropy is observed in c-lengths ranging from 33 to 50 Å. Regarding AFt, the excess entropy is significantly higher in metaettringite when compared to ettringite [Honorio et al., 2021a], which is in agreement with the lack of well-defined peaks in metaettringite XRD patterns and the visual analysis of the snapshots in Figure 2.2.

H-bond statistics are other structural descriptors relevant to analyze hydrated systems. H-bond network statistics inform on the degree of constraints in water and have been linked to changes in the thermo-mechanical properties of layered silicates [Honorio et al., 2020f, Honorio et al., 2021b]. Various criteria exist to define an H-bond [Marcus, 2009], for example, a distance of donor-acceptor oxygens $2.5 < d_{HB} < 3.5$ Å and angle between O-H and O...O vectors $\theta_{HB} \le 100^{\circ}$. Changes in H-bonds statistics are analyzed in Section 3.3 for C-S-H with varying slit pore size thickness. As discussed in Section 3.2, H-bonds changes can be associated with the hysteresis under sorption in microporous materials [Chen et al., 2018], including ettringite [Honorio et al., 2021a].

## 2.2 Elastic constants

Two approaches are generally adopted to get the full stiffness tensor $C_{ijkl}$:

* Finite differences. Based on the definition [Allen and Tildesley, 1989]:

$$ C_{ijkl} = \left. \frac{\partial \sigma_{ij}}{\partial \epsilon_{kl}} \right|_{T, \epsilon_{kl}} \quad (2.6) $$

in which $\sigma$ and $\epsilon$ are the second-rank stress and strain tensors, and T is the temperature. The finite difference method can be deployed using energy minimization or under finite temperature (the former strategy is expected to lead to a more rigid system). With this approach, the simulation box is slightly deformed (within the elastic domain) in each one of the six axial and tangential directions, followed by an NVT run in which the ensemble-averaged pressure tensor is computed. Generally, the average of both negative and positive deformations is considered in the computations to account for asymmetric under loading.
* Fluctuation formula. $C_{ijkl}$ can be computed under finite temperature from a simulation in $N\sigma T$ ensemble using the fluctuation of the strain tensor with [Parrinello and Rahman, 1982]:

$$C_{ijkl} = \frac{kT}{\langle V \rangle_{N\sigma T}} [\langle \delta\epsilon_{ij} \otimes \delta\epsilon_{kl} \rangle_{N\sigma T}]^{-1}.$$

where $\langle \delta. \rangle_{ens}$ operator denotes ensemble average, and $\langle \delta A \star \delta B \rangle_{ens} = \langle A \star B \rangle_{ens} - \langle A \rangle_{ens} \star \langle B \rangle_{ens}$ denotes the fluctuation (co-variance) in ensemble $ens$ of random variables $A$ and $B$ according to a product $\star$ to be defined, in the case above " $\otimes$ " is the dyadic product. Even though a single simulation in $N\sigma T$ ensemble is needed with this approach, convergence must be carefully checked, and often long simulations are required.

The general definition of the Lagrangian deformation tensor $\epsilon_{ij}$ for a triclinic system is obtained from the matrices associated with the metric of the reference $h_0 = [a_0, b_0, c_0]_{x,y,z}$ and deformed $h = [a, b, c]_{x,y,z}$ supercell (i.e., the columns of $h$ matrices correspond to the components of the lattice vectors in an orthogonal frame $x, y, z$):

$$\epsilon_{ij} = \frac{1}{2} [(h_0^{-1})^T h^T h h_0^{-1} - \mathbf{I}], \qquad (2.7)$$

where $\mathbf{I}$ is the identity matrix, and the superscript $[\star]^{T}$ is a matrix transpose. This definition of deformation tensor is general because it can cope with the lowest symmetry system (triclinic) and finite deformations, but of course, the definition also holds for systems with higher symmetry (up to isotropy) and infinitesimal deformations.

Once the full stiffness tensor is obtained, the "isotropized" estimates of the elastic constants (bulk, shear, and Young moduli, as well as Poisson ratio) can be computed using the Voigt-Reuss-Hill (VRH) approximation [Hill, 1952]. The indentation modulus can then be computed from VRH moduli providing data for direct comparison with indentation experiments.

In my works, the elastic constants have been computed for the following phases in cement systems:

* **Clinker minerals.** Tricalcium silicate [Barbosa and Honorio, 2022] and tricalcium aluminate [Honorio et al., 2020e] have been studied. Full elastic tensor of the triclinic polymorph T1 of $C_3S$, which is the stable phase at ambient temperature of pure tricalcium silicate [Taylor, 1997], was provided for the first time. The values obtained for the T1 polymorph are similar to the ones found for M1 and M3 polymorphs (more common polymorphs found in industrial cement where substitutions occur in $C_3S$ structure) [Claverie et al., 2021].

* **C-(A)-S-H.** Full elastic tensors have been computed for C-S-H (S. Ait Hamadouche, Master thesis) and (no cross-linked) C-A-S-H [Hawchar and
Honorio, 2022]. The elastic constants obtained for C-A-S-H were close to the ones reported for C-S-H (based on defective tobermorite structures also with a Ca/Si of 1.6 to 1.7) [Shahsavari et al., 2010] but lower than the values reported for cross-linked C-A-S-H [Geng et al., 2018, Zheng et al., 2020]. The existence of (interlayer) cross-linking is C-A-S-H is under debate in the literature, with some authors arguing that no experimental evidence clearly indicates the presence of cross-linked C-(A)-S-H [Kunhi Mohamed et al., 2020].

<figure><img src="image_6.png" /><figcaption>**Figure 2.3:** Components of stiffness tensor (in Voigt notation) of Friedel's salt as a function of the temperature [Honorio et al., 2022a]. The main components are shown in (A)-(C), whereas (D) shows the components with a value close to 0 with respect to the standard deviation of the simulations (on the order of a few GPa). Savitzky-Golay filter was deployed to smooth simulation data in (A)-(D). (E) Elastic moduli and (F) Poisson ratio obtained from VRH approximation.</figcaption></figure>
* **AF-phases.** Full elastic tensors have been computed for ettringite, and metaettringite [Honorio et al., 2020e, Honorio, 2022b], monosulfoaluminate [Honorio et al., 2020e], and Friedel's salt [Honorio et al., 2022a] using AFFF in all cases. Simulations were performed for various water contents for ettringite, and monosulfoaluminate [Honorio et al., 2020e]. The changes in the elastic constants in the transition of ettringite into collapsing ettringite followed by transformation into metaettringite were also quantified [Honorio, 2022b]. In the cases of metaettringite and Friedel's salt, the molecular simulations were the first estimates of the elastic constants in the literature. (Note that structural features and interfacial processes in Friedel's salt has been extensively studied with molecular simulations [Kalinichev and Kirkpatrick, 2002, Kalinichev et al., 2000, Kirkpatrick, 2005, Kirkpatrick et al., 1999]). Friedel's salt exhibits a low temperature (LT), and a high temperature (HT) polymorph with a transition occurring at approximately 35°C [Terzis et al., 1987, Renaudin et al., 1999]. The simulation showed that some elastic constants, when plotted as a function of the temperature, display a significant change in the vicinity of 35°C in Friedel's salt, as shown in Figure 2.3.

* **ASR products.** Full elastic tensors were computed for the first time for crystalline ASR products based on shlykovite structure [Honorio et al., 2020f]. Isomorphic substitutions of the charge balance ions Li, Na, K, and Cs were considered in the simulations. The elastic constants obtained were affected by the substitution of charge balance ions. The VRH moduli of Na and K products (as well as the upscaled values considering the gel scale as a higher level - see next chapter for details on upscaling strategies) compare well with the nanoindentation data.

Most of the studies above deploy the finite difference approach (fluctuation approach has been used in ref. [Honório et al., 2019]).

## 2.3 Thermal properties

### 2.3.1 Heat capacity

Molecular simulations of the heat capacity of fluids can be performed in a full classical framework, whereas the heat capacity of solids might require a multi-scale approach in which quantum corrections are introduced.

In a **full classical approach**, the heat capacity can be computed using the finite difference or fluctuations approaches:

* *Finite difference* where the definitions of isochoric and isobaric heat capacity are used:

$$C_V = \left. \frac{\partial E_i}{\partial T} \right|_V; \quad C_P = \left. \frac{\partial H_p}{\partial T} \right|_P \qquad (2.8)$$
where V and P are the volume and pressure, $E_i$ is the internal energy, and $H_p = E_i + PV$ is the enthalpy.

* Fluctuation formula are given by [Allen and Tildesley, 1989]:

$$C_V = \frac{\langle \delta E_i^2 \rangle_{NVT}}{k_B T^2} \qquad (2.9)$$

$$C_P = \frac{\langle \delta H_p^2 \rangle_{NPT}}{k_B T^2} \qquad (2.10)$$

In a semi-classical approach, quantum effects are accounted for, but the vibrational density of states are obtained using classical methods [Dickey and Paskin, 1969, Berens et al., 1998]. In this approach, the heat capacity at constant volume is [Winkler and Dove, 1992]:

$$C_V = \frac{dE_{vib}}{dT} = 3Nk_B \int_0^\infty \left( \frac{\hbar\omega n(\omega, T)}{kT} \right)^2 \exp\left[\frac{\hbar\omega}{kT}\right] g(\omega) d\omega \qquad (2.11)$$

where $3Nk_B$ corresponds to the Dulong-Petit law, $\omega$ is the angular frequency, $\hbar$ is the reduced Planck constant, $n(\omega, T) = \frac{1}{\exp[\frac{\hbar\omega}{kT}] - 1}$ is Bose-Einstein distribution, and $g(\omega)$ is the vibrational density of the states (VDOS) [Dickey and Paskin, 1969, Berens et al., 1998]:

$$g(\omega) = \frac{1}{\pi} \int_{-\infty}^{\infty} \sum_{j} \frac{\langle m_j \mathbf{v}_j(t) \cdot \mathbf{v}_j(0) \rangle}{\langle m_j \mathbf{v}_j^2(0) \rangle} \exp[i\omega t] dt \qquad (2.12)$$

The VDOS is obtained from the mass-weighted velocity auto-correlation function $\langle m_j \mathbf{v}_j(t) \cdot \mathbf{v}_j(0) \rangle$, with $m_j$ being the mass of particle $j$; and $\mathbf{v}_j(t)$ being the velocity vector of particle $j$ at time $t$ (see the supporting material in [Honorio and Brochard, 2022] for calculation details). Various authors show that the incorporation of quantum corrections is often required to reproduce the experimental heat capacity of crystalline materials correctly [Winkler and Dove, 1992, Oganov et al., 2000, Qomi et al., 2015, Claverie et al., 2021]. Due to the high phonon frequency in crystals (i.e., when $\hbar\omega \sim k_B T$ in ambient conditions), quantum effects become non-negligible, leading to a heat capacity that should be lower than the classical estimate may yield. Semi-classical simulations are to be performed in the NVE ensemble to avoid the fictitious forces from the thermostats. Timesteps of at least 0.1 fs are needed to capture the oscillatory allure of mass-weighted velocity autocorrelation functions properly. To illustrate this approach, Figure 2.4 (A) and (B) display the normalized total velocity auto-correlation $\langle \mathbf{v}_j(t) \cdot \mathbf{v}_j(0) \rangle / \langle \mathbf{v}_j(0)^2 \rangle$, and normalized mass-weighted velocity auto-correlation $\langle m_j \mathbf{v}_j(t) \cdot \mathbf{v}_j(0) \rangle / \langle m_j \mathbf{v}_j(0)^2 \rangle$ functions of ettringite. These functions were used to calculate the VDOS $g(\omega)$ (Figure 2.4 (C)).

A hybrid approach combining the semi-classical method for the dry mineral and a classical approach for the water has been shown to best capture the heat capacity of hydrated minerals (with molecular water), including swelling clays [Honorio and Brochard, 2022], C-A-S-H [Hawchar and Honorio, 2022] and AFms [Honorio
<figure><img src="image_7.png" /><figcaption>Figure 2.4: (A) Normalized total velocity auto-correlation $\langle \mathbf{v}_j(t) \cdot \mathbf{v}_j(0) \rangle / \langle \mathbf{v}_j(0)^2 \rangle$, (B) normalized mass-weighted velocity auto-correlation $\langle m_j \mathbf{v}_j(t) \cdot \mathbf{v}_j(0) \rangle / \langle m_j \mathbf{v}_j(0)^2 \rangle$, and (C) VDOS $g(\omega)$ per RH at ambient temperature (298 K) [Honorio, 2022b].</figcaption></figure>

et al., 2022a]. A way of computing the heat capacity via this hybrid approach is first to compute the classical approach for the entire system ($C_V^{class}$) and then subtract the contribution of the dry mineral (the classical approach applied to the dry mineral yields Dulong-Petit law ($C_{V-min}^{class} = 3N_s k_B$ [Honorio and Brochard, 2022]). The contribution of water in the classical approach is therefore $C_{V-w}^{class} = C_V^{class} - 3N_s k$. The total heat capacity is then obtained as the sum of the inter-layer, and dry mineral ($C_{V-min}^{semi-class}$) contributions [Honorio and Brochard, 2022]:

$$C_V = C_{V-w}^{class} + C_{V-min}^{semi-class} = C_V^{class} - 3N_s k_B + C_{V-min}^{semi-class} \quad (2.13)$$

where $C_{V-min}^{semi-class}$ is computed from the semi-classical approach detailed above applied only to the dry mineral.

For comparison with experimental data, it is helpful reporting the heat capacities in terms of specific heat capacities $c_V$, $c_P$ (in a unit of energy per mass.Kelvin) The specific heat capacities can be computed from the total heat capacity $C_V$, $C_P$ using $c_V = C_V / (m_{total})$ and $c_P = C_P / (m_{total})$, with $m_{total}$ being the mass.

The difference between $c_P$ and $c_V$ can be computed using Mayer's relation (e.g. [Wallace, 1998]):

$$c_P - c_V = \frac{T}{\rho} C_{ijkl} \alpha_{ij} \alpha_{kl}, \quad (2.14)$$

where $\rho$ is the density, and $\alpha_{ij}$ is the tensor of thermal expansion.

In my works, the heat capacity has been computed for the following phases in cement systems:

* **Clinker minerals.** The heat capacities of triclinic tricalcium silicate [Barbosa and Honorio, 2022] were computed, for the first time, using full classical and
semi-classical approaches. The semi-classical yielded values similar to the one obtained for the monoclinic polymorphs, with which the T1 structure shares a superstructure. These results corroborate previous studies [Claverie et al., 2021, Winkler and Dove, 1992, Oganov et al., 2000, Qomi et al., 2015] showing that the estimates accounting at least partly for quantum effects are required to capture the heat capacity of crystals properly.

*   **C-(A)-S-H.** The heat capacities were computed for C-A-S-H [Hawchar and Honorio, 2022] using the full classical, semi-classical, and hybrid approaches. The latter two approaches yielded results closer to the available data from the literature [Roosz et al., 2018].

*   **AF-phases.** The heat capacities were computed for ettringite and metaettringite [Honorio, 2022b] and Friedel's salt [Honorio et al., 2022a] using AFFF. The hybrid approach returned results more consistent with the available data from the literature. Figure 2.5 compares full classical, semi-classical, and hybrid approaches estimates for ettringite and metaettringite [Honorio, 2022b].

<figure><img src="image_8.png" /><figcaption>**Figure 2.5:** A) Comparison of the specific heat capacity $c_V$ at constant volume of ettringite and metaettringite obtained using full classical methods (when rigid or flexible SPC/E water models are adopted), semi-classical and hybrid approaches [Honorio, 2022b]. Dulong-Petit estimate ($3Nk_B$) is shown for comparison. (B) Heat capacity $c_P$ at constant pressure, according to the hybrid approach, of ettringite (and metaettringite) as a function of the RH compared with experimental data from refs [Ederova and Satava, 1979, Baquerizo et al., 2016]. (C) $c_P - c_V$ of ettringite (and metaettringite) function of the RH computed from the full rigidity tensor and tensor of thermal expansion (only the hybrid approach is shown). All results refer to ambient temperature (298 K).</figcaption></figure>
* **ASR products.** The heat capacities of crystalline ASR products based on the shlykovite structure [Honorio et al., 2020f] were computed using the full classical approach for the first time. The values obtained were affected by the charge balance ions substitutions (Li, Na, K, and Cs). Estimates with the semi-classical and hybrid approaches are yet to be worked out.

### 2.3.2 Thermal conductivity

The thermal conductivity tensor can be computed from molecular simulations using equilibrium and non-equilibrium methods. Non-equilibrium methods include the direct thermostating method [Ikeshoji and Hafskjold, 1994], and the reverse perturbation method [Müller-Plathe, 1997]. The main equilibrium method is the Green-Kubo formalism, detailed in the following.

Using Green-Kubo formalism, the thermal conductivity is the time integral of the ensemble average of the heat flux $J_i$ ($i = x, y, z$) auto-correlation (e.g., [Boon and Yip, 1991]):

$$ \lambda_{ij} = \frac{V}{k_B T^2} \int_0^\infty \langle J_i(0) \otimes J_j(t) \rangle dt \qquad (2.15) $$

The terms $\Phi_{ij} = \langle J_i(0) \otimes J_j(t) \rangle$ are the components of the heat flux auto-correlation function (HFACF). The heat flux vector is computed from the viral stress, kinetics, and potential energy per atom:

$$ \mathbf{J}_p = \frac{1}{V} \left[ \sum_p E_p \mathbf{v}_p - \sum_p \sigma_p \mathbf{v}_p \right] \qquad (2.16) $$

where $E_p$ is the sum of the kinetics and potential energy, $\mathbf{v}$ is the velocity vector, and $\sigma_p$ is the stress tensor for atom $p$.

The conductivity $\lambda_{ij}$ computed from Eq. 2.15 might be prone to fluctuations and convergence problems as $t \to \infty$ in the integration. A solution to overcome this problem is to compute the thermal conductivity tensor using the frequency-dependent response of the complex conductivity [Termentzidis et al., 2011]:

$$ \lambda_{ij}(\omega) = \frac{V}{k_B T^2} \int_0^\infty \langle J_i(0) \otimes J_j(t) \rangle \exp[i\omega t] dt \qquad (2.17) $$

The static thermal conductivity $\lambda_{ij}^0$ can be obtained from fitting, in the low-frequency range, the modulus $|\lambda_{ij}(\omega)|$ of the complex conductivity with the expression $|\lambda_{ij}(\omega)| = \lambda_{ij}^0 + A\omega^\xi$, where $\lambda_{ij}^0$, $A$ and $\xi$ are free parameters in the fitting procedure [Qomi et al., 2015].

In my works, the thermal conductivity tensor has been computed for the following phases in cement systems:

* **Clinker minerals.** The full thermal conductivity tensor for T1 tricalcium silicate was computed for the first time [Barbosa and Honorio, 2022]. The calcula-
tion was made using the Green-Kubo formalism accounting for the frequency-dependent response of the complex conductivity. Compared to previous results on the monoclinic polymorph [Qomi et al., 2015], our results suggest improvement in molecular modeling of the thermal conductivity of the C₃S since we obtain a volumetric thermal conductivity in better agreement with the experimental values reported for the cement powder of 1.55 W/(m.K) [Bentz, 2007].

* **AF-phases.** The full thermal conductivity tensor of ettringite and metaettringite [Honorio, 2022b] were computed for the first time using AFFF. Estimates obtained using standard Green-Kubo formalism yielded similar results to the ones obtained from the complex frequency-dependent thermal conductivity. On average, the diagonal components of the tensor of thermal conductivity of ettringite (in a Cartesian frame) follow $\lambda_{xx} > \lambda_{yy} > \lambda_{zz}$ ($3.09 \pm 0.90 > 2.83 \pm 0.85 > 2.41 \pm 0.88$ W/(m.K)). The out-of-diagonal components $\lambda_{ij}$ ($i \neq j$) are close to zero (considering the standard deviation of the simulations), suggesting that the Cartesian frame approximately corresponds to the principal frame. In the ettringite domain, volumetric thermal conductivity $\lambda_v = \text{Tr}(\lambda_{ij})/3 = \sum \lambda_{ii}/3$ is $2.76 \pm 0.51$ W/(m.K). This value can be compared with the conductivity of ettringite recently reported in the literature of $2.19 \pm 0.4$ W/(m.K) [Du and Ge, 2021]. The values obtained at low RH are: for collapsing ettringite $\lambda_v = 4.45 \pm 0.73$ W/(m.K) (a quite high value considering the typical conductivity of hydrated minerals; further studies are thus needed to validate this result), and for metaettringite $\lambda_v = 2.06 \pm 0.65$ W/(m.K).

The thermal conductivity of C-S-H, jennite, and tobermorites were computed using molecular simulations [Qomi et al., 2015, Sarkar and Mitra, 2021, Hong et al., 2020]. The thermal conductivity of ASR products and C-A-S-H is yet to be worked out.

## 2.4 Mass transport properties

### 2.4.1 Diffusion

**Bulk fluids.** The dynamics of bulk fluids can be quantified in terms of the Mean Squared Displacement (MSD):

$$ \langle r_i(t)^2 \rangle \approx D_i t^\alpha \quad \text{with } i = x, y \text{ or } z \qquad (2.18) $$

where $D_{ii}$ is the self-diffusion coefficient, $r_i(t)$ are the components of the displacement vector $\mathbf{r}(t)$ at time $t$, and $\alpha$ is the anomalous diffusion exponent: $\alpha = 1$ for Fickian diffusion; anomalous regimes are classified in subdiffusion for $0 < \alpha < 1$ and superdiffusion for $\alpha > 1$. Fickian diffusion corresponds to a random walk, while anomalous diffusion results from a break of ergodicity in the system's behavior, which can have different physical origins such as fluid confinement, different
particle sizes, and aging phenomena [Metzler et al., 2014]. The MSD can be directly computed from an MD simulation for a group of particles, and the diffusion regimes can be discussed as a function of the particle type.

In isotropic systems showing Fickean dynamics, the self-diffusion coefficient can be computed using Einstein's equation for 3D systems [Allen and Tildesley, 1989]:

$$D^{\text{MSD}} = \frac{1}{6} \lim_{t \to \infty} \frac{\langle |\mathbf{r}(t)|^2 \rangle}{t} \qquad (2.19)$$

where the operator $| \cdot |$ denotes the module of a vector.

Self-diffusion coefficient can also be computed using the velocity auto-correlation (VAC) function via the Green-Kubo relation [Allen and Tildesley, 1989]:

$$D^{\text{VAC}} = \frac{1}{3} \int_{0}^{\infty} \langle \mathbf{v}(0) \cdot \mathbf{v}(t) \rangle dt \qquad (2.20)$$

where the velocity $\mathbf{v}(t)$ is directly sampled for a group of atoms with MD simulations.

MD simulations with periodic boundary conditions (PBC) introduce artifacts that render the computation of self-diffusion coefficients dependent on the simulation box size. In a cubic simulation box of length $L$, the corrected diffusion coefficient $D(L \to \infty)$ can be obtained from the diffusion coefficient $D^{\text{PBC}}$ computed with periodic boundary conditions using [Yeh and Hummer, 2004]:

$$D(L \to \infty) = D^{\text{PBC}} + \frac{k_B T \xi}{6 \pi \eta L} \qquad (2.21)$$

where $\xi \approx 2.837298$ is the self-term for a cubic lattice at room temperature, and $\eta$ is the viscosity of the fluid, which in the case of SPC/E water is reported to be insensitive to simulation cell size [Yeh and Hummer, 2004].

MSD and VAC methods were used to compute the self-diffusion coefficients of water and ions in bulk aqueous solutions mimicking the composition of the pore solutions in cement-based materials [Honorio et al., 2019a]. It is relevant to investigate bulk solutions because the most significant part of the dynamics of water and ions in pore solutions within cement-based materials can be described by the dynamics of the fluid portions less affected by pore walls. The ionic transport occurring in the pore zones affected by the interface (Stern layer) is reported to be negligible when compared to transport in the bulk pore, especially for aqueous solutions with a concentration similar to pore solutions in cement-based materials [Samson et al., 2005, Revil, 1999].

The compositions of the systems studied are experimental data [Lothenbach et al., 2008] for ordinary Portland cement (PC) systems in a formulation with a water-to-cement $w/c$ ratio of 0.4. Seven scenarios PC1-PC7 corresponding to various ages are chosen to be simulated. The species with a concentration below 5 mM are not taken into account. The ionic concentrations were adjusted to comply with the electroneutrality of the simulated systems. The resulting aqueous solutions compositions are within the variability experimentally observed in the pore solution
<figure><img src="image_9.png" /><figcaption>**Figure 2.6:** Evolution of the concentrations of the main ionic species (concentrations exceeding 5 mM) in ordinary Portland cement (PC) pore solutions [Honorio et al., 2019a, Honorio et al., 2020a]. The scenarios PC1-PC7 are retained (snapshots of atomic configurations at the bottom) to study the age dependency of the pore solution composition. PC1-PC2 corresponds to *very early-age*, PC3-PC4 to *early-age* and PC5-PC7 to *late ages* typical compositions for OPC system with a *w/c* of 0.4 [Vollpracht et al., 2016].</figcaption></figure>

of cement systems [Vollpracht et al., 2016]. Figure 2.6 shows the evolution of the concentration of the ionic species plotted as a function of age. Three main stages are identified regarding the ionic composition of the pore solution for a given cement system [Honorio et al., 2019a]:

* *very early age* up to setting (first hours): PC1 and PC2;

* *early-age* related to property development (within the first days): PC3 and PC4;

* *late stage* associated with the service life of cement-based materials: PC5 to PC7.

Such classification has been useful in linking the dynamics of ions to the electrical properties of cement-based materials [Honorio et al., 2020b]. A snapshot of the atomic configurations in a system equilibrated in an NVT simulation at 300 K is shown for each age.

Table 2.1 gathers the self-diffusion coefficients computed using Eq. 2.19. These results agree with data from the literature regarding the self-diffusion of water and ions according to simulation, and experiments [Tinnacher et al., 2016, Honorio et al., 2019b, Hammond, 1997, Mark and Nilsson, 2001, Zangi and Engberts, 2005]. Calcium ions present slower dynamics, which can be associated with their size and higher valence than other ions. Hydroxide and sulfate ions also exhibit dynamics slower than water. The ionic compositions in PC2 and PC3 slightly decelerate water dynamics, while the ionic compositions in PC7 slightly accelerate water dynamics;
<table><caption>**Table 2.1:** Self-diffusion coefficients $D_i$ (in $10^{-9}$ m²/s) at 300 K and 1 atm of species *i* computed using Eq. 2.19 [Honorio et al., 2019a]. For comparison, we provide data from the literature: <sup>a</sup> [Tinnacher et al., 2016] (Sim.), <sup>b</sup> [Honorio et al., 2019b](Sim.), <sup>c</sup>[Hammond, 1997](Exp. infinitely diluted), <sup>d</sup>[Mark and Nilsson, 2001] (Sim.), <sup>e</sup> [Zangi and Engberts, 2005] (Sim.)</caption><thead><tr><th></th><th>Water</th><th>Na<sup>+</sup></th><th>K<sup>+</sup></th><th>Ca<sup>2+</sup></th><th>SO<sub>4</sub><sup>2-</sup></th><th>OH<sup>-1</sup></th></tr></thead><tbody><tr><th>PC1</th><td>2.83</td><td>1.32</td><td>1.42</td><td>0.34</td><td>0.93</td><td>1.77</td></tr><tr><th>PC2</th><td>2.22</td><td>1.67</td><td>1.50</td><td>0.88</td><td>1.55</td><td>0.88</td></tr><tr><th>PC3</th><td>2.52</td><td>2.24</td><td>2.56</td><td>-</td><td>-</td><td>0.76</td></tr><tr><th>PC4</th><td>2.67</td><td>2.77</td><td>2.37</td><td>-</td><td>-</td><td>1.57</td></tr><tr><th>PC5</th><td>2.92</td><td>3.76</td><td>4.09</td><td>-</td><td>2.27</td><td>1.69</td></tr><tr><th>PC6</th><td>2.93</td><td>2.28</td><td>3.19</td><td>-</td><td>1.73</td><td>1.56</td></tr><tr><th>PC7</th><td>3.14</td><td>2.01</td><td>2.58</td><td>-</td><td>1.71</td><td>1.37</td></tr><tr><th>Data<br>from the<br>literature</th><td>2.68<sup>a</sup><br>2.70, 3.02<sup>d</sup><br>2.3<sup>c</sup></td><td>1.08<sup>a</sup><br>1.334<sup>c</sup></td><td>1.957<sup>c</sup></td><td>0.85<sup>a</sup><br>0.792<sup>c</sup><br>0.797<sup>b</sup></td><td>1.065<sup>c</sup><br>0.752<sup>b</sup></td><td>1.38<sup>e</sup><br>5.273<sup>c</sup></td></tr></tbody></table>

water dynamics is similar to the dynamics of pure water for the other pore solutions. The dynamics of ions are generally much more variable depending on the pore solutions. Sodium ions, for example, exhibit dynamics similar to that of a dilute solution at the very early-age (PC1 and PC2) and faster dynamics at PC3 and later ages. Sulfate and potassium also show dynamics acceleration with the evolution of the composition of the solution with age. The differences in solution dynamics and structure were also evaluated by analyzing the normalized velocity auto-correlation functions and non-Gaussian parameters, water orientation time correlation functions, and H-bonds counting in ref. [Honorio et al., 2019a].

**Confined fluids.** Self-diffusion coefficients under confinement are subject to anisotropic response. In the case of C-S-H (or other layered phases), the tensor of diffusion has the form:

$$ \mathbf{D}(z) = \begin{pmatrix} D_{\parallel}(z) & 0 & 0 \\ 0 & D_{\parallel}(z) & 0 \\ 0 & 0 & D_{\perp}(z) \end{pmatrix} \qquad (2.22) $$

where the parallel $D_{\parallel}(z)$ and perpendicular $D_{\perp}(z)$ components are inhomogeneous along *z* direction, and the components along *x* and *y* direction were verified to be similar $D_x = D_y \equiv D_{\parallel}$ [Honorio et al., 2022b].

The components of the tensor of diffusion can be computed using (equilibrium) molecular dynamics simulations, as detailed in the following.

* *Diffusion parallel to the walls in a slit pore.* Using Einstein's formula, the parallel components of the diffusion tensor can be computed by (e.g., [Boon and Yip, 1991]):

$$ D_{\parallel} = \frac{1}{4} \lim_{t \to \infty} \frac{\langle [r_x(t)]^2 + [r_y(t)]^2 \rangle}{t} \qquad (2.23) $$
where $r_i$ ($i = x, y$) are the components of the displacement vector. As for bulk fluids, corrections periodic boundary conditions artifacts are proposed in the literature [Simonnin et al., 2017].

* *Diffusion perpendicular to the walls in a slit pore.* The diffusion perpendicular to the wall in a slit pore involves collisions with the walls and the effects of the confinement potential $V_c(z)$ resulting in an overall non-Fickean dynamics. At short timescales, $t \ll H^2/D_\perp$, the dynamics of a single particle should be Fickean because the particle has not yet experienced any collision on the confining walls or has been subjected to the restoring force from the confinement potential [Burov et al., 2010]. At larger timescales, a subdiffusive regime takes place. The MSD along z-direction $\langle [r_z(t)]^2 \rangle$ can be linked to diffusion coefficient $D_\perp$ using a formalism based on Smoluchowski equation, in which the collisions and, in a simplified fashion, the confinement potential $V_c(z)$ (assumed as constant) are taken into account [Kusumi et al., 1993]:

$$ \langle [r_z(t)]^2 \rangle = \frac{H^2}{6} - \frac{16H^2}{\pi^4} \sum_{n=1(odd)}^{\infty} \frac{1}{n^4} \exp \left[ -D_{\perp} \left( \frac{n\pi}{H} \right)^2 t \right] \quad (2.24) $$

Within this formalism, the diffusion model is equivalent to a free Brownian diffusion within an infinitely high square well potential [Kusumi et al., 1993]. In agreement with what was discussed above, at times $t = 0$, the slope of the MSD profile is proportional to $D_\perp$, as in the Fickean case; then, the MSD goes asymptotically to $H^2/6$, as a manifestation of a subdiffusive regime [Kusumi et al., 1993]. Finite-size corrections have also been proposed for $D_\perp$ following the form of Eq. 2.21 [Botan et al., 2015].

The effect of the temperature in both parallel and perpendicular self-diffusion coefficients of water confined in C-S-H are shown in Figures 2.7 (A) and (C), respectively. The simulations were performed according to five temperatures. PBC corrections are applied to the parallel self-diffusion coefficients. The self-diffusion of water confined in C-S-H can be fitted by the expression [Honorio et al., 2022b]:

$$ D_i(d, T) = D_{ref} \left( 1 - \exp \left[ -\frac{d-20}{16.7} \right] \right) \exp \left[ -\frac{E_a}{R} \left( \frac{1}{T} - \frac{1}{300} \right) \right] \quad (2.25) $$

where $R$ is the gas constant, and $i$ corresponds to the parallel or perpendicular direction. For the parallel component, $D_{ref} = D_{bulk}$ is the bulk value of water self-diffusion for a reference temperature (here, $T=300$ K). For the perpendicular component, and since the dynamics are reduced in that direction, $D_{ref} \approx D_{bulk}/14$. The activation energy $E_a(d)$ depends on the pore-size, as shown in Figure 2.8, and can be fitted according to the expression (for $d$ in Å) [Honorio et al., 2022b]:

$$ E_a(d) = -A_1 (1 - \exp[-A_2(d - A_3)])^2 + A_4 \quad (2.26) $$
where the fitting parameter are ($A_1, A_2, A_3, A_4$) = (4.18 kJ/mol, 0.054 1/Å, 30.64 Å, 22.11 kJ/mol) for the parallel component, and ($A_1, A_2, A_3, A_4$) = (7.01 kJ/mol, 0.061 1/Å, 27.00 Å, 22.40 kJ/mol) for the perpendicular component.

<figure><img src="image_10.png" /><figcaption>**Figure 2.7:** Self-diffusion coefficient of water confined in C-S-H as a function of the interlayer spacing for various temperatures [Honorio et al., 2022b]: (A) parallel and (C) perpendicular components; (B) and (D) respective Log-Log plots showing a dynamic transition near $d \approx 26$ Å. Full lines in (A) and (B) correspond to fitting with Eq. 2.25.</figcaption></figure>

Unlike the diffusion coefficient, no significant anisotropy is observed in the activation energy. The activation energy exhibits a non-monotonous dependence on the pore size with a maximum at $d \approx 28$Å. This value in which the peak appears is close to the one observed for the dynamic transition observed in Figures 2.7 (B) and (D) showing that two regimes can be observed with $D_{xy}$ and $D_z$ scaling as $\propto d^\beta$ with $\beta = 12-13$ for $d \lesssim 26$ Å, and $\beta = 0.6-0.7$ for $d \gtrsim 26$ Å with a slight dependence of the system temperature. In Section 3.3, it is discussed that for interlayer distances $d \gtrsim 26$ Å, the radial distribution functions (RDF) of confined water are very similar to the one of bulk water [Honorio et al., 2021b]. For smaller basal spacings, though, a peakier RDF is observed, which indicates a significant change in water structure. Peakier RDF and slow dynamics [Gado et al., 2004] can be related to a transition from a liquid to an arrested (or glassy [Youssef et al., 2011]) state [Honorio et al., 2021b]. Eq. 2.25 shows that the temperature dependence of
the diffusion of confined water approximately follows an Arrhenius thermal activation. The associated activation energy exhibits a non-monotonous dependence on the pore size with a maximum at $d \approx 28$ Å. Note that the expressions above are only valid for $H > 9.5$ Å, not capturing, therefore, the strongly decelerated dynamics in smaller pores (for $H < 9.5$ Å the Eq. 2.25 yields non-physical negative self-diffusion coefficients).

<figure><img src="image_11.png" /><figcaption>**Figure 2.8:** Activation energy of the self-diffusion of water confined in C-S-H according to the [Honorio et al., 2022b] (A) parallel and (B) perpendicular self-diffusion coefficients.</figcaption></figure>

Further developments are needed to link self-diffusion to the corrected/collective diffusion, and macroscopic diffusion [Coasne, 2016].

### 2.4.2 Viscosity

The most usual way of computing the viscosity of fluids in a molecular simulation is from the auto-correlation of the (out-of-diagonal components) of pressure tensor based on the Green-Kubo formalism. For the shear viscosity [Medina et al., 2011, Allen and Tildesley, 1989]:

$$ \eta = \frac{1}{a_n} \sum_{\alpha,\beta} \lim_{t \to \infty} \frac{V}{k_B T} \int_{t=0}^{\infty} \langle P_{\alpha\beta}(t) P_{\alpha\beta}(0) \rangle dt \quad (2.27) $$
some authors consider that $P_{\alpha\beta}$ are only three out-of-diagonal components of pressure tensor $P_{xy}$, $P_{xz}$, and $P_{zy}$ and thus $a_n=3$ [Allen and Tildesley, 1989]; other authors also include the differences between the diagonal elements $(P_{xx} - P_{yy})/2$ and $(P_{yy} - P_{zz})/2$, thus $a_n=5$ [Medina et al., 2011, Allen and Tildesley, 1989]. For the bulk viscosity [Allen and Tildesley, 1989]:

$$ \eta_{bulk} = \frac{1}{3} \sum_{\alpha,\beta} \lim_{t \to \infty} \frac{V}{k_B T} \int_{t=0}^{\infty} \langle P_{\alpha\beta}(t) P_{\alpha\beta}(0) \rangle dt \quad (2.28) $$

where $P_{\alpha\beta}$ are the diagonal components of pressure tensor $P_{xx}$, $P_{yy}$, and $P_{zz}$. To get the properties of confined fluids, $V$ in the equations above should be the volume effectively occupied by the fluid, which can be computed using the sum of the Voronoï volume of water molecules. To define the upper integration bound in the equations above, the decay of the stress auto-correlations must be checked (for confined water, a choice of a few picoseconds is generally enough). Contrary to self-diffusion, the viscosity of SPC/E water is insensitive to simulation cell size [Yeh and Hummer, 2004].

Other strategies can be deployed to compute the shear viscosity of confined water [Honorio, 2022a]. Non-equilibrium Molecular Dynamics (NEMD) of Poiseuille flow gives velocity profiles (see next Section) that can be fitted to get the viscosity. Another strategy is an indirect calculation: if the self-diffusion coefficient of confined water is known, the shear viscosity $\eta$ can be calculated using Einstein-Stokes equation (e.g., [Batchelor and Batchelor, 2000]):

$$ \eta = \frac{k_B T}{6\pi D a} \quad (2.29) $$

where $a$ is a coefficient corresponding to the radius of an equivalent spherical fluid particle, and $D$ is the self-diffusion coefficient. The radius of a water molecule can be approximated by the van der Waals radius of oxygen so that $a = 1.52$ Å.

Based on a recent work [Honorio, 2022a], Figure 2.9(A) shows the shear viscosity of water confined in C-S-H using the three strategies. Direct computation (using auto-correlations) yields a non-negligible standard deviation. The three approaches are reasonably in agreement converging to the shear viscosity of SPC/E bulk water for pore sizes $H$ above approximately 3-4 nm. The values obtained for the viscosity in the larger pores are in close agreement with the values reported in the literature for SPC/E water $\eta =0.64-0.72$ cP(=mPa.s=10⁻³ Pa.s) (see [Tazi et al., 2012] for a collection of various references). Figure 2.9(B) compares the viscosity of water confined in C-S-H with data from the literature on other nanoporous materials (clays [Low, 1976] and silica [Haneveld et al., 2008, Persson et al., 2007, Thamdrup et al., 2007, Jarlgaard et al., 2006]). The results are normalized with respect to the corresponding shear viscosity of bulk water since SPC/E shows a smaller viscosity than experiments in water. The global evolution of viscosity of confined water in C-S-H with the pore size follows a similar trend to the other data.
<figure><img src="image_12.png" /></figure>

**Figure 2.9:** Shear viscosity of water confined in C-S-H as a function of the effective pore size *H* [Honorio, 2022a]: (A) Values computed for C-S-H from fitting NEMD velocity profiles, using Einstein-Stokes and (parallel to the pore wall) self-diffusion coefficient [Honorio et al., 2022b], and direct auto-correlation functions (Eq. 2.27); (B) Shear viscosity normalized with respect to the shear viscosity of bulk water ($\eta_{bulk}$, not to be confused with the bulk viscosity) comparison of results in (A) with experimental data from the literature on other nanoporous materials. The results of Low [Low, 1976] refer to interlayer distance in clays. The experimental studies [Haneveld et al., 2008, Persson et al., 2007, Thamdrup et al., 2007, Jarlgaard et al., 2006] refers to silica nanochannels. In both graphics, the value of the shear viscosity of bulk water is shown for comparison.

### 2.4.3 Permeability

Following usual hydrodynamics, the (intrinsic) permeability tensor $\kappa$ (with each component with dimension [$L^2$]) can be obtained from the average of the velocity:
$$ \bar{\mathbf{v}} = \frac{1}{H} \int_{-H/2}^{H_p/2} \mathbf{v} dz = \frac{\kappa}{\eta} \nabla \mathbf{p} \qquad (2.30) $$

where $H$ is the effective thickness of the pore, $z$ is the orthogonal distance from the pore central plane, and $\nabla \mathbf{p} = [\nabla p_x, \nabla p_y, 0]_{x,y,z}$ is the pressure gradient. Molecular simulations can be deployed to get velocity profiles in nanochannels resulting from pressure gradient-induced flows accounting for the relevant stick/slip boundary conditions [Bocquet and Barrat, 1993, Bocquet and Barrat, 1994].

The solution of the dynamics of a Newtonian fluid described by the Navier-Stokes equation in the case of laminar flow due to a pressure gradient is [Bazant and Vinogradova, 2008]:

$$ \mathbf{v}(z) = -\frac{H^2}{4\eta} \left\{ \frac{1}{2} \left[ 1 - \left( \frac{2z}{H} \right)^2 \right] \mathbf{I} + \mathbf{A}_{pf} + \left( \frac{2z}{H} \right) \mathbf{B}_{sf} \right\} \nabla \mathbf{p} \qquad (2.31) $$

Equation 2.31 reduces to a scalar form in the case of isotropic identical lower and upper boundary conditions [Předota et al., 2007]. The first term is the usual parabolic profile from Poiseuille flow in the direction of the pressure gradient; the second term is a slip-driven plug flow in the $\mathbf{A}_{pf} \nabla \mathbf{p}$ direction; and the third term is a linear shear flow in the $\mathbf{B}_{sf} \nabla \mathbf{p}$ direction, which vanishes when $b^+ = b^-$ [Bazant and Vinogradova, 2008]. The tensors $\mathbf{A}_{pf}$ and $\mathbf{B}_{sf}$ can be diagonalized when $\mathbf{b}^\pm$ is symmetric. The corresponding principal values are the pairs $A_s(b_1^+, b_1^-)$, $A_s(b_2^+, b_2^-)$ and $B_s(b_1^+, b_1^-)$, $B_s(b_2^+, b_2^-)$, respectively, with

$$ A_s(b^+, b^-) = \frac{b^+ + b^- + 4b^+b^-/H}{H + b^+ + b^-} \qquad (2.32) $$

$$ B_s(b^+, b^-) = \frac{b^+ - b^-}{H + b^+ + b^-} \qquad (2.33) $$

In the case of anisotropic textures, $\mathbf{A}_{pf}$ and $\mathbf{B}_{sf}$ are to be rotated appropriately.

Figure 2.10(A) depicts a scheme showing the (parabolic) velocity profile from a Poiseuille flow confined between two surfaces with different textures described by the boundary conditions $b^+$ (upper surfaces) and $b^-$ (lower surface). These values of the slip boundary can assume negative, null, or positive values according to the fluid-wall interactions (Figure 2.10(B)) [Zhan et al., 2020]. Hydrophilic surfaces generally lead to negative slip length due to the layers of adsorbed fluid. Hydrophobic surfaces lead to positive slip lengths. An intermediary midwetting surface could be found so that no slip is observed; only in this specific condition slip boundary conditions is not relevant.

In the case of C-S-H, the density and velocity profiles obtained for mesopores for a forced flow along in-plane directions are shown in Figure 2.11. The parabolic profiles characteristic of Poiseuille flow are obtained in the central portion of mesopores. The expression in Eq. 2.31 is used to fit the profiles using as independent fitting parameters the viscosity $\eta$, upper $b^+$, and lower $b^-$ slip lengths. The values obtained
<figure><img src="image_13.png" /><figcaption>**Figure 2.10:** (A) Scheme showing a (parabolic) velocity profile from a (Poiseuille) flow confined between two surfaces with arbitrary textures described by the boundary conditions $b^+$ (upper) and $b^-$ (lower) (adapted from [Bazant and Vinogradova, 2008]). Each texture can be defined in their principal directions (here noted, 1 and 2); the textures are aligned when each pair {$b_1^+, b_1^-$} and {$b_2^+, b_2^-$} is colinear. (B) Negative slip, no-slip, and positive slip boundary conditions associated with hydrophilic, midwetting, and hydrophobic surfaces (adapted from [Zhan et al., 2020]).</figcaption></figure>

for the viscosity in the larger pores are in close agreement with those reported in the literature for SPC/E water $\eta = 0.64-0.72$ cP [Tazi et al., 2012] (the pore-size dependence of water viscosity was discussed in Section 2.4.2).

In the case of the Poiseuille flow in Eq. 2.31, the integral in Eq. 2.30 yields [Bazant and Vinogradova, 2008]:

$$ \kappa = \frac{H^2}{12} \kappa^* = \frac{H^2}{12} (\mathbf{I} + 3\mathbf{A}_{pf}) \qquad (2.34) $$

where $\kappa^* = (\mathbf{I} + 3\mathbf{A}_{pf})$ is the dimensionless permeability. The expression above shows that the permeability increases with the slip-driven plug flow (according to the direction $\mathbf{A}_{pf} \nabla \mathbf{p}$) but is unaffected by the slip-driven shear flow (associated with $\mathbf{B}_{sf}$) [Bazant and Vinogradova, 2008]. For upper and lower aligned textures of the pore surfaces (even if the textures are different), the permeability tensor will also be aligned with the textures. In this case, the principal values of $\kappa^*$ [Bazant and Vinogradova, 2008] are $\kappa_s^*(b_1^+, b_1^-)$ and $\kappa_s^*(b_2^+, b_2^-)$ with

$$ \kappa_s^*(b^+, b^-) = \frac{H + 4(b^+ + b^-) + 12b^+b^-/H}{H + b^+ + b^-} \qquad (2.35) $$
<figure><img src="image_14.png" /></figure>

**Figure 2.11:** Profiles of (A) density and (B) velocity for mesopores in C-S-H at 300 K for a pressure gradient applied along the x-direction [Honorio, 2022a]. The vertical dashed and solid lines indicate the position for the (negative) slip boundary; the distance between two vertical solid lines is the effective pore size according to the density profile of water oxygens. Density profiles are shifted one unit upwards for readability. Fittings (solid lines following the symbols in velocity profiles) using Eq. 2.31 are provided for each velocity profile.

The intrinsic permeability at the channel scale is shown in Fig. 2.12(A), considering no-slip and slip boundary conditions. In the latter, bilinear (defining micro- and meso-pore domains) and exponential fits are considered; and the average of upper and lower slip lengths is used as input. Stick boundary condition leads to larger permeability than (negative) slip boundary conditions. The bilinear fits of the slip lengths lead to a singularity in the vicinity of $H_p=20$ Å. The exponential fits ensure continuity but with a non-monotonous evolution (a minimum is also observed in the vicinity of $H=20$ Å).
<figure><img src="image_15.png" /><figcaption>**Figure 2.12:** (A) Intrinsic permeability calculated using Eq. 2.34 as a function of the interlayer space in slit channels with stick and slip boundary conditions [Honorio, 2022a]. The bilinear (defining micro- and meso-pore domains) and exponential fits of the slip length are considered. (B) Corrected intrinsic permeability $\kappa^c = f_\eta f_H \kappa$ calculated in slit channels with stick and slip boundary conditions. Only the exponential fit of slip lengths is considered. The inset zooms in the small micropores range.</figcaption></figure>

Two correction factors are proposed for the intrinsic permeability: $\kappa^c = f_\eta f_H$, where $f_\eta = \eta_{bulk}/\eta(H)$ is the correction regarding the confinement dependence of the viscosity accounting for the decrease in viscosity with decreasing pore size; and $f_H = H/(H + d_{sol})$ is the correction regarding the volume occupied by the channel (in the layered configuration studied, the entire volume of the system is $(H+d_{sol})S_s$, where $S_s$ is the surface area parallel to the pore plane, and $d_{sol}$ is the thickness of the solid layer; for the system studied here $d_{sol}=4.11$ Å). Both corrections are only significant for small $H$. The corrected intrinsic permeability is shown in Fig. 2.12(B) considering stick and slip boundary conditions. The stick boundary conditions lead to larger corrected permeability than the slip boundary. The corrected permeability, considering the slip length, for a nanochannel thickness of 10 Å (comparable to the
hydraulic radius reported by Powers for C-S-H gel as 0.5 nm [Powers, 1958]) is on the order of 10⁻²² m², which is in agreement with the measurements on C-S-H.

## 2.5 Electromagnetic properties

Electromagnetic methods, particularly dielectric and electrical measurements, have been deployed to assess the water content of porous materials, including cement-based materials [Klysz and Balayssac, 2007, Rhim and Buyukozturk, 1998, Laurens et al., 2002, Cataldo et al., 2018]). The success of such methods relies on the high contrast between the electrical conductivity and dielectric permittivity of water solutions with ionic species and the solid phases.

### 2.5.1 Dielectric permittivity

**Bulk aqueous solutions.** The theoretical background for calculating the dielectric response of aqueous solution is presented in the sequel. The frequency-dependent dielectric permittivity $\epsilon(f)$ can be computed from the correlations of polarization at equilibrium [Köhler and Findenegg, 1972, Rinne et al., 2014a]:

$$ \epsilon(f) - 1 = \chi(f) = -\frac{1}{3Vk_BT\epsilon_0} \int_0^\infty e^{-2\pi i f t} \langle \mathbf{P}(0).\dot{\mathbf{P}}(t) \rangle dt \quad (2.36) $$

where $\dot{\mathbf{P}}(t)$ is the time-derivative of the total polarization $\mathbf{P}(t) = \sum_{i=1}^n q_i \mathbf{r}_i(t)$ of a system with $n$ particles $i$ with charge $q_i$, position $\mathbf{r}_i(t)$ (e.g. [Praprotnik and Janežič, 2005]). The total system polarization $\mathbf{P}$ is related to the electric field $\mathbf{E}$ via the complex frequency-dependent dielectric susceptibility $\chi(f) = \chi'(f) - i\chi''(f)$ by $\mathbf{P}(f) = \chi(f)\epsilon_0\mathbf{E}(f)$ where $\epsilon_0$ is the vacuum permittivity.

For a salt aqueous solution, the total polarization is the sum of the contributions of the water $\mathbf{P}_W$ and ionic $\mathbf{P}_I$ polarizations: $\mathbf{P} = \mathbf{P}_W + \mathbf{P}_I$. The ionic current $\mathbf{J}_I$ is related to the ionic polarization by $\mathbf{J}_I(t) = \dot{\mathbf{P}}_I(t)$, and to the electric field $\mathbf{E}$ via the frequency-dependent ionic conductivity $\sigma_e(f) = \sigma'_e(f) - i\sigma''_e(f)$ by $\mathbf{J}_I(f) = \sigma_e(f)\epsilon_0\mathbf{E}(f)$. The microscopic relaxation behavior is described by the correlations functions associated with water-water, ion-water, and ion-ion contributions, respectively [Rinne et al., 2014a, Rinne et al., 2014b]:

$$ \phi_W(t) = \frac{\langle \mathbf{P}_W(0).\mathbf{P}_W(t) \rangle}{3Vk_BT\epsilon_0} \quad (2.37) $$

$$ \phi_{IW}(t) = \frac{\langle \mathbf{P}_W(0).\mathbf{J}_I(t) - \mathbf{J}_I(0).\mathbf{P}_W(t) \rangle}{3Vk_BT\epsilon_0} \quad (2.38) $$

$$ \phi_I(t) = \frac{\langle \mathbf{J}_I(0).\mathbf{J}_I(t) \rangle}{3Vk_BT\epsilon_0} \quad (2.39) $$

From which the regularized susceptibility $\Delta\chi(f)$ can be decomposed into three separate contributions $\Delta\chi(f) = \chi_W(f) + \chi_{IW}(f) + \Delta\chi_I(f)$ with
<table><caption>**Table 2.2:** Cole-Cole fit parameters from MD simulations of pore solution and SPC/E water, and data from the literature on water. <sup>a</sup> [Honorio et al., 2020a]. <sup>b</sup> Buchner et al. [Buchner et al., 1999] experiments. <sup>c</sup> Rinne et al. [Rinne et al., 2014a] MD simulations.</caption><thead><tr><th></th><th>ε<sub>CC</sub></th><th>τ<sub>CC</sub></th><th>α</th></tr></thead><tbody><tr><td>PC1</td><td>73.6 ± 2.9</td><td>9.97 ± 0.97</td><td>0.000± 0.060</td></tr><tr><td>PC2</td><td>67.1 ± 2.6</td><td>10.26± 0.97</td><td>0.009 ± 0.060</td></tr><tr><td>PC3</td><td>73.9 ± 3.2</td><td>12.55 ± 1.33</td><td>0.000± 0.061</td></tr><tr><td>PC4</td><td>60.7 ± 2.7</td><td>9.33 ±1.14</td><td>0.000± 0.050</td></tr><tr><td>PC5</td><td>56.0 ± 2.1</td><td>9.80 ± 0.96</td><td>0.040 ± 0.056</td></tr><tr><td>PC6</td><td>55.2 ± 2.4</td><td>9.33 ± 1.27</td><td>0.098 ± 0.055</td></tr><tr><td>PC7</td><td>60.4 ± 3.0</td><td>10.30 ± 1.52</td><td>0.052 ± 0.063</td></tr><tr><td>SPC/E <sup>a</sup></td><td>68.2 ± 2.7</td><td>8.35 ± 0.94</td><td>0.033 ± 0.055</td></tr><tr><td>SPC/E <sup>b</sup></td><td>78.4</td><td>8.27</td><td>0.000</td></tr><tr><td>SPC/E <sup>c</sup></td><td>69.9</td><td>10.72</td><td>0.014</td></tr></tbody></table>

$$ \chi_W(f) = \phi_W(0) - i2\pi f \int_0^\infty e^{-2\pi i f t} \phi_W(t) dt \quad (2.40) $$

$$ \chi_{IW}(f) = -2 \int_{0}^{\infty} e^{-2\pi i f t} \phi_{IW}(t) dt \quad (2.41) $$

$$ \Delta\chi_I(f) = -\frac{i}{2\pi f} \int_0^\infty (e^{-2\pi i f t} - 1) \phi_I(t) dt \quad (2.42) $$

The dielectric response for the bulk aqueous solutions (Figure 2.6) mimicking the pore solutions in cement-based materials (as done for the self-diffusion coefficients [Honorio et al., 2019a]) was also determined. Most high-frequency electromagnetic remote sensing method is carried out in a frequency range of around 1 GHz. In this frequency range, the effect of free water (or bulk water) is predominant. The contribution of adsorbed water is expected to manifest in low frequency range [Miura et al., 1998, Hager and Domszy, 2004]. Thus, interesting information can be obtained from the analysis of bulk solutions representative of the pore solutions.

The dielectric spectra obtained for the pore solutions with MD simulation can be fitted using the Cole-Cole relaxation model [Cole and Cole, 1941]:

$$ \epsilon(f) = \frac{\epsilon_{CC} - \epsilon_{\infty}}{1 + (i2\pi f \tau_{CC})^{1-\alpha_{CC}}} + \epsilon_{\infty} \quad (2.43) $$

where $\epsilon_{CC}$ and $\tau_{CC}$ are, respectively, the Cole-Cole amplitude and characteristic relaxation time; and, $\alpha_{CC}$ is an exponent parameter ranging from -1 to 1 ($\alpha_{CC}=0$ corresponding to Debye model). An $\epsilon_{\infty} = 1$ is adopted since the force fields employed in MD simulation do not comprise atomic polarization. To fit simultaneously the real and imaginary parts of the dielectric spectrum, the function $\Delta\chi_i(f)\Delta\chi_i(f)$ (where
<figure><img src="image_16.png" /><figcaption>**Figure 2.13:** Frequency-dependent dielectric permittivity $\epsilon = \Delta\chi + 1$ of the pore solution [Honorio et al., 2020a]: comparison of the mean value computed from dielectric spectra of PC1-PC7 with the experimental data from the literature [Guihard et al., 2019] on real pore solutions. The real $\epsilon'$ and imaginary $\epsilon''$ parts are shown. The light blue zones depict the standard deviation associated with the MD results.</figcaption></figure>

$\Delta\chi_i(f)$ refers to the spectrum data to be fitted and $\bar{\epsilon}$ denotes the complex conjugate) is minimized. The parameters of Cole-Cole fits are gathered in Tab. 2.2. Simulation and experimental data on bulk water from the literature are also provided. The static dielectric permittivity of SPC/E water is known to be lower than the experimental value (e.g., [Zarzycki and Gilbert, 2019]). We obtain 69.0 for SPC/E water, which agrees with other values from the literature (e.g.[Smith and Dang, 1994] 67 ± 10).

Figure 2.13 compares the frequency-dependent dielectric permittivity $\epsilon = \Delta\chi + 1$ of the pore solution obtained from MD simulation with experimental data [Guihard et al., 2019] on synthesized pore solution based on the composition of a CEM I concrete. A good agreement is observed with respect to the imaginary part of the dielectric spectra. The difference between the experimental part can be attributed to the well-known deviation from the experimental value of the static dielectric constant SPC/E water (e.g., [Zarzycki and Gilbert, 2019]). In agreement with our Cole-Cole fits that yields low $\alpha$, Guihard et al. [Guihard et al., 2019] observe that Debye fits can describe reasonably well the dielectric spectra of the pore solution.

**Confined water.** To compute the dielectric response of confined water, the microscopic relaxation process described by Eq. 2.37 must be put in the tensorial form [Ballenegger and Hansen, 2005, Bonthuis et al., 2011]:

$$ \phi_{W,ij}(t) = \frac{\langle \mathbf{P}_{W,i}(0) \cdot \mathbf{P}_{W,j}(t) \rangle}{3V_{cav}k_B T \varepsilon_0} \qquad (2.44) $$
where i, j are the x, y or z directions in a Cartesian frame, and $V_{cav}$ is the volume of the cavity occupied by water. The water self-contribution to the susceptibility and dielectric response is computed as in the case of the bulk solution but accounting for the possible anisotropies due to the cavity geometry. As for the bulk solutions, the frequency-dependent components of the dielectric permittivity tensor are given as the Fourier transform of $\phi_{W,ij}(t)$ [Froltsov and Klapp, 2007].

In the case of a slit pore with an infinite slab of fluid confined in the z direction by two infinite layers separated by a distance L, the dielectric permittivity tensor can be written, as in the case of self-diffusion, as a function of the parallel $\epsilon_\parallel$ and perpendicular $\epsilon_\perp$ components:

$$ \epsilon(f, z) = \begin{pmatrix} \epsilon_\parallel(f, z) & 0 & 0 \\ 0 & \epsilon_\parallel(f, z) & 0 \\ 0 & 0 & \epsilon_\perp(f, z) \end{pmatrix} \qquad (2.45) $$

The frequency-dependent components of the dielectric permittivity tensor reads [Froltsov and Klapp, 2007]: $\epsilon_\parallel(f) = 1 + \chi_\parallel(f)$ and $\frac{1}{\epsilon_\perp(f)} = 1 - \chi_\perp(f)$.

In the case of C-S-H, only the static dielectric response of water confined in C-S-H, at the layer level, has been quantified in the literature using simulations [Masoumi et al., 2019]. The frequency-dependent complex response was provided for the first time in the results shown in Figure 2.14 for the in-plane components [Ait Hamadouche et al., 2022]. The perpendicular component is negligible due to the translational invariance. The water response for interlayer distance d larger than 5 nm is similar to that of bulk water. Significant changes in the dielectric response with the d are observed in smaller pores. Again, Cole-Cole [Cole and Cole, 1941] model fits reasonably well the response: the obtained relaxation time increases significantly with d decrease, and the $\alpha_{CC}$ show a non-monotonous evolution with d for values of d larger than 4 nm ($\alpha_{CC}$ converges to a positive value close to zero, in this case).

### 2.5.2 Electrical conductivity

**Bulk aqueous solutions.** The static electrical conductivity of an ionic system is given by the integral of the ionic current auto-correlation function (e.g., [Caillol et al., 1986]):

$$ \sigma_{PS} = \frac{1}{3Vk_BT} \int_0^\infty \langle \mathbf{J}_I(0) \cdot \mathbf{J}_I(t) \rangle dt \qquad (2.46) $$

where $\langle \cdot \rangle$ stands for statistical ensemble averages. The ionic current $\mathbf{J}$ auto-correlation function $\langle \mathbf{J}_I(0) \cdot \mathbf{J}_I(t) \rangle = \langle \sum_{i=1}^N q_i \mathbf{v}_i(0) \cdot \sum_{i=1}^N q_i \mathbf{v}_i(t) \rangle$ can be written in terms of velocity $\mathbf{v}$ auto- and cross-correlation functions. Using the VAC formula for the self-diffusion coefficient (Eq. 2.20), the total static electrical conductivity of the electrolyte can be decomposed into two terms:
<figure><img src="image_17.png" /><figcaption>**Figure 2.14:** Frequency-dependent dielectric permittivity of confined water in C-S-H for various interlayer distances *d*: (A), (B) real and (D), (E) imaginary parts [Ait Hamadouche et al., 2022]. The (isotropic) response of bulk water is shown for comparison (experimental values according to Buchner et al. [Buchner et al., 1999] and simulation values for SPC/E water according to Rinne et al. [Rinne et al., 2014a]). The microscopic response obtained from MD simulation was fitted with a simplex exponential decay (e.g., Debye model) and with a high degree exponentially decaying polynomial capturing local oscillations. The functions regarding each interlayer distance are shifted 50 units upwards for readability (the horizontal gray line marks the zero for each curve); bulk water values are shifted 500 units upwards. To quantify the (an)isotropy in the parallel components of the dielectric tensor, the difference $\epsilon_{yy} - \epsilon_{xx}$ is calculated for the (C) real and (D) imaginary parts of the dielectric response for selected slit pore sizes.</figcaption></figure>
$$ \sigma_{PS} = \sigma_{PS}^{NE} + \sigma_{PS}^{X} = \sigma_{PS}^{NE}(1 - \Delta) \quad (2.47) $$

where

$$ \sigma_{PS}^{NE} = \frac{e^2}{k_B T} \sum_{i}^{N} z_i^2 c_i D_i \quad (2.48) $$

$$ \sigma_{PS}^{X} = -\frac{e^2}{k_B T} \sum_{i}^{N} \sum_{j>i}^{N} z_i z_j c_i c_j D_{ij} \quad (2.49) $$

with $c_j$ being the molar concentration of ion $i$ whose charge number is $z_i$; and the coefficient $D_{ij} = \frac{1}{3} \int_0^\infty \langle \vec{v}_i(0) \cdot \vec{v}_j(t) \rangle$ accounts for the velocity correlations of ions $i$ and $j$, where the self-diffusion coefficient (Eq. 2.20) is retrieved for $D_{ii} = D_i$. Equation 2.48 is the Nernst-Einstein relation. The analysis of the ion pair-states in the electrolyte provides critical information regarding the relative contribution of $\sigma_{PS}^{NE}$ and $\sigma_{PS}^{X}$ terms to the overall conductivity. Figure 2.15 shows a schematic representation of the main ion pair states in aqueous solutions (e.g., [Rinne et al., 2014a]): Contact Ion-Pairs (CIP), single Solvent-separated Ion Pairs (SIP), Doubly Solvent-separated Ion Pairs (DSIP), and Free Ions (FI). CIP is ion pairs in contact: no water molecule is in-between CIP. SIP and DSIP are ion pairs in which one and two, respectively, water molecules are found in between the ions. If three or more water molecules separate a given ion pair, these ions are considered as FI. CIP exhibits correlated velocities, which translate in a not vanishing $D_{ij}$, whereas the ion velocity auto-correlation of other ion pairs, notably FI, tends to 0. CIP involving ions with opposite charges $a+$ and $a-$ does not carry a net charge, therefore, not contributing to overall Direct Current (DC) conductivity. In molecular dynamics computations of the ion self-diffusion coefficient $D_i$, these CIP that are not carrying charge are accounted for, though. The Nernst-Einstein relation describes the conductivity of an electrolyte works when ions cannot interact (for example, FI in a diluted system). In a system where the ions can interact (CIP, SIP, DSIP), the deviation parameter $\Delta$ is always positive (apart from the possibility of super-ionicity) [Harris and Kanakubo, 2015]. Consequently, the Nernst-Einstein relation overestimates the conductivity.

In the study of bulk solution mimicking pore solutions in cement systems, ion pair states were assessed using radial distribution functions [Honorio et al., 2019a, Honorio et al., 2020a]. The simulations show that cation-sulfates pairs are mostly CIPs, except for sodium-sulfate pair that is a DSIP at very early-ages pore solution (age <1 day) [Honorio et al., 2020a]. These CIP involve single-charged cations ($K^+$, $Na^+$) in contact with (doubly charged) sulfates ($SO_4^{2-}$). Such CIPs carry a charge and therefore contribute to the DC-conductivity. In the simulations, only a small fraction of cation-hydroxide pairs is a CIP [Honorio et al., 2020a]. When the total fraction of CIPs is small, and these occur between ions with dissimilar charges, then the main contribution in Eq. 2.47 is expected to come from the Nernst-Einstein term (Eq. 2.48). The deviations from the Nernst-Einstein relation are expected to be still
less significant in pore solution from cement systems with larger w/c due to a dilution effect. Indeed, previous studies report that the Nernst-Einstein equation seems sufficiently accurate to capture the experimentally measured electrical properties of concrete [Lu, 1997, Garboczi and Bentz, 1996, Snyder and Marchand, 2001], and other geomaterials [Revil, 2013]. Following the Nernst-Einstein relation but using as input self-diffusion computed from MD simulations, it was shown that the electrical conductivity is proportional not to the ionic strength $I = \frac{1}{2} \sum_{j=1}^{n} c_j z_j^2$ but instead to the factor $z_i^2 c_i D_i(c_i)$, in which the self-diffusion coefficients of the ions appear.

<figure><img src="image_18.png" /><figcaption>**Figure 2.15:** (a) Schematic representation of Contact Ion-Pairs (CIP), single Solvent-separated Ion Pairs (SIP), Doubly Solvent-separated Ion Pairs (DSIP), and Free Ions (FI), where a and b are the cation and anion charges, respectively [Honorio et al., 2020c]. (b) Atomic configuration of a bulk solution mimicking the pore solution of a cement paste at the very early-age in references [Honorio et al., 2019a, Honorio et al., 2020a]: water molecules are depicted as transparent sticks. CIPs observed among potassium and sulfates, and sodium and hydroxide are identified.</figcaption></figure>

As for the dielectric response, the frequency-dependent complex electrical conductivity can also be calculated from molecular simulations. The frequency-dependent ionic conductivity can be computed using the cross-correlations [Caillol et al., 1986, Rinne et al., 2014a]:

$$ \sigma_e(f) = \frac{1}{3Vk_BT} \int_0^\infty e^{-2\pi i ft} \langle \mathbf{J}_I(0).\dot{\mathbf{P}}(t) \rangle dt \quad (2.50) $$

The frequency-dependent ionic conductivity can be decomposed into the two terms $\sigma_e(f) = \sigma_{IW}(f) + \sigma_I(f)$ with

$$ \sigma_{IW}(f) = -i2\pi f \epsilon_0 \int_0^\infty e^{-2\pi i ft} \phi_{IW}(t) dt \quad (2.51) $$

$$ \sigma_I(f) = \epsilon_0 \int_0^\infty e^{-2\pi i ft} \phi_I(t) dt \quad (2.52) $$
Therefore, the static conductivity can be computed via $\sigma_e(f = 0) = \sigma_I(f = 0) = \epsilon_0 \int_0^\infty \phi_I(t)dt$.

The anisotropic electrical conductivity of confined electrolytes (e.g., in C-S-H) can also be computed. This aspect is yet to be worked out.

## 2.6 Coupling properties

### 2.6.1 Thermoelasticity: coefficient of thermal expansion

Finite difference and fluctuation approaches can be used to calculate the full thermal expansion tensor:

* *Finite differences* is performed based on the definition in an anisotropic framework (.e.g [Cagin et al., 1992, Ray, 1983]):

$$ \alpha_{ij} = \left. \frac{\partial \epsilon_{ij}}{\partial T} \right|_{\sigma_{ij}} \qquad (2.53) $$

* *Fluctuation formula.* The thermal expansion tensor is obtained from the cross fluctuations of the enthalpy $H_p$, and the deformation tensor in the N$\sigma$T ensemble [Cagin et al., 1992, Ray, 1983, Honório et al., 2019]:

$$ \alpha_{ij} = \frac{\langle H_p \epsilon_{ij} \rangle_{N\sigma T} - \langle H_p \rangle_{N\sigma T} \langle \epsilon_{ij} \rangle_{N\sigma T}}{k_B T^2}. \qquad (2.54) $$

In my works, the thermal expansion tensor has been computed for the following phases in cement systems using the fluctuation approach:

* **Clinker minerals.** The thermal expansion tensor was calculated for tricalcium silicate [Barbosa and Honorio, 2022]. The value at the principal basis reveals a significant anisotropy regarding the thermal expansion, with the expansion along the first principal direction being roughly twice the value in the third principal direction. The triclinic polymorph presents a smaller thermal expansion than the monoclinic polymorphs [Claverie et al., 2021, Qomi et al., 2015].

* **C-(A)-S-H.** For the first time, the full thermal expansion tensor has been computed for C-A-S-H [Hawchar and Honorio, 2022]. The thermal expansion of C-A-S-H shows a significant anisotropy, as expected in layered phases. The average volumetric thermal expansion $\alpha_v/3 = \text{Tr}(\alpha_{ij})/3$, is similar to the values reported for C-S-H in the literature including inverse analysis with micromechanics [Ghabezloo, 2011] and molecular simulations [Qomi et al., 2015, Krishnan et al., 2016].
* **AF-phases.** The thermal expansion tensors have been computed for ettringite, metaettringite [Honorio et al., 2020e, Honorio, 2022b], and Friedel's salt [Honorio et al., 2022a] using AFFF. As for the elastic constants, the simulations were performed for various water contents for ettringite, and monosulfoaluminate [Honorio et al., 2020e]. In the cases of metaettringite and Friedel's salt, the molecular simulations were the first estimates of this property in the literature. In the case of Friedel's salt, the diagonal components (and $\alpha_v$) peak at around 35°C, which can be associated with the LT→HT transition (displacive phase transitions involving symmetry changes can be accompanied by significant changes in the coefficients of thermal expansion or discontinuities in the case of second-order transitions [Testardi, 1975]).

* **ASR products.** Full thermal expansion tensors were computed, also for the first time, for crystalline ASR products based on shlykovite structure [Honorio et al., 2020f]. This property also exhibited ion-specificity, in particular in the component of the thermal expansion along a-direction (the coefficient decreases with the ion size: $\alpha_a$ followed the order Li > Na > K > Cs).

## 2.6.2 Poroelasticity: Biot coefficient

The apparent Biot coefficient $b$ quantifies how pore pressure $P_w$ is transferred to the solid skeleton in a porous medium (e.g., [Honorio et al., 2017a]):

$$b_i = \left. \frac{\partial \sigma}{\partial P_w} \right|_V \qquad (2.55)$$

where $P_w$ is the water pressure. Whenever the concept of Terzaghi effective stress holds, the apparent Biot coefficient is strictly 1. Poromechanics' usual interpretation [Coussy, 2010] posits that Terzaghi's effective stress principle holds whenever the solid skeleton is poorly compressible and the pores are saturated. At the nanoscale, even for rigid solid layers and saturated pores, apparent Biot coefficients may differ from 1. This result was observed in tobermorite [Honorio, 2019] and swelling clays [Honorio et al., 2017a], indicating that a poromechanics framework accounting for nanoscale effects is required to explain the behavior of these materials.

## 2.6.3 Drying shrinkage

Volume changes at the molecular scale due to RH decrease were evaluated for ettringite [Honorio et al., 2021a], and C-S-H [Masara et al., 2023].

**Ettringite.** The tensor of (crystal) shrinkage $\epsilon^{dry}(RH)$ was computed from the lattice parameters obtained from the GCMC-N$\sigma$T simulations (the GCMC stage enacts water exchange with an infinite reservoir of water with fixed temperature and chemical potential, which can be directly related to the RH; the N$\sigma$T stage allows for the changes in simulation box dimensions). Figure 2.16(A)-(D) shows the changes in lattice parameters and volume of ettringite due to desorption. Simulations show
a good agreement with the experimental data of [Baquerizo et al., 2016, Zhou et al., 2004].

According to the hexagonal crystal symmetry of ettringite, $\epsilon^{sh.}(RH)$ has three independent components in the form:

$$ \epsilon^{dry}(RH) = \begin{pmatrix} \epsilon_{xx}(RH) & \epsilon_{xy}(RH) & 0 \\ \epsilon_{xy}(RH) & \epsilon_{yy}(RH) & 0 \\ 0 & 0 & \epsilon_{zz}(RH) \end{pmatrix} \qquad (2.56) $$

where the components $\epsilon_{xz} = \epsilon_{yz} = 0$ and (due to the symmetry) $\epsilon_{xy} = \epsilon_{yx}$. The hypothesis of infinitesimal strain is incompatible with the significant dimension changes in ettringite. Therefore, the analysis should adopt the Lagrangian finite strain tensor (Eq. 2.7, which accounts for large strains).

<figure><img src="image_19.png" /><figcaption>**Figure 2.16:** Anisotropic volume changes of ettringite under desorption: comparison of (A)-(C) lattice parameters (a, b and c) and (D) volume changes obtained from simulations against the experimental data of Baquerizo et al. [Baquerizo et al., 2016], Zhou et al. [Zhou et al., 2004], and Gard (apud [Zhou et al., 2004]). (E) Density changes as a function of the water content comparison against Baquerizo et al. [Baquerizo et al., 2016] and the trend (dashed line) proposed by same authors: $\rho = 2.600 - 0.02925(n + 6)$ in g/cm$^3$ (adapted to the notation used here).</figcaption></figure>

The components of the Lagrangian tensor of crystal shrinkage $\epsilon^{dry}$ of ettringite are shown in Figure 2.17 (A) and (C) as a function of the RH and the water content, respectively. The fit of the lattice parameters from Figure 2.16 is used to compute the deformations, as detailed in [Honorio et al., 2021a]. In the ettringite domain, the infinitesimal strain hypothesis approximates the results obtained as shown in Figure 2.17 (B) and (D) as a function of the RH and water content, respectively. The components of the Lagrangian tensor of crystal shrinkage of ettringite are proportional to the water content.
<figure><img src="image_20.png" /><figcaption>**Figure 2.17:** Components of strain tensor obtained from simulations using the data from Fig. 2.16 as a function of (A) the RH (with (B) zoom in the small strain domain), and (C) the water content (with (D) zoom in the small strain domain) [Honorio et al., 2021a]. Linear fitting of the small strain domain is provided for both cases.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 2.18:** Changes in the basal spacing of C-S-H quantified in GCMC-NVT simulations [Masara et al., 2023] compared against experimental data from the literature for C-S-H with different Ca/Si ratios during desorption [Gutteridge and Parrott, 1976, Smith and Bayliss, 1972, Bayliss, 1973].</figcaption></figure>

**C-S-H.** GCMC-NVT simulations were performed to get the effective interactions at various RH (see Section 3.1). The results of these simulations provide a way of
quantifying the changes in basal spacing due to the RH variations. Simulation results display a larger basal spacing compared to experiments, which is associated with the presence of Ca(OH)₂ clusters (micoportlandite) in the interlayer C-S-H space of the C-S-H model adopted [Kunhi Mohamed et al., 2018]. These clusters are known to lead to larger basal spacings [Abdolhosseini Qomi et al., 2021] when compared with systems with lower interlayer calcium content. The rate of decrease of the stable basal distance in our simulation is less pronounced than in the experimental data. The effects of layer deformability (neglected in our simulations), might be one of the aspects at the origin of this difference. Ongoing GCMC-NσT simulations on C-S-H (as done for ettringite) can shed light on the effects of layer deformability.

## 2.7 Acoustic response

The acoustic behavior at the molecular level may be assessed using the thermoelastic properties as detailed in the following. The volumetric sound velocity [Bhowmick and Shenoy, 2006]:

$$S_v = \frac{1}{\sqrt[3]{\frac{1}{3}\left(\frac{1}{S_l^3} + \frac{2}{S_t^3}\right)}} \quad (2.57)$$

is computed from the average of longitudinal sound velocity (according to each Cartesian axial direction) $S_l = \sqrt{\frac{C_{11}+C_{22}+C_{33}}{3\rho}}$ and the transverse sound velocity $S_t = \sqrt{\frac{C_{44}}{\rho}}$, both obtained from the elastic constants moduli and the density $\rho$.

The mean free path of phonons $l_m$ can be estimated in the framework of kinetic theory using the thermal properties and acoustics parameters [Ziman, 2001]:

$$l_m = \frac{3\lambda_v}{\rho C_V S_v} \quad (2.58)$$

In the same framework, the relaxation time of phonons $\tau_m$ is given by [Bhowmick and Shenoy, 2006]:

$$\tau_p = \frac{3\lambda_v}{\rho C_V S_v^2} \quad (2.59)$$

These parameters have been computed for triclinic tricalcium silicate at ambient conditions (300 K) for the first time, as reported in Table 2.3. The values can be compared with the ones computed for monoclinic M3 C₃S [Qomi et al., 2015].

# 3 Sorption processes and behavior of adsorbed water

Nanoscale processes are at the heart of sorption in microporous adsorbing materials such as C-S-H, AF-phases, ASR products, etc. Understanding these processes
<table><caption>**Table 2.3:** Volumetric $S_v$, transverse $S_t$ and longitudinal $S_l$ acoustic velocities; and the mean free path $l_m$ and relaxation time $\tau_p$ of phonons of triclinic tricalcium silicate at ambient conditions (300 K) [Barbosa and Honorio, 2022]. <sup>a</sup> Computed here from the thermoelastic properties presented in [Qomi et al., 2015].</caption><thead><tr><th></th><th>This work</th><th>M3 <sup>a</sup></th></tr></thead><tbody><tr><td>$S_v$ (km/s)</td><td>5.5</td><td>4.4</td></tr><tr><td>$S_t$ (km/s)</td><td>5.0</td><td>3.9</td></tr><tr><td>$S_l$ (km/s)</td><td>8.36</td><td>7.1</td></tr><tr><td>$l_m$ (Å)</td><td>4.4</td><td>9.5</td></tr><tr><td>$\tau_m$ (fs)</td><td>81</td><td>216</td></tr></tbody></table>

requires a description at the molecular scale. Recent advances in understanding these processes have been obtained with molecular simulations [Abdolhosseini Qomi et al., 2021]. However, important questions remain unanswered (or, at least, poorly understood) regarding the behavior of adsorbed water (surface water or confined water) under relative humidity and temperature changes. This section deals with the recent results shedding light on the molecular processes associated with sorption and the behavior of adsorbed water in cement systems applications.

## 3.1 Disjoining pressure and effective interactions

The disjoining pressure can be directly computed from GCMC simulations that reproduce osmotic equilibrium conditions between the porous phase and an infinite fluid reservoir. As detailed below, the effective interactions between layers can be calculated from disjoining pressure isotherms.

In porous materials, the thermodynamic potential associated with drained experimental conditions is written as a function of the Helmholtz free energy $F$, the number of the fluid molecules $N_w$, and the fluid's chemical potential $\mu$ (e.g., [Honorio et al., 2017b]):

$$ \Lambda = F - \mu N_w \qquad (2.60) $$

This potential gives the "effective interactions" or "potential of mean force" (PMF) between two nanolayers under the osmotic equilibrium.

Given a system with a controlled interlayer distance $d$, $\Lambda$ reads [Honorio et al., 2021b, Honorio et al., 2017b]:

$$ \Lambda(d) = \Lambda_0(d_0, \mu, T) - S_s \int_{d_0}^{d} P \, d(d) \qquad (2.61) $$

where $S_s$ is the surface area (constant in our case) and $P$ is the confining pressure.

The confining pressure $P$ is associated with the thermodynamic potential $\Lambda$ through [Honorio et al., 2021b, Honorio et al., 2017b]:
$$P = \frac{1}{S_s} \left| \frac{\partial \Lambda}{\partial d} \right|_{\mu,T} \qquad (2.62)$$

The confining pressure can be directly calculated from the usual virial formula along the z- direction. The disjoining pressure Π is the pressure difference between the pore and the external bulk water reservoir: Π=P-P<sub>bulk</sub>.

In adsorbing layered materials, typical confining pressure isotherms are oscillating curves with decaying amplitudes when plotted as a function of the interlayer distance [Brochard et al., 2017]. The dehydrated state corresponds to the first branch of the curve in which the confining pressure strongly increases with the decrease of d due to steric repulsion between the solid layers as well as the counterions in the interlayer space. The first oscillation corresponds to the monohydrated state; the subsequent i oscillations correspond to the i-hydrated states up to the pore water state. Accordingly, the swelling energies also exhibit a multi-well configuration when plotted versus basal spacing. This is a critical element of hydration transitions. The (meta)stable domains, stable basal spacings, and discrete water layering in C-S-H nanopores can be identified from the confining pressures and PMF profiles.

Disjoining pressure and effective interaction profiles also inform on the (e.g., [Abdolhosseini Qomi et al., 2021]):

* **Cohesive pressure $P_{co}$** =min Π(d): the pressure needed to disjoin two solid layers, which corresponds to the global minimum in the pressure isotherms.

* **Surface energy $\gamma_s$** = -1/(2S<sub>s</sub>) ∫<sub>d<sub>eq</sub></sub><sup>∞</sup> Pd(d): the energy required to separate two layers placing them infinitely apart.

* **Out-of-pore-plane elastic modulus $E_{zz}$** can be computed from the slope in the pressure isotherms at a given interlayer distance $d_i$: $E_{zz} = -d_i \frac{\partial P}{\partial d}|_{d_i}$, where we have choose $d_i = d_{min}$.

These properties can be compared to experimental data [Abdolhosseini Qomi et al., 2021].

In my work, disjoining pressures and effective interaction were calculated for the systems below:

* **Swelling clays.** Disjoining pressure and effective interaction profiles of Na- and Ca-montmorillonite were provided [Honorio et al., 2017b, Honorio and Brochard, 2022]. The stable interlayer distances (minimum in the PMFs) agree with experimental data and previous molecular simulations [Smith et al., 2006, Tambach et al., 2006, Karaborni et al., 1996]. The analysis of disjoining pressure and PMF profiles enables explaining the coexistence of hydration states and the stability of clays at the scale of the stack of layers, giving new insights to explain the thermo-poro-mechanical behavior of clays bottom-up [Honorio et al., 2017b, Brochard et al., 2017].
<figure><img src="image_22.png" /><figcaption>**Figure 2.19:** Disjoining pressure isotherms (top) and free energy profiles (bottom) as a function of the interlayer distance for various RH at 300 K during adsorption and desorption [Masara et al., 2023].</figcaption></figure>

* **C-S-H.** The effects of temperature, chemical potential, and pressure on the disjoining pressure and effective interactions were studied for tobermorite [Honorio, 2019]. The disjoining pressure isotherms of tobermorite are negligibly affected by changes in water pressure under liquid saturated conditions. Using a realistic C-S-H model, the disjoining pressure and effective interactions exhibited a strong RH and pore size dependence [Masara et al., 2023], as shown in Figure 2.19. In C-S-H, the total disjoining pressure is mostly "joining" than disjoining, and Ca counterions play an important role in that, as shown in previous works [Gmira et al., 2004; Pellenq et al., 2008; Honorio
et al., 2021b]; calcium counterions are the primary source of cohesion in C-S-H while water has a disjoining behavior.

<figure><img src="image_23.png" /><figcaption>**Figure 2.20:** Crystalline ASR products (Na- and K-shlykovite) controlled on the RH [Honorio et al., 2020d]: (a)-(b) water content, (c)-(d) disjoining pressure Π and (e)-(f) free energy Λ profiles as a function of the basal spacing for Na- and K-shlykovite (left and right, respectively) at 296 K and under an imposed water pressure $P_W(\mu_W)$ of 1 atm. The insets in (c)-(f) zoom in on the oscillatory portions of the pressure and energy profiles.</figcaption></figure>

* **ASR products.** Disjoining pressure and effective interaction profiles were calculated for Na- and K-shlykovite [Honorio et al., 2020d]. Figure 2.20 shows the water content, disjoining pressure Π and free energy Λ profiles as a function of the interlayer distance for Na- and K-shlykovite under ambient conditions. The equilibrium basal distance in K-shlykovite is larger than in Na-shlykovite, which is expected since the radius of potassium is larger than sodium's. The cohesive pressures are similar for both cases. The surface energy for Na-
shlykovite is slightly larger than for K-shlykovite. As detailed in the next chapter, the PMF profiles are used as input in coarse-grained simulations.

## 3.2 Hysteresis, cavitation and capillary condensation

Hysteresis under sorption is explained as resulting from the metastability of the adsorbed fluid during adsorption, desorption, or both [Abdolhosseini Qomi et al., 2021]. Both desorption mechanisms (pore blocking and cavitation) and capillary condensation are associated with metastable transitions [Coasne et al., 2013, Pellenq et al., 2009]. Sorption hysteresis loop disappears in small nanopores, and capillary condensation becomes continuous and reversible at a critical temperature [Coasne et al., 2013]. Critical temperatures and pore sizes can be defined for the disappearance of hysteresis ($T_{ch}$ and $H_{ch}$), and for the reversibility/continuity of capillary condensation transition ($T_{cc}$ and $H_{cc}$). For a temperature $T$ and a pore size $H$ in the ranges $T_{ch} < T < T_{cc}$, $H_{cc} < H < H_{ch}$ capillary condensation transition occurs reversibly without hysteresis; whereas in the ranges at $T_{bulk} < T < T_{ch}$, $H_{ch} < H$, capillary condensation is hysteretic [Neimark et al., 2003].

Desorption in nanoporous materials occurs according to *pore blocking* or *cavitation* [Coasne et al., 2013]. With pore blocking, the pores empty when the constrictions isolate them from the neighboring bulk liquid phase. Pore blocking occurs if the liquid in the pore resists the depression caused by drying until the liquid in the pore evaporates. With cavitation, a gas bubble nucleates (and grows) within the nanopore while the constrictions remain saturated with liquid. Cavitation occurs if the liquid in the pore cannot resist the depression caused by drying, which leads to spontaneous bubble nucleation. The occurrence of both mechanisms concomitantly in the same pore material has been reported [Kikkinides et al., 2003, Reichenbach et al., 2011]. Cavitation is pointed out as a critical phenomenon in the drying of cement-based materials [Maruyama et al., 2018, Rastogi et al., 2022]. Cavitation can be either homogeneous or heterogeneous, and both cases seem relevant for cement-based materials [Rastogi et al., 2022].

Despite important developments in the last decade deploying molecular simulation techniques [Bonnaud et al., 2012, Pinson et al., 2015, Zhou et al., 2019], sorption hysteresis in microporous phases in cement systems is yet to be fully understood.

In my works, hysteresis under sorption was studied for the following systems:

*   **C-S-H.** GCMC simulations were performed for various RH under ambient temperature [Masara et al., 2023]. Figure 2.21 (A) and (B) show the desorption and adsorption isotherms, respectively, of the slit pores. For C-S-H interlayer pores of widths $d \ge 20.3$ Å, the desorption isotherms are continuous until a rapid decrease in the water content is observed. This steep decrease in the moisture content at a given RH has been discussed in the literature as being usually associated with a single-step process cavitation in slit-shaped pores [Sarkisov and Monson, 2001, Vishnyakov and Neimark, 2003, Coasne
and Pellenq, 2004]. For pores > 21.3 Å, cavitation occurs at an RH higher
than that of bulk water (i.e., at RH=60 % for the largest pore of 25.3 Å
interlayer distance). This observation might be attributed to the existence of
hydrophobic sites [Rasmussen et al., 2012] in an overall hydrophilic surface
(as is the case of C-S-H). The critical pore sizes associated with hysteresis
loop disappearance ($H_{ch} \approx 21.3$ Å) and reversibility of capillary condensation
transition ($H_{cc} \approx 19.3$ Å) were identified for the C-S-H for the first time for
ambient temperatures. The value obtained for $H_{cc}$ is on the order of the value
estimated from theoretical considerations [Evans et al., 1986, Coasne, 2016].

*   **ASR products.** GCMC-NσT simulations show that the water content and volume of Na and K-shlykovite remain stable for RH ranging from 20 to 100% [Honorio et al., 2020d]. The water content on shlykovite does not change significantly in a system under RH control, which corroborates that crystalline ASR products do not swell in the same way that swelling clays [Shi et al., 2019]. No hysteresis was observed in Na and K-shlykovite with the simulations.

*   **Ettringite.** Sorption isotherms obtained from GCMC-NσT simulations are shown in Figure 2.22 in comparison with the experimental data [Skoblinskaya and Krasilnikov, 1975]. The simulated desorption branch agree with the experimental data. For the cycles with rewetting from RH of 30 and 1%, close loop reversible hysteresis is observed by both experiments and simulations. However, for the cycle with rewetting from 0.1% RH, the reversibility of the adsorption branch at high RH (> 60%) is not captured by the simulations. The direct simulations in the osmotic ensemble on microporous phases exhibiting significant volume changes under sorption generally cannot capture experimental reversibility at high RH (a review of simulations in the osmotic ensemble and associated difficulties can be found in ref. [Coudert et al., 2011]). This observation was reported for various microporous materials such as MOFs (i.e., metal-organic frameworks, see [Bousquet et al., 2012, Ghoufi et al., 2012]). Therefore, provided the significant volume changes associated with ettringite desorption, it is expected that the reversibility in the adsorption cycle is not observed in the simulations at high RH, as is the case in the experiments. Other better-suited simulation techniques (e.g., Wang-Landau free energy [Bousquet et al., 2012] or using force fields fine-tuned to impose the flexibility of the crystal framework [Yang et al., 2011]) have been proposed to capture the reversibility at high RH. Upon drying, H-bonds between OH groups from neighboring columns and between OH groups and sulfates are formed, making it difficult to water molecules ingress in the structure upon rewetting. In the ettringite domain, H-bonds breaking and formation are reversible and follow the same path as desorption. Following the discussion in Chen et al. [Chen et al., 2018], hysteresis under sorption originates from a complex coupling between the fluid grand free energy and the microporous solid energy landscapes. Even if H-bonds are not the origin *per se* of sorption hysteresis, H-bonds are directly
<figure><img src="image_24.png" /><figcaption>**Figure 2.21:** (A) Water desorption isotherms for C-S-H basal distances between 13.7 and 25.3 Å and bulk water versus RH. The 13.7 Å basal distance pore is the equilibrium interlayer pore. The water content at a given RH ($N_w$) is normalized to the number of water molecules at 100 % RH (N). (B) Adsorption isotherms for C-S-H interlayer distances between 13.7 and 25.3 Å. (C) Sorption isotherms for the same range of C-S-H interlayer spacings and bulk water. The graphs are moved upward for clarity [Masara et al., 2023].</figcaption></figure>
<figure><img src="image_25.png" /></figure>

**Figure 2.22:** Sorption isotherms of ettringite obtained from GCMC-NσT simulations [Honorio et al., 2021a] compared against experimental data [Skoblinskaya and Krasilnikov, 1975] ($T = 291$ K): (A) Desorption branch compared to the experimental sets "A"-"E" as labeled in [Skoblinskaya and Krasilnikov, 1975]. (B) and (C) Sorption cycles with close loop reversible hysteresis from 30 and 1 % RH (dark yellow labels), respectively. (D) Sorption cycle with open-loop hysteresis: simulation cannot capture the reversibility under adsorption. The simulated desorption and adsorption branches are fitted with Langmuir isotherms. The standard deviations of water contents obtained from simulations are on the order of the symbol size and are not shown for conciseness.

related to the fluid grand free energy, being therefore an important ingredient to understand hysteresis under sorption.

## 3.3 Structure of water in microporous hydrates

Confinement and surface adsorption induces changes in the water structure, which is expected to result in different physical properties. In my works, structural changes of water were analyzed for the systems below.

* **C-S-H.** Water structural changes are quantified using the radial distribution function in Figure 2.23. For systems with interlayer spacing as large as approximately 3-5 nm, the structure of confined water retrieves the one observed in bulk conditions. The first peak in Ow-Ow RDF is slightly shifted to the left and is narrower with the confinement. Figure 2.23(C) shows the position of the first peak (the argument of $g_{Ow-Ow}(r)$ when it reaches its maximum) in
Ow-Ow RDF as a function of the c-length (here, the c-length is the interlayer
distance plus 13.7 Å). From c-lengths above 4 nm, the position of the peak is
> 2.7 Å, close to the values corresponding to bulk water. For smaller c-lengths,
a decrease in this peak position shows that the confined water densifies with
the confinement. Densification with decreasing pore sizes is also observed in
the Hw-Hw RDFs. Regarding *H-bond statistics*, the total H-bonds *HB* are
computed as the sum of the contribution of water oxygen *HB*^WW, hydroxide
groups *HB*^HH, water oxygen and hydroxide groups *HB*^WH, layer oxygen and
water oxygen *HB*^LW, and layer oxygen and hydroxide *HB*^LH [Honorio et al.,
2021b]. The *HB*^WW retrieves the value associated with bulk water for inter-
layer distances exceeding about 5 nm. A larger number of H-bonds is observed
when one accounts for the contributions of oxygens out of water molecules, and
this contribution is more significant the smaller the pore. This result indicates
that confined fluids are more constrained in what concerns H-bonds than the
bulk system. More constraints can explain water's higher bulk modulus in
ultra-confined systems.

<figure><img src="image_26.png" /></figure>

**Figure 2.23:** Radial distribution function of (A) Ow-Ow and (B) Hw-Hw pairs of water confined in C-S-H compared to bulk water [Honorio et al., 2021b]. The peaks indicating the nearest, second (first peak), second, third, etc., nearest correlations are labeled. The bifurcation of the second peak in Ow-Ow RDF for a c-length below 4 nm is also indicated. The positions of (C) the first peak in Ow-Ow RDF and (B) second peak in Hw-Hw RDF (the first peak is imposed by the equilibrium bond and angle in the respective harmonic potentials used in flexible SPC/E water).
* **Ettringite.** Calculation of *H-bond statistics* for ettringite shows that the contributions to the total H-bonds *HB* can be divided into H-bonds between water oxygen *HB*<sup>WW</sup>, columns OH groups *HB*<sup>CC</sup>, water oxygen and columns OH groups *HB*<sup>WC</sup>, sulfate oxygen and water *HB*<sup>SW</sup>, and sulfate oxygen and columns OH groups *HB*<sup>SC</sup> [Honorio et al., 2021a]. Due to the decrease in water content, *HB* as well as the H-bond involving molecular water (*HB*<sup>WW</sup>, *HB*<sup>SW</sup>, *HB*<sup>WC</sup>) decrease with desorption. The increase in *HB*<sup>CC</sup> and *HB*<sup>SC</sup> with desorption shows that direct interaction between OH groups from neighboring columns and between OH groups and sulfates becomes more significant with drying. These newly formed H-bonds contribute to the cohesion of the dried material, making it difficult for water molecules to ingress in the structure upon rewetting. In the ettringite domain, H-bonds breaking and formation are reversible and follow the same path as desorption, as indicated by the adsorption curves with rewetting at RH of 30 and 1%. The total H-bonds vary linearly with the water content irrespective of the sorption branch considered.

* **Crystalline ASR products.** The contributions to the total H-bonds *HB* can be divided in water-water *HB*<sup>WW</sup>, layer OH groups *HB*<sup>LL</sup> and water, and layer OH groups self *HB*<sup>WL</sup> contributions [Honorio et al., 2020f]. The number of *HB*<sup>WW</sup> and *HB*<sup>WL</sup> (and therefore *HB*) are ion specific following the leading order Cs > K > Na > Li. In all cases, the total H-bonds remain lower than the value of bulk SPC/E water (3.2 to 3.4 H-bonds per water molecule [Marcus, 2009]). The individual contributions of *HB*<sup>WW</sup> and *HB*<sup>WL</sup> follow the same trends. Properties of Na- and K-shlykovite can be linked to *HB* according to the overall trends: the elastic moduli (except for K-shlykovite) and heat capacity increases with *HB*, and the coefficient of thermal expansion decreases (except for α<sub>b</sub>) with *HB*. One possible explanation for these trends is that an increase in the *HB* per water molecule increases the constraints in the system, which translates into an increase in the rigidity and reduction of the coefficient of thermal expansion. This observation shows that one of the primary mechanisms through which charge balance ions change the thermo-mechanics of shlykovite is through the modification of H-bond networks.

## 3.4 Properties of confined and surface water

The physical properties of confined water can be computed using the following fluctuations formula in the grand canonical ensemble (μVT or hybrid μVT-NVT adapted to microporous materials as discussed in Section 3.1):

* Bulk modulus *K* (reciprocal of the isothermal compressibility β<sub>T</sub> = 1/*K*) [Allen and Tildesley, 1989]:

$$K = \frac{\langle N_w \rangle_{\mu VT}^2 k T}{\langle \delta N_w^2 \rangle_{\mu VT} V} \qquad (2.63)$$
* Specific heat capacity at constant volume $c_V$ [Shang, 2013, Schofield, 1966]:

$$c_V = \frac{1}{mkT^2} \left( \langle \delta E^2 \rangle_{\mu VT} - \frac{\langle \delta E \delta N_w \rangle_{\mu VT}^2}{\langle \delta N_w^2 \rangle_{\mu VT}} \right) \quad (2.64)$$

* Volumetric coefficient of thermal expansion $\alpha$ [Allen and Tildesley, 1989]:

$$\alpha = \frac{\langle P \rangle_{\mu VT}}{TK} - \frac{\langle \delta U \delta N_w \rangle_{\mu VT}}{\langle N_w \rangle_{\mu VT} k T^2} + \frac{\langle U \rangle_{\mu VT} \langle \delta N_w^2 \rangle_{\mu VT}}{\langle N \rangle_{\mu VT}^2 k T^2} \quad (2.65)$$

or generalizing for an anisotropic frame as [Honorio et al., 2021b]:

$$\alpha_{ij} = \frac{\langle P_{ij} \rangle_{\mu VT} \beta_T}{T} - \frac{\langle \delta U \delta N \rangle_{\mu VT}}{\langle N \rangle_{\mu VT} k T^2} + \frac{\langle U \rangle_{\mu VT} \langle \delta N^2 \rangle_{\mu VT}}{\langle N \rangle_{\mu VT}^2 k T^2} \quad (2.66)$$

* Molar incremental enthalpy $h_a$ [Nicholson, 1982]:

$$h_a = \frac{\partial H_p}{\partial N_w} \Big|_{T, P_w} = \frac{\langle \delta E \delta N_w \rangle_{VT}}{\langle \delta N_w^2 \rangle_{\mu VT}} \quad (2.67)$$

where $N_w$ is the number of water molecules; $m$ is the mass; $E$ and $U$ are the total and potential energies, respectively; $H_p$ is the enthalpy. The molar incremental enthalpy $h_a$ of the adsorbed phase describes how $H_p$ varies upon adsorption of $N_w$ molecules of fluid at a fixed $T$ and external fluid pressure $P_w$. This quantity is directly associated with the isosteric heat of adsorption $q_{st} = h_{bulk} - h_a$, with $h_{bulk}$ being the molar enthalpy of the bulk fluid [Abdolhosseini Qomi et al., 2021].

Figure 2.24 shows the properties of **confined water** in C-S-H under liquid saturated conditions at 300 K [Honorio et al., 2021b]. The fluctuation formula above was used to get the effective property of confined water but without accounting for the heterogeneity of properties along the pore thickness. The thermomechanical properties of water confined in C-S-H depend on pore size, following the example of transport and electromagnetic properties [Manzano et al., 2012, Honorio, 2019, Masoumi et al., 2019, Youssef et al., 2011]. For the bulk modulus, heat capacities, and coefficients of thermal expansion, the bulk values are retrieved for interlayer distances as small as 16.3-26.3 Å (i.e., c-lengths of approximately 3-4 nm). In smaller pores, the properties of water are more solid-like, in agreement with observations of transitions liquid-to-glassy and glassy-to-solid, also reported for self-diffusion (as discussed in Section 2.4.1).

To assess properties of **surface water**, the response at the water layer scale is assumed to be close to grand canonical conditions: the water molecules in each layer can exchange with the water molecules in the neighboring layers [Honorio et al., 2022d]. This assumption is valid if the water layers thickness is much smaller than the total pore thickness so that this pore be considered as a large reservoir with
<figure><img src="image_27.png" /><figcaption>**Figure 2.24:** Apparent density $\rho^*$; bulk modulus $K$, heat capacities $C_V$ and $C_P$; coefficient of thermal expansion $\alpha_{ii}$ and $\alpha_v = \text{Tr} \, \alpha_{ii}/3$; isosteric heat of adsorption $q_{st}$ as a function of the interlayer distance in C-S-H [Honorio et al., 2021b] compared to bulk SPC/E water and experimental data on bulk water.</figcaption></figure>

respect to each layer, approximating, therefore $\mu$VT conditions. Thus, a simulation box treated in the NVT ensemble allows computing the local behavior (per layer) by approximating the behavior of each layer as in $\mu$VT ensemble. Figure 2.25(A) shows a sketch of the local-$\mu$VT conditions; in (B) the exchange of the molecules among
<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

**Figure 2.25:** (A) Sketch of the local-μVT conditions used to compute the properties of surface water using fluctuation formula in μVT ensemble [Honorio et al., 2022d]. (B) Initial layering (at time $t=0$) in C-S-H following a color code per layer and how the system evolves to a mixed state in a timescale of 100 ps. Layers were defined according to the density profiles of water oxygens Ow.

the layers is highlighted by showing the initial layering (at time $t=0$) following a color code and how the system evolves to a mixed state (especially to the layers towards the bulk zone) in a timescale of 100 ps.

Adopting the local-μVT approximation per layer, the same fluctuation formula can be deployed to compute properties in the μVT ensemble. The layers were defined according to the density profiles of water oxygen: each oscillation in these
<figure><img src="image_30.png" /><figcaption>**Figure 2.26:** Properties per water layer according to two drained poromechanical conditions [Honorio et al., 2022d]: potential energy *U*, isochoric specific heat capacity *c*<sub>v</sub>, molar incremental enthalpy *h*<sub>a</sub>, bulk modulus *K*, and (volumetric) coefficient of thermal expansion *α*.</figcaption></figure>

profiles determines a layer. Figure 2.26 shows how the potential energy *U*, isochoric specific heat capacity *c*<sub>v</sub>, molar incremental enthalpy *h*<sub>a</sub>, bulk modulus *K*, and (volumetric) coefficient of thermal expansion *α* averaged per water layer evolve for both poromechanical conditions considered. All these parameters show a marked dependence on the distance from the C-S-H adsorbing surface. The first three layers are directly influenced by surface roughness and the presence of microportlandite clusters. For all the properties considered, water retrieves the bulk values after layers 4 or 5. Property data per layer can be useful in revisiting the sorption models generally used for cement-based materials, especially with respect to the assumption made regarding mono and multi-layer structuring, adsorption energy, and boundary between bulk and surface water.
# 4 Conclusion

This chapter showed how molecular simulations are a toolbox to compute missing or hard to (experimentally) determine properties of phases present in cement systems. I presented how these simulations have been used in my works to provide property data on various phases in cement systems: C-(A)-S-H [Honorio, 2019, Honorio et al., 2021b, Honorio et al., 2022b, Hawchar and Honorio, 2022], AF-phases [Honorio et al., 2020e, Honorio et al., 2021a, Honorio, 2022b, Honorio et al., 2022a], ASR products [Honorio et al., 2020d, Honorio et al., 2020f, Honorio, 2021], bulk pore solutions [Honorio et al., 2019a, Honorio et al., 2020a], and clinker minerals [Honorio et al., 2020e, Barbosa and Honorio, 2022]. And that in a multi-physics framework, including elastic properties, thermal properties, mass transport properties, electromagnetic properties, and coupling properties. For most properties, the anisotropy in the response at the unit cell level was quantified. In some cases, the results obtained from simulations contributed to the mission of completing property data databases for the various phases that can be found in cement systems [Bullard et al., 2019].

The second part of this chapter coped with the molecular mechanisms of sorption. The properties of confined and surface water were quantified using the adapted techniques from statistical physics and molecular simulations. The quantification of effective interaction provides a critical input for simulations dealing with the mesoscale associated with the atomistic-to-continuum bridge. A strategy to quantify microporous phases sorption-induced deformations was discussed in applications for C-S-H [Masara et al., 2023, Honorio et al., 2022d], AF-phases [Honorio et al., 2021a] and ASR products [Honorio et al., 2020d]. Physical insights on sorption processes regarding water layering, the role of (heterogeneous) cavitation in drying, and the origins of hysteresis under sorption will allow future work to revisit the sorption models used in cement-based materials research.
# Chapter 3

## Upscaling properties: from nano to macro

### Contents

<table><tr><td>1</td><td>Introduction</td><td>62</td></tr><tr><td>2</td><td>Mesoscale(s)</td><td>64</td></tr><tr><td>3</td><td>Coarse-grained simulations</td><td>64</td></tr><tr><td>3.1</td><td>ASR gels</td><td>64</td></tr><tr><td>3.2</td><td>C-S-H gels</td><td>68</td></tr><tr><td>4</td><td>Micromechanics</td><td>75</td></tr><tr><td>4.1</td><td>Representation of the microstructure</td><td>75</td></tr><tr><td>4.2</td><td>Conductivity-like properties</td><td>78</td></tr><tr><td>4.3</td><td>Elastic properties</td><td>86</td></tr><tr><td>4.4</td><td>Aging linear viscoelasticity</td><td>89</td></tr><tr><td>4.5</td><td>Thermal expansion</td><td>92</td></tr><tr><td>4.6</td><td>Heat capacity</td><td>93</td></tr><tr><td>4.7</td><td>Local information in elasticity and strength estimates</td><td>95</td></tr><tr><td>5</td><td>Conclusion</td><td>100</td></tr></table>
# 1 Introduction

In multiscale multi-physics modeling, upscaling techniques are required to link composition to effective properties at the macroscale or the scale of applications. Micromechanics and Statistical physics are powerful tools in this framework.

Micromechanics is the study of the mechanical, thermal, electromagnetic, and transport behavior of materials with a microstructure (or heterogeneous materials) [Milton, 2002, Torquato, 2002]. This theory provides estimations of macroscopic properties of heterogeneous materials in a multi-physics framework using as input the system's composition, assumptions on the microstructure, and component property (i.e., the properties of each phase present in the material). With micromechanical models, it is possible to associate the relevant physical phenomena to the appropriate scale and then model the macroscopic behavior of the materials, properly accounting for their multiscale nature. The theory can be extended to upscale complex phenomena, including aging viscoelasticity [Sanahuja, 2013, Honorio et al., 2016b, Honorio, 2017], mechanical strength [Pichler et al., 2009, Pichler and Hellmich, 2011, Pichler et al., 2013, Honorio et al., 2018b], as well as complex frequency-dependent dielectric permittivity spectra [Hashin and Shtrikman, 1963, Guihard et al., 2020, Honorio et al., 2020a] in cement-based materials.

However, the atomistic-to-continuum bridged for some key phases in cement systems is not trivial due to the existence of a critical intermediary level that must be treated independently. This is the case of C-(A)-S-H and ASR products (and also clays), for which a "granular"/colloidal (meso)scale is present. This nano-to-micro mesoscale is too big for full atomistic simulation but too small to direct apply (without validation) some results from usual continuum mechanics. Full atomistic molecular dynamics simulations are generally limited to simulation box sizes on the order of tens to hundreds of nanometers and simulated times on the order of the microsecond. Representative Elementary Volumes (REV) of the nano-to-micro mesoscale in the case of layered silicates (e.g., C-S-H, clays) are on the order of at least a few hundreds of nanometers. Nanoscale phenomena such as thermal fluctuations and couplings between chemical transformations and mechanical deformation are significant at this mesoscale [Ioannidou et al., 2022]. In this context, coarse-grained (CG) simulations arise as a powerful tool to overcome time and length scale limitations in full atomistic simulations and provide information on the nano structuration of phases with colloidal particles, and that can be made using input from the molecular scale. These simulations follow the principles of Statistical Mechanics, whose main goal "is to relate the macroscopic properties of a system of many particles (atoms, molecules, spins, etc.) to its microscopic properties, which include the interparticle interactions as well as the spatial statistics of the particles" [Torquato, 2002]. This theoretical framework enables describing the mesoscale systems rigorously [Ioannidou et al., 2022], justifying the choice for coarse-graining the interactions and the matter granularity, as well as the upscaling of information from the molecular scale. For calculation of properties such as elasticity, the strategy
<table><caption>**Figure 3.1:** Upscaling strategies in a multiscale multi-technique framework.</caption><thead><tr><th></th><th>Molecular<br>Scale</th><th>Colloidal<br>Mesoscale</th><th>Continuum Micro and<br>Mesoscales</th></tr></thead><tbody><tr><th scope="row">Strategy 1</th><td></td><td></td><td>Micromechanics</td></tr><tr><th scope="row">Strategy 2</th><td>Molecular Simulations</td><td></td><td>Micromechanics</td></tr><tr><th scope="row">Strategy 3</th><td>Molecular Simulations</td><td>Coarse-Grained Sim.</td><td>Micromechanics</td></tr><tr><th scope="row">Strategy 4</th><td></td><td>Coarse-Grained Sim.</td><td>Micromechanics</td></tr></tbody></table>

presented in the last chapter remains valid for coarse-grained simulations.

The upscaling strategies in a multiscale multi-technique framework are schema- tized in Figure 3.1. Strategy 1 corresponds to using Micromechanics and available data on constituent properties. The available data here can be properties measured at the micrometer scale (e.g., by indentation) or obtained from "pure" phases tested macroscopically. Strategy 2 corresponds to the use of molecular simulations (as de- scribed in Chapter 2) to get property data and then use this data directly as input in Micromechanics. Strategy 3 corresponds to the use of molecular simulations to get effective interactions to be used as input in coarse-grained simulations, which finally gives property data to be used as input in Micromechanics. Finally, strategy 4 is to perform coarse-grained simulations using general interaction potentials or interaction tuned to capture some behavior; the results can be then used as input in Micromechanics to get the properties at the macroscale. Strategies 1 to 3 have been deployed in my works, as presented next.

The chapter discusses the mesoscale(s) of interest in the passage atomistic-to- continuum and also in the realm of continuum mechanics. It is shown how coarse- grained simulations give information on the nano structuration of phases with col- loidal particles. Then, it is discussed how Micromechanics, a continuum-medium- based theory, and the Statistical Physics language of coarse-grained simulation pro- vide consistent information about the "granular" mesoscale in colloidal phases such as C-S-H. It is presented how Micromechanics theory provides estimates of effective properties in a multi-physics context and also information about local fields. Up- scaling strategies of aging viscoelasticity related properties, compressive strength, transport properties, thermal properties, and electromagnetic properties are pre- sented in applications to cement-based materials.
## 2 Mesoscale(s)

Various intermediary scales - or mesoscales - between atoms and structural applications are relevant for cement and concrete research. According to the applications, different mesoscales are targeted. In concrete modeling, "mesoscale" generally refers to the scale associated with the REV, in which the heterogeneities have a characteristic size on the order of tens of µm up to a few mm. The heterogeneities at the mesoscale comprise aggregates (coarse and/or fine), a cement paste or mortar matrix, an interfacial transition zone (ITZ), and in some cases, a macroporosity (e.g., air voids). This definition is to be confronted with the usual scale definitions based on Micromechanics textbooks (e.g., [Milton, 2002]) in which the microscale refers to the scale of the heterogeneity itself, the mesoscale to the scale of the REV and the macroscale to the scale of the structure for which the behavior can be assumed to be homogeneous. The mesoscale of the cement hydration community is that of microstructure development (from tenths up to thousands of meters) [Ioannidou et al., 2022]. These two mesoscales are relevant for micromechanics modeling.

Finally, the nano-to-micro "colloidal" mesoscale, already mentioned above, is of particular interest because of its significance for the comprehension of gel structure and properties. Techniques allowing assessing this mesoscale are discussed in the sequel.

## 3 Coarse-grained simulations

Coarse-grained simulations on ASR and C-S-H gels are presented following recent work [Ait-Hamadouche and Honorio, 2022b, Ait-Hamadouche and Honorio, 2022a].

### 3.1 ASR gels

To date, most of the studies on ASR focus on experimental characterization in the macroscopic scale, or in a scale in which ASR products appear mixed with other phases in cement-based materials. More recent developments provide atomic-level understanding of the structure of crystalline [Shi et al., 2019, Honorio et al., 2020d] and disordered [Dupuis and Pellenq, 2022] ASR products. But concerning the intermediary gel scale, there is a knowledge gap regarding ASR products' structural features and properties.

Coarse-grained simulations can help by providing estimates of gel properties that remain elusive. The viscous behavior is acknowledged to be critical in understanding ASR development [Struble and Diamond, 1981, Ichikawa and Miura, 2007, Gholizadeh Vayghan et al., 2016, Bažant and Steffens, 2000, Bažant and Rahimi-Aghdam, 2017]. However, few works focus on measuring the viscous properties of ASR gels. The picture of ASR gels as viscous fluids that, upon formation, immediately flow into the neighboring pores is often evoked to explain the distribu-
tion of ASR products in cement-based materials [Gholizadeh Vayghan et al., 2016]. The viscosity of ASR gels and how it evolves with gelation is yet to be quantified.

<figure><img src="image_31.png" /><figcaption>**Figure 3.2:** Snapshots of the mesostructures of Na-ASR and K-ASR gels as a function of the packing density $\eta_g$: networks (at top) and isosurfaces (at the bottom).</figcaption></figure>

In a recent work [Ait-Hamadouche and Honorio, 2022b, Ait-Hamadouche and Honorio, 2022a], the capability of CG simulations to capture gel mesotexture formation and properties of ASR products is tested. The PMFs obtained from full atomistic simulations for Na- and K-ASR products reported in [Honorio et al., 2020d]
are used as input in the CG simulations. Using PMFs obtained from atomistic sim-
ulations has the advantage of better describing the details relevant to the system
studied instead of relying on general function forms. The functional forms of the
PMFs in Figure 2.20 deviates from the potentials often deployed in CG simulation of
calcium silicates (i.e., Mie or generalized Lennard-Jones potentials with or without
Yukawa term [Masoero et al., 2012, Ioannidou et al., 2014, Ioannidou et al., 2016])
by showing a non-monotonous increase in the energy in the attractive portion. This
locally oscillatory behavior leads to local changes in the forces correlations.

<figure><img src="image_32.png" /><figcaption>**Figure 3.3:** Eigenstresses evolution with densification of Na-and K-ASR gels for various packing densities [Ait-Hamadouche and Honorio, 2022b]. The PDFs are depicted in a log-scale for better identification of the behavior of the distribution tails.</figcaption></figure>

In the coarse-grained simulations, grains are modeled as spherical particles with
diameter 20 Å consistent with experiments [Gaboriaud et al., 1999]. The use of
monodisperse spheres limits the maximum packing density η_g that can be obtained
under zero stress (e.g., [Torquato, 2002]): η_g = 0.64 for the maximum random
packing density and η_g = 0.74 for the close-packing for face-centered cubic (fcc) or
hexagonal close-packed (hcp) spheres. Monodisperse spheres are sufficient to fully
describe the range of packing density reported in the literature for ASR gels: these
gels would present a porosity f_g (= 1 − η_g) ranging from 0.4 to 0.8 [Geng et al.,
2021] (i.e., a packing density ranging from 0.2 to 0.6).
<figure><img src="image_33.png" /><figcaption>**Figure 3.4:** Evolution of Young *E* and shear *G* moduli of Na-ASR (top) and K-ASR (bottom) as a function of the packing density at zero stress η₀ [Ait-Hamadouche and Honorio, 2022b]. The self-consistent estimate, using the Voigt-Reuss-Hill (VRH) bulk and shear moduli obtained from molecular simulations on Na and K-shlykovite as input, is provided for comparison.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>**Figure 3.5:** Shear (left) and bulk (right) viscosities as a function of the packing density η_g for Na- and K-ASR gels [Ait-Hamadouche and Honorio, 2022b]. The experimental data of [Yang et al., 2008] on static shear viscosity of sodium silicate gels (the gray curves) are shown for comparison.</figcaption></figure>
The gelation and densification processes are compared according to various struc-
tural features, including the pore size distribution (PSD), specific (iso-)surface area,
radial distribution functions, and scattering intensity. Figure 3.2 shows snapshots of
the mesostructures of Na-ASR and K-ASR gels as a function of the packing density.

Coarse-grained simulation provides a description of ASR gel structural features
and properties in agreement with the available experimental data [Ait-Hamadouche
and Honorio, 2022b]. The simulations yield numerical samples of ASR gels that are
structurally consistent with SAXS results [Gaboriaud et al., 1999]. The elastic prop-
erties and viscosity computed from simulations reasonably agree with the available
experimental data, which builds confidence in using coarse-grained simulations to
unveil details of the structure and behavior of ASR gels that remain elusive or hard
to probe experimentally. As demonstrated in the case of C-S-H, CG simulations
should be seen as a valuable tool for investigating ASR gels at the mesoscale.

The ASR gel mesotexture is shown to be alkali-specific. The analysis of gel
structure shows a clear dependence on the effective interactions (or PMF), which
are a function of the charge balance alkali in ASR products. The Na-ASR PMF is
multistable while the K-ASR only shows one clear minimum; also, Na-ASR PMF
goes slightly slower to zero than K-ASR PMF. Microstructural descriptors affected
by the charge balance alkali are the pore size distribution (PSD), specific surface
area, particle neighbor arrangements, and cluster size. Na-ASR shows a coarser
mesotexture with larger clusters than in K-ASR gel. Both PSD are bi-modal in
denser systems with a pore class associated with micropores akin to interlayer pores
and another pore class associated with mesopores (or gel pores), as also observed
in C-S-H [Ioannidou et al., 2022]. K-ASR shows a larger specific surface area than
Na-ASR. These results and the slightly larger gel pores in Na-ASR suggest that the
ion specificity may play a role in the transport properties and swelling of the gel.

The nanomechanics of ASR gel morphology is also ion-specific. The eigenstresses
build-up (Fig. 3.3), elastic properties (Fig. 3.4), and viscosity (Fig. 3.5) also exhibit
a dependence on the charge balance alkali in ASR products. The Na-ASR gel
shows more considerable variability in the eigenstresses developed upon densifica-
tion. Na-ASR products show larger elastic Young and shear moduli than K-ASR,
and consistently also a larger (bulk and shear) viscosity than K-ASR. Both viscosi-
ties increase exponentially with the packing density, in particular after percolation.
Both gels exhibit a viscosity akin to honey or ketchup at the packing densities just
after percolation.

These results might be helpful in future work focusing on the viscoelastic be-
havior of ASR gels. Further works may cope with morphologies more complex than
monodisperse spheres and other ions relevant to ASR products aging.

## 3.2 C-S-H gels

For C-S-H gels, a multi-scale strategy is also adopted but following an approach in
which layers are explicitly represented and their flexibility is accounted for. First, the
general theory of nanolayers flexibility is presented, followed by the coarse-grained simulation of C-S-H *per se*.

### 3.2.1 Theory of flexibility of nanolayers

The general form of the free energy associated with the flexibility of a nanolayer is given by Landau-Helfrich bending free energy [Helfrich, 1973]:

$$U_{LH} = \int_S \left( \frac{1}{2} D_i^f (\bar{\kappa} - \kappa_0)^2 + k_g K_g + \gamma \right) dS \quad (3.1)$$

where $D_i^f$ is the bending modulus; $\bar{\kappa} = (\kappa_1 + \kappa_2)$ is the mean curvature in terms of principal curvatures $\kappa_1$ and $\kappa_2$; $\kappa_0$ is the spontaneous curvature; $K_g = \kappa_1\kappa_2$ is the Gaussian curvature; $k_g$ is the saddle-splay modulus; $\gamma$ is an external tension; and the integral is performed over the layer surface $S$. The curvature in a direction $a$ contained in the xy plane is related to the displacement $\delta(x, y)$ orthogonal to xy plane by $\kappa_a = \frac{\partial^2\delta(x,y)}{\partial a^2}$. This bending free energy remains valid provided that the in-plane lengths $L$ of the membrane are much larger than the atomic granularity (related to the typical size of an atom, i.e., a few Å) and much smaller than the persistence length. The persistence length $l_p$ is a fundamental characteristic of flexibility, frequently used in colloid and polymer science. Using de Gennes and Taupin definition [De Gennes and Taupin, 1982]:

$$l_p = a \exp \left[ 2\pi D_i^f / (k_B T) \right] \quad (3.2)$$

where $a$ is the monomer size. The persistence length is a characteristic size defining when a chain or membrane can be treated as a rigid particle (if size $L \ll l_p$) or as a crumpled random landscape due to thermal fluctuations (if $L \gg l_p$). In other words, $l_p$ can be viewed as the minimal length above which the chain or 2D-particle crumples under thermal agitation [Suter et al., 2007]. From Eq. 3.2, $l_p$ increases exponentially with the bending stiffness. Theoretical considerations show that in the case of foil-like particles, when their average distance $d < l_p$, the particles are expected to stack with increasing packing density [De Gennes and Taupin, 1982].

In the case of C-S-H, no spontaneous curvature is expected ($\kappa_0 = 0$) since the layers are overall symmetric with respect to their mid-plane (This is not the case of other phyllosilicates such as chrysotile, and potentially M-S-H [Bernard, 2022]). It can be reasonably assumed that no topology changes occur due to thermal fluctuations, so the term on $k_g$ remains constant [Helfrich, 1973]. Finally, we may adopt zero external surface tension. Thus, only the term related to the bending modulus in Eq. 3.1 is relevant for C-S-H. C-S-H flexibility is expected to be on the order of $10^{-17}$ J as smectite [Suter et al., 2007, Honorio et al., 2018d] and tobermorite [Honorio and Brochard, 2017] layers. Thus, $2\pi D_i^f / (k_B T) \approx 10^4$, meaning that the persistence length in C-S-H should be several orders of magnitude larger than an Angstrom. For comparison, the persistence length of smectite layers is on the order of 100 nm [Suter et al., 2007]. Furthermore, considering that the C-S-H layer length
$L$ is much larger than its thickness $h_s$ and the atom size, thin plate (Kirchhoff-Love) theory can be used (Figure 3.6(A)), and the free energy associated with bending takes the form:

$$U_{LF} = \frac{1}{2} \sum_{i,j,k,l \in \{x,y\}} D_i^f \frac{\partial^2 u_z}{\partial i \partial j} \frac{\partial^2 u_z}{\partial k \partial l} \quad (3.3)$$

where $u_z$ is the displacement orthogonal to the plane of the plate. For a thin plate of thickness $h$, the bending modulus in direction $i$ is given by:

$$D_i^f = \frac{h^3 E_i}{12(1 - \nu_{ij}\nu_{ji})} \quad (3.4)$$

with $h$ being the layer thickness, $E_i$ the Young's modulus in the $i$ direction and $\nu_{ij}$ the Poisson's ratio in the $i$ direction under a load in the perpendicular in-plane direction $j$.

The texture of nanolayered adsorbing minerals also depends on the interplay between flexibility and effective interactions. Figure 3.6(B) depicts two extreme types of micro textures of nanolayered materials: stacked or "brick-like" and entangled or "spaghetti-like" types [Honorio et al., 2018d, Zabat et al., 1997]. Rigid layers with strong interparticle cohesion forces would favor a "brick-like" texture; flexible layers with less pronounced interparticle cohesion would favor a "spaghetti-like" texture. The implication to transport and mechanical properties are significant [Zabat et al., 1997]. Brick-like texture exhibits larger pore sizes and fewer contact points between stacks, leading to larger permeability (which is linear with the porosity but generally scales with the square of pore size) and more brittle failure. Conversely, spaghetti-like texture exhibits smaller pores and a much larger contact surface between particles, which lead to lower permeability and more ductile failure. The existence of bent configurations in the mesoscale is also to be assessed as a function of the particle aspect ratio and concentration regime. Long thin particles are more prone to show bent configurations [Honorio et al., 2018d]. As discussed above, in concentrated regimes (not very flexible), nanolayers are expected to stack coherently [De Gennes and Taupin, 1982].

The characteristic length $l_{trans}$ (Figure 3.6(C)) associated with the interplay between the flexibility and the effective interactions reads [Honorio et al., 2018d]:

$$l_{trans} = \left( \frac{(\Delta d)^2 D_i^f}{\Delta \Lambda_P / S} \right)^{1/4} \quad (3.5)$$

where $\delta d$ is the variation in the interlayer distance, $\Delta \Lambda_P / S$ is the associated variation in the effective interactions per surface area $S$. For $l_{trans} \gg L$, the nanolayer is comparatively too rigid, and very little bending is observed. Conversely, for $l_{trans} \ll L$, the nanolayer is comparatively very flexible, deforming significantly within a short distance with respect to the layer dimension when interspersed inbetween layers not sharing translational or orientational correlations. Considering
A

**Single layer scale**

<figure><img src="image_35.png" /></figure>

**Particle scale**

<figure><img src="image_36.png" /></figure>

B

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

C

<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /><figcaption>**Figure 3.6:** (A) Separation of scales in nanolayers and stack of layers scales: conditions under which the assumptions of homogeneous plate theory holds [Honorio et al., 2018d]. (B) Representation of the "extreme" nanotextures of nanolayers materials (e.g., clays, C-S-H): (left) stacked or "brick-like" and (right) entangled or "spaghetti-like" types [Honorio et al., 2018d, Zabat et al., 1997]. (C) At left, sketch of a simple stack of clay layers undergoing a transition of hydration state from (n + 1)W to nW state. At right, a sketch of the transition between two hydration states and the superposition of solid layers: *l*<sub>trans</sub> characterizes the distance over which a layer bends in-between stack due to the interplay of flexibility and hydration forces [Honorio et al., 2018d].</figcaption></figure>
the energy barriers associated with the effective interaction in C-S-H, the length over which the layers bend during a hydration transition or mechanical load is expected to be in the order of a few nanometers.

### 3.2.2 Coarse-grained simulations

Various morphologies are generally identified for C-S-H including [Richardson, 2004, Zhang et al., 2018]: isotropic globules, foil-like (2D) particles, and fibrillar (1D) particles. Strong experimental evidence supports the 2D foil-like structures as one of the most relevant to C-S-H gel [Richardson, 2000, Richardson, 2004, Richardson et al., 2010, Scrivener et al., 2019]. Herein, 2D foil-like morphology is considered here since it is the most prevalent in C-S-H. Simulations at the molecular scale are used to get potentials of mean force (PMFs) describing:

* Intraparticle interactions: the uniaxial elasticity in the layer plane and layer's flexibility (2- and 3-body interactions, respectively).

* Interparticle interactions: layer face-to-face effective interactions under liquid saturated drained conditions based on grand potential interactions obtained in [Honorio et al., 2021b] for a realistic C-S-H model [Kunhi Mohamed et al., 2018] with Ca/Si of 1.67.

These PMFs are used as input in coarse-grained simulations at the gel scale. The configurations obtained from these simulations are then used to assess the structural features and properties at the gel scale. We compare the simulation results with experimental data from the literature whenever possible.

To compute the intraparticle PMFs, we calculate the rigidities $k_{2b}$ and $k_{3b}$. We consider the C-S-H layer as a structure discretized in spherical mass points, with 2-body interactions defining the stretch potential between mass points i and j ($U_{2b}$). Similarly, 3-body interactions define bending potential between mass points i, j, and m ($U_{3b}$). Then, this approach is combined with linear elastic beam theory by equating the harmonic potentials used in CG simulations with the free energy expressions of beam theory [Keremides et al., 2018]. Therefore $k_{2b}$ and $k_{3b}$ take the forms:

$$k_{2b} = EA/L_a; \quad k_{3b} = 12EI_a/L_a \qquad (3.6)$$

where E is the (in-plan) Young's modulus, $L_a = 12.85$ Å (=σg) is the distance between the center of two mass points (taken as the diameter of each spherical mass), $A = πσ_g^2/4$ is section of each mass point, and $I_i = \frac{1}{4}π(σ_g/2)^4$ is the moment of geometric inertia. The bending modulus $D_i^f$ and rigidities obtained with this approach are displayed in Table 3.1.

These bending modulus values are in the same order as that found in a previous study on the flexibility of tobermorite [Honorio and Brochard, 2017]. Simulations on a single fibril (chain of spherical grains) are performed to verify this approach. In these simulations, a flexion effort was generated through an imposed displacement
<table><caption>**Table 3.1:** Bending modulus and rigidities ($k_{2b}$ and $k_{3b}$) of C-S-H layers.</caption><thead><tr><th>D<sub>i</sub><sup>f</sup> [N.m]</th><th>1.66×10<sup>-17</sup></th></tr></thead><tbody><tr><td>k<sub>2b</sub> [kcal/mol.Å<sup>2</sup>]</td><td>80.3</td></tr><tr><td>k<sub>3b</sub> [kcal/mol.rad<sup>2</sup>]</td><td>9950</td></tr></tbody></table>

$\delta_i = \delta_{i-1} + \Delta\delta$, with $\Delta\delta = 0.02$ Å at the distance of one-quarter from each edge of the chain. Based on the Kirchhoff-Love plate theory, for this configuration, the elastic free energy is a quadratic function of the displacement $\delta : U = 6(\frac{4}{L})^3 D_i\delta^2$. Our values agree with Kirchhoff-Love plate theory, and also the expression proposed by [Prathyusha et al., 2018] linking $D_i^f$ and $k_{3b} : D_i^f = k_{3b} \times 2\frac{\sigma}{L_{tot}}$ (with $L_{tot}$ being the length of the spherical grain chain).

To investigate the effect of layers flexibility on the meso-structuration of C-S-H gel, we launched three simulations named F01, F1, and F10, meaning that the rigidities used as input in coarse-grained simulations were multiplied respectively by 0.1, 1, and 10. The grain size was chosen based on experimental values and in agreement with the PMF. The grain's diameter with a value of 12.85 Å corresponds to the arg[$U(r) = 0$] value in the interparticle PMF. The disk thickness is 308.4 Å, corresponding to the average diameter reported for a disk shape in C-S-H particle [Maruyama et al., 2021].

In a simulation box of 1000 nm, 400 disks are randomly placed. Then, we adopt the method proposed by [Masoumi et al., 2020] using high-frequency pressure oscillations and annealing in MD simulations to bring the system to the most favorable packing density. We use the Nosé-Hoover thermostat and barostat. We begin our simulation with small timesteps to avoid numerical instabilities, then we perfom a significant number of cycles, with a cycle being composed of two stages: (i) NVT runs at a temperature $T = 3000$ K and then at ambient conditions ($T = 300$ K); and (ii) NPT runs at $T = 300$ K with at first a pressure $P = 200$ atm, followed by a gradual decrease to 1 atm.

The aggregation of C-S-H layers forming the C-S-H gel is displayed in Fig. 3.7(a). The three configurations, F01, F1, and F10, composed of stacks of layers, have differences that can be visually identified. For example, the layers in F01, being less stiff, appears in more pronounced bent/curved configurations than those in F1 and F10. Moreover, for F1, there is predominately face-to-face contact between the layers, while for F10, there is also face-to-edge contact. Some experimental and simulation studies have investigated the number of C-S-H layers in stacks: experiments found an average number of layers $n = 10.9$ [Chiang et al., 2012], simulations found clusters containing up to forty stacked layers at $\eta_g = 70\%$ [Masoumi et al., 2020]. In our simulations, it can be visually identified stacks constituted of at least 10 layers. Finally, Fig. 3.7 (b) gives a better view of the space occupied by the solid part and the remaining void.

We also compute the effective Young modulus of C-S-H as a function of the packing density at zero stress $\eta_0$ according to the flexibility of the layers (Fig. 3.8).
<figure><img src="image_41.png" /><figcaption>**Figure 3.7:** Layers aggregation (a) and isosurface (b) of C-S-H gel according to the flexibility of the layers for $\eta_g = 0.25$ [Ait-Hamadouche and Honorio, 2022a].</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>**Figure 3.8:** Effective Young modulus of C-S-H a function of the packing density at zero stress $\eta_0$ according to the flexibility of the layers. The estimate of a self-consistent scheme for ellipsoidal inclusion with an aspect ratio of $a_r = 12.85/154.2 = 0.0833$ is provided for comparison. [Ait-Hamadouche and Honorio, 2022a]</figcaption></figure>

Higher flexibility leads to a slightly larger Young modulus, but the results are within the standard deviations and are therefore inconclusive. For comparison, we also provide the estimates of the self-consistent scheme for ellipsoidal inclusion with an aspect ratio of $a_r = 12.85/154.2 = 0.0833$. The Voigt-Reuss-Hill (VRH) moduli of C-S-H calculated at the molecular scale are used as input in the homogenization. The percolation threshold of the self-consistent scheme with this oblate inclusion is $\eta_p = 0.28$, which is in agreement with previous work [Sanahuja et al., 2007]. CG simulation leads to a slightly earlier percolation. Other packing densities are yet to
be assessed. For this, longer simulations (on the order of the µs [Masoumi et al., 2020]) are needed to get configurations at larger packing densities.

Perspectives include using polydisperse disks to represent the C-S-H gel to reach higher packing fractions associated with low and high-density C-S-H. Also, future work can deal with the computation of the other properties of gels.

# 4 Micromechanics

Strategies deployed in my works to get effective properties using Micromechanics in a multi-physics context are briefly described. The application of these strategies in the study of cement systems made are presented.

## 4.1 Representation of the microstructure

In the choice of the representation of the microstructure, the scale separation criterion must be considered (e.g., [Zaoui, 2002]):

$$d_0 \ll d_i \ll l_{REV} \ll l_{str} \qquad (3.7)$$

where $d_0$ is the size under which continuum mechanics is no more valid, $d_i$ is the characteristic length of the inhomogeneities, $l_{REV}$ is the size of the REV, and $l_{str}$ is the characteristic size of the structure in which the behavior of the material is to be considered as homogeneous.

<figure><img src="image_43.png" /><figcaption>**Figure 3.9:** Representation of the microstructure of concrete adopted in [Honorio et al., 2016b] to upscale aging linear viscoelastic properties: an approach combining analytical and numerical homogenization.</figcaption></figure>
<figure><img src="image_44.png" /><figcaption>**Figure 3.10:** Cement paste microstructure considered for micromechanics upscaling schemes in [Honorio et al., 2022c]: Input volume fractions are obtained from Powers, Königsberger-Hellmich-Pichler (KHP) or Tennis and Jennings (TJ) models. SC or MT are deployed to upscale cement paste elastic properties. For 2-scale representations, C-S-H gel effective properties are upscaled using the SC scheme. All the other phases, including pores, are considered as spherical inclusions. SC scheme is deployed to upscale cement paste elastic properties. For the representations with ellipsoidal inclusions: with the Powers model, hydrates are considered as prolate particles $a_r = 10$; with the TJ model, C-S-H LD and HD are considered as prolate particles $a_r = 10$, and CH and AFm as oblate particles $a_r = 0.2$.</figcaption></figure>

Various representations of the microstructure of cement-based materials exist with varying number of scales, assumptions on system morphology and inclusion shape, and the use of different models to describe the volume fraction of the constituents. Figures 3.9 and 3.10 show the microstructures representation adopted in