# Modelling the poromechanical behaviour of class G cement paste: A multiphysics approach from early age to hardened state

Maxime Pierre, Marcos Samudio, Siavash Ghabezloo, Patrick Dangla

► To cite this version:

Maxime Pierre, Marcos Samudio, Siavash Ghabezloo, Patrick Dangla. Modelling the poromechanical behaviour of class G cement paste: A multiphysics approach from early age to hardened state. Cement and Concrete Research, 2025, 193, pp.107852. 10.1016/j.cemconres.2025.107852. hal-05079783

HAL Id: hal-05079783

https://enpc.hal.science/hal-05079783v1

Submitted on 22 May 2025

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

# Modelling the poromechanical behaviour of class G cement paste: A multiphysics approach from early age to hardened state

Maxime Pierre<sup>D</sup>, Marcos Samudio, Siavash Ghabezloo<sup>D</sup>,*, Patrick Dangla<sup>D</sup>

Navier, ENPC, Institut Polytechnique de Paris, Univ Gustave Eiffel, CNRS, Marne-la-Vallée, France

## ARTICLE INFO

### Keywords:

Oil-well cement

Poromechanics

Multiphysics couplings

Cement hydration

Constitutive modelling

## ABSTRACT

Modelling cement-based materials from the early-age to the hardened state is crucial in numerous applications such as deep well cementing or 3D printing, which require comprehensive modelling of multiphysics couplings. To answer these requirements, a thermodynamically consistent time-dependent constitutive model based on the extent of hydration is developed in the framework of thermoporomechanics of partially saturated materials. Using minimal fitting, complex undrained oedometric tests on hydrating cement paste, combining effects of hydration progress, pore pressure evolution, elastic, viscous, and plastic deformations, are well reproduced numerically. In particular, the impact of early-age loading on the behaviour at a subsequent age, paramount in oil-well applications to understanding the consequences of pressurising the casing when the cement sheath is partially hydrated, is explained and quantitatively reproduced.

## 1. Introduction

The behaviour of cement-based materials at the early-age and its evolution with cement hydration is complex, due to the progressive hardening of the material and numerous couplings between mechanical, hydraulic, thermal and chemical aspects. A good understanding of this behaviour is necessary in various industrial applications, from the integrity assessment of cement-sheath in deep wells [1] to more recent construction technologies, such as 3D printing with cement-based materials.

The term “early-age” is to be understood here as “right after mixing”, at the very beginning of the fluid-to-solid transition and until the hardened state is reached. It therefore concerns a time frame prior to that which is often referred to in the literature with the expression “early-age”.

In the structure of deep wells (oil & gas, geothermal, CO₂ storage) a cement slurry is used to fill the annular space between the metal tubular (casing) and the rock formation. It forms a cement sheath that provides zonal isolation of different fluids along the well, as well as mechanical support and protection against corrosion for the casing. Loss of integrity of the cement sheath can lead to leakage and compromise the stability of an oil well. Early-age shrinkage of the cement paste linked with the decrease in pore pressure with hydration is believed to be a major factor in the leakage of oil wells [2–4]. In addition, pressure testing of the casing, which can be performed at pressures up to 80 MPa for high-pressure reservoirs, can lead to irreversible deformations of the cement

sheath, creating voids in the annulus where fluids and gas can migrate. The poromechanical behaviour of fresh cement paste and its evolution with hydration are therefore of prime interest for the oil industry.

The present work is the result of extensive experimental investigations of the poromechanical behaviour of oil-well cement from the early-age to the hardened state, along with micromechanical modelling and macroscale constitutive models. Ghabezloo et al. [5] first investigated the poromechanical behaviour of hardened class G cement paste using isotropic compression experiments, allowing to determine poromechanical parameters as well as showing evidence of a degradation of elastic properties stemming from microcracking. This work was extended to temperature effects by studying the thermal pressurisation of cement paste [6]. Vu et al. [7] revisited experimental results from [5] in the framework of viscoplasticity, providing a 1D model capable of capturing the time-dependent behaviour of the hardened paste and showing its temperature dependency. More recently, other researchers such as Lima et al. [8] provided additional tests on the effects of confining pressure on the stress–strain response of hardened oil well cement paste, namely confirming the transition to a more ductile behaviour at higher confining pressures. Similar methodologies were applied to hydrating cement paste by Agofack et al. [9,10], who investigated the poromechanical behaviour of class G cement paste starting right after mixing until the hardened state, bringing to light the effects of stress level, temperature and loading history on the paste’s behaviour at early ages. In addition, a predictive model based on hydration kinetics for estimating the time of creation of residual stresses (CRS) was formulated

<footnote>* Corresponding author.
E-mail address: siavash.ghabezloo@enpc.fr (S. Ghabezloo).</footnote>
List of symbols

<table><tr><td>A</td><td>Set of reactants</td></tr><tr><td>A<sub>s</sub></td><td>Set of solid reactants</td></tr><tr><td>σ</td><td>Stress tensor (Pa)</td></tr><tr><td>ε</td><td>Strain tensor</td></tr><tr><td>T</td><td>Temperature (K)</td></tr><tr><td>S</td><td>Lagrangian density of entropy (J·K<sup>-1</sup>·m<sup>-3</sup>)</td></tr><tr><td>Ψ</td><td>Lagrangian density of Helmholtz free energy (J·m<sup>-3</sup>)</td></tr><tr><td>μ<sub>α</sub></td><td>Chemical potential of species α (J·kg<sup>-1</sup>)</td></tr><tr><td>M<sub>α</sub></td><td>Lagrangian mass flow vector of species α (kg·m<sup>-2</sup>·s<sup>-1</sup>)</td></tr><tr><td>m<sub>α</sub></td><td>Lagrangian density of mass of species α (kg·m<sup>-3</sup>)</td></tr><tr><td>ρ<sub>α</sub></td><td>Density of species α (kg·m<sup>-3</sup>)</td></tr><tr><td>s<sub>α</sub></td><td>Entropy per unit of mass of species α (J·K<sup>-1</sup>·kg<sup>-1</sup>)</td></tr><tr><td>φ<sub>α</sub></td><td>Volume fraction of species α</td></tr><tr><td>p<sub>α</sub></td><td>Pressure of species α (Pa)</td></tr><tr><td>S<sub>s</sub></td><td>Lagrangian density of solid entropy (J·K<sup>-1</sup>·m<sup>-3</sup>)</td></tr><tr><td>ν<sub>α,i</sub></td><td>Stoichiometric coefficient for species α in reaction i</td></tr><tr><td>M<sub>α</sub></td><td>Molar mass of species α (kg·mol<sup>-1</sup>)</td></tr><tr><td>ε<sub>i</sub></td><td>Lagrangian density of extent of reaction i (mol·m<sup>-3</sup>)</td></tr><tr><td>φ</td><td>Porosity</td></tr><tr><td>φ<sub>0</sub></td><td>Initial porosity</td></tr><tr><td>φ</td><td>Mechanical contribution to porosity variation (φ − φ<sub>0</sub>)</td></tr><tr><td>S<sub>i</sub></td><td>Lagrangian saturation degree (i ∈ {l, g})</td></tr><tr><td>ν<sub>s</sub></td><td>Solid volume fraction</td></tr><tr><td>ΔV<sub>s,i</sub></td><td>Molar solid volume change in reaction i ()</td></tr><tr><td>p<sub>c</sub></td><td>Capillary pressure (Pa)</td></tr><tr><td>p<sup>*</sup></td><td>Equivalent pore pressure (Pa)</td></tr><tr><td>U</td><td>Lagrangian density of interface surface energy per unit of initial porosity (J·m<sup>-3</sup>)</td></tr><tr><td>□<sup>ir</sup></td><td>Irreversible part of variable □</td></tr><tr><td>C</td><td>Stiffness tensor (Pa)</td></tr><tr><td>R</td><td>Lagrangian density of recoverable energy (J·m<sup>-3</sup>)</td></tr><tr><td>L</td><td>Lagrangian density of locked energy (J·m<sup>-3</sup>)</td></tr><tr><td>b</td><td>Biot's effective stress coefficient</td></tr><tr><td>b<sub>i</sub></td><td>Partial Biot coefficient (i ∈ {l, g})</td></tr><tr><td>α<sub>s</sub></td><td>Solid thermal expansion (K<sup>-1</sup>)</td></tr><tr><td>L</td><td>Latent heat of hydration (J·mol<sup>-1</sup>)</td></tr><tr><td>C<sub>σs</sub></td><td>Volumetric heat capacity of the solid at constant stress (J·K<sup>-1</sup>·m<sup>-3</sup>)</td></tr><tr><td>V̄<sub>α</sub></td><td>Molar volume of species α (m<sup>3</sup>·mol<sup>-1</sup>)</td></tr><tr><td>A</td><td>Chemical affinity (J·mol<sup>-1</sup>)</td></tr><tr><td>ξ<sub>∞</sub></td><td>Theoretical maximum Lagrangian density of hydration extent (mol·m<sup>-3</sup>)</td></tr><tr><td>α</td><td>Hydration degree (ξ/ξ<sub>∞</sub>)</td></tr><tr><td>R<sub>H</sub></td><td>Relative humidity</td></tr><tr><td>σ'</td><td>Mean effective stress (Pa)</td></tr><tr><td>τ</td><td>Deviatoric stress (Pa)</td></tr><tr><td>σ<sub>c</sub></td><td>Yield stress in compression (Pa)</td></tr></table>

<table><tr><td>σ<sub>i</sub></td><td>Yield stress in traction (Pa)</td></tr><tr><td>k<sub>i</sub></td><td>Liquid permeability (m<sup>2</sup>)</td></tr><tr><td>ρ<sub>w</sub></td><td>Water density (kg·m<sup>-3</sup>)</td></tr><tr><td>η<sub>w</sub></td><td>Water dynamic viscosity (Pa·s)</td></tr><tr><td>k<sub>T</sub></td><td>Thermal conductivity (W·m<sup>-1</sup>·K<sup>-1</sup>)</td></tr><tr><td>E</td><td>Young's modulus of cement paste (Pa)</td></tr><tr><td>K</td><td>Bulk modulus of cement paste (Pa)</td></tr><tr><td>E<sub>g</sub></td><td>Chemical stiffening modulus (Pa)</td></tr><tr><td>h<sub>g</sub></td><td>Chemical hardening modulus (Pa)</td></tr><tr><td>h<sub>ε</sub></td><td>Isotropic hardening coefficient</td></tr></table>

and validated against experiments. Ghabezloo [11] proposed a two-scale micromechanical model to estimate the poroelastic parameters of hardened cement paste of different water-to-cement ratios. It was extended to different cement compositions using a volume fraction model [12]. Using the same micromechanical framework as for the poroelastic parameters, Ghabezloo [13] calculated the macroscopical thermal properties of cement pastes based on their water-to-cement ratio, extrapolating previously acquired experimental results to other mix designs.

Rich of these experimental results and micromechanical models, a comprehensive macromechanical model encompassing the evolution of the poro-elastoplastic behaviour of fresh cement pastes with hydration was developed by Samudio [14]. A similar approach by Agofack et al. [1] was proved able to explain both the early shrinkage and stress-strain behaviour of early-age class G cement paste hydrated under high pressures. In recent years, other research groups also presented modelling approaches based on poromechanics aiming at understanding early-age shrinkage and integrity of the cement sheath in deep wells. Sasaki et al. [15] looked at water migration from the rock formation towards the fresh cement paste using a simple hydromechanical model taking into account consumption of water by the hydration reaction through a sink rate fitted on experimental data. This model allowed to determine a threshold permeability for which water from the surrounding rock formation can seep into the cement sheath and reduce its shrinkage. Meng et al. [16] also proposed a model based on poromechanics, considering in addition the thermal aspect, in order to predict the integrity of the cement sheath. The model relies on experimental values of the stress drop at the interphase with the steel casing as input for simulations, which then check for stress-based criteria for debonding and cracking. Neither of these models introduce a variable related with chemical processes explicitly, relying instead on proxy measures of stress drop or volume change during the curing phase, which limits their generalisability, in particular with respect to temperature which is subject to vary significantly along the vertical profile of the well. Petersen and Ulm [17,18] have looked at cement sheath integrity in the framework of fracture mechanics, using a chemo-poromechanical model to describe the cement paste. Poroelastic parameters are evaluated using a three-level self-consistent homogenisation scheme. The model relies on expressions for hydration kinetics as well as water consumption from the hydration reaction, capturing the underlying physics behind the stress and pressure evolutions in the annulus. However, this model only accounts for the case of a water-saturated cement paste, although cement pastes can reach partial saturation in undrained conditions even under loads of several MPa, as evidence by Agofack et al. [1].

The thermo-poromechanical approach developed herein lies in the continuity of the work of Coussy on general poromechanics [19,20] and of Ulm and Coussy on cement-based materials [21,22]. This initial approach to cement-based materials introduced the idea of chemically hardening plasticity, but did not incorporate the fluid explicitly, taking into account cement shrinkage as an additional inelastic strain. Follow-up applications to other fields in cement and concrete research, such
as sprayed concrete (commonly known as shotcrete), can be found in studies by Sercombe et al. [23] and Hellmich et al. [24], for instance. They applied the chemomechanical hardening plasticity framework, enriching it with viscosity to model creep of shotcrete supporting tunnel walls. The porous nature of the material is hidden in this model, as pore pressure is not taken into account as a state variable. It is however paramount to understanding the behaviour of oil well cement, as water transport from the surrounding rock formation is expected to greatly impact shrinkage [15]. Water transport is also a prime source of concern in other applications such as 3D printing, with drying having negative actions on both shrinkage and development of hardened material properties [25].

This work further extends the above chemomechanical formalism to a strong coupling involving all physics (chemistry, hydraulics, thermals and mechanics), unifying aspects of different previous models under a unique, general theoretical framework for solidifying porous media. Not limited to cement-based materials, the theoretical part of this paper could also address other similar material such as geopolymers, for instance. A careful derivation of the model from thermodynamic principles is conducted, taking into account an arbitrary number of different chemical reactions and associated mechanisms. This paves the way towards better modelling of materials with complex chemistries such as cement-based materials incorporating accelerating admixtures. New to this study, also, is the proof of consistency of the poromechanical state equations when extended before the time of solid skeleton percolation.

Most of the ingredients of the model are taken over from Agofack et al. [1], including the determination of poroelastic parameters from a two level homogenisation process. However, substantial improvements have been carried out. Equations have been re-written in order to make their origins in the chemical reactions more apparent, resulting in a parameter set which can be mostly determined from physicochemical properties of the different phases. The formulation has been extended from poromechanical to thermo-poromechanical by considering the temperature as a field variable. The early-age creep of the material is also taken into account based on the viscoelastic model proposed by Carol and Bažant [26]. The plasticity model, based on the Cam-Clay yield surface, has also been improved with different hardening laws and the introduction of viscoplasticity. Account of the partially saturated state was improved by introducing a dependency of the water retention curve (WRC) to hydration. Finally, the model is presented as an open system rather than a closed one, adding classical transport laws for fluid and heat along with their evolution with certain internal variables according to relations from the literature.

After a detailed thermodynamics-based derivation of the material model with a description of all its components, the determination of model parameters for the class G cement paste of interest in the remaining of the study is presented. Then, applications to the simulation of several oedometric experiments from Agofack et al. [9] on class G cement paste in the first 144 h of hydration at 7 °C are shown. In particular, the complex undrained oedometric tests designed to investigate the effect of loading at the early age on the poromechanical response at later ages are successfully reproduced numerically for the first time. This shows the ability of the model, which relies on simple relationships in each of its unit components and minimal fitting, to capture the complex hardening mechanisms at play in hydrating cement paste. Following these results, a short discussion on current limitations and future perspectives is given before some concluding remarks.

## 2. Theoretical framework

The present model is an extension to chemically active media of the work of Coussy [19,20] on thermo-poromechanics, from which we borrow the core of the notation. Readers are referred to the aforementioned literature for more details.

The following description of cementitious materials is based on the mechanics of porous solids, therefore assuming the existence of a

connected solid skeleton. This assumption does not, in general, hold for materials at the very early age. Prior to the solid skeleton's percolation, that is to say, prior to the solid forming a connected network, the poromechanics description cannot hold. In the following section, we will first consider a percolated system. The constitutive equations for nonpercolated systems as well as the definition of the percolation threshold will be discussed thereafter in Section 2.8.

### 2.1. Basic coupled Thermo-Hydro-Mechanical (THM) framework with internal chemical reaction

Let us consider a porous medium in which chemical reactions occur: let A be the set of all reactants and products involved in these reactions, and A_s the subset of A comprised of only solid phases (namely clinker components and hydration products). With l, v and a representing the liquid, vapour and air respectively, A_s can be rewritten as A_s = A \ {l, v, a}.<fnref n="1" /> In the following, ≡ designates the Cauchy stress tensor (with the convention of positive stresses in tension), ε the linearised strain tensor, T the absolute temperature, S the Lagrangian density of entropy of the system and Ψ its Lagrangian density of Helmholtz free energy per unit of initial volume dΩ₀ of porous material. For any species α ∈ A, let M_α be its mass flow vector and μ_α its chemical potential.

#### 2.1.1. Dissipation of the solid

As per [19], Clausius-Duhem's inequality regarding the combined dissipation of the solid and the chemical reaction Φ_s→ gives:

$$ \Phi_{s \to} = \bar{\sigma} : \frac{d\bar{\epsilon}}{dt} - \sum_{a \in A} \mu_a \bar{\nabla}_a \cdot \bar{M}_a - S \frac{dT}{dt} - \frac{d\Psi}{dt} \geq 0. \quad (1) $$

Let us now consider N chemical reactions occurring in the material. Let $\dot{m}_{\to a,i}$ be the mass variation of species α ∈ A in reaction $i \in [1:N]$, and $\dot{m}_{\to a} = \sum_{i \in [1:N]} \dot{m}_{\to a,i}$ the total mass variation of species α due to all chemical processes. The mass conservation equation for species α is:

$$ \frac{dm_a}{dt} = -\nabla_a \cdot \bar{M}_a + \dot{m}_{\to a}. \quad (2) $$

Replacing (2) into (1):

$$ \Phi_{s \to} = \bar{\sigma} : \frac{d\bar{\epsilon}}{dt} + \sum_{a \in A} \mu_a \left( \frac{d\dot{m}_a}{dt} - \dot{m}_{\to a} \right) - S \frac{dT}{dt} - \frac{d\Psi}{dt} \geq 0. \quad (3) $$

The dissipation associated with chemical reactions, namely $\Phi_{\to} = -\sum \mu_a \dot{m}_{\to a}$ must satisfy the inequality $\Phi_{\to} \geq 0$ independently from the skeleton dissipation $\Phi_s = \Phi_{s \to} - \Phi_{\to}$:

$$ \Phi_s = \bar{\sigma} : \frac{d\bar{\epsilon}}{dt} + \sum_{a \in A} \mu_a \frac{d\dot{m}_a}{dt} - S \frac{dT}{dt} - \frac{d\Psi}{dt} \geq 0. \quad (4) $$

As we are interested in the behaviour of the solid skeleton, we must separate its free energy from that of the fluid phases occupying the pore space. Let us now remove the free energy of the fluid, composed of the liquid, vapour and air, and introduce the free energy $\Psi_s$ relative to the solid skeleton only:

$$ \Psi_s = \Psi - m_l \Psi_l - m_v \Psi_v - m_a \Psi_a. \quad (5) $$

with $m_f$ the Lagrangian density of mass and $\Psi_f$ the free energy per unit of mass of the different fluid constituents $f \in \{l, v, a\}$. We recall the following state equations for the Helmholtz free energy of a fluid:

$$ \psi_f = \psi_f \left( \frac{1}{\rho_f}, T \right), \quad p_f = - \frac{\partial \psi_f}{\partial \left( \frac{1}{\rho_f} \right)}, \quad s_f = - \frac{\partial \psi_f}{\partial T}, \quad (6) $$

where $\rho_f$ is the fluid's density, $p_f$ its pressure and $s_f$ its specific entropy. We would like now to express the Helmholtz free energy of the

<footnote n="1">¹ A \ {l, v, a} = {x \in A | x \notin \{l, v, a\}}</footnote>
fluid in terms of certain external variables of interest, in particular pressure and temperature, and other state variables such as the porosity, entropy and chemical potential, so as to make the different $\psi_f$ vanish from our expression of the dissipation. With $\phi_f$ as the volume fraction occupied by the different fluid constituents $f \in l, v, a$, differentiation of $m_f \psi_f$ using the above relations (6) along with $\mu_f = \psi_f + p_f / \rho_f$ gives:

$$ 
\begin{aligned}
d(m_f \psi_f) &= \psi_f dm_f + m_f \frac{p_f}{\rho_f^2} d\rho_f - m_f s_f dT \\
&= \psi_f dm_f + \phi_f \frac{p_f}{\rho_f} d\rho_f - m_f s_f dT \\
&= \left( \psi_f + \frac{p_f}{\rho_f} \right) dm_f - p_f d\phi_f - m_f s_f dT \\
&= \mu_f dm_f - p_f d\phi_f - m_f s_f dT,
\end{aligned} 
\quad (7) $$

which holds for the liquid, vapour and air.

In the air-vapour gas mixture, denoted by the subscript 'g', air and vapour occupy the same volume, namely $\phi_v = \phi_a = \phi_g$, and the pressure of the mixture is $p_g = p_v + p_a$. The separation between solid skeleton and fluid phases also allows for the introduction of the entropy of the solid skeleton only, $S_s = S - m_f s_l - m_v s_v - m_a s_a$, which naturally appears as the dual variable from the temperature $T$. Indeed, introducing it in (4) yields the following where the specific entropies of the fluid phases vanish:

$$ \Phi_s = \sigma : \frac{d\epsilon}{dt} + p_g \frac{d\phi_g}{dt} + p_l \frac{d\phi_l}{dt} + \sum_{\alpha \in A_s} \mu_\alpha \frac{dm_\alpha}{dt} - S_s \frac{dT}{dt} - \frac{d\Psi_s}{dt} \geq 0. \quad (8) $$

### 2.1.2. Mass balance in the chemical reactions

We now want to link the mass evolution of the different phases to their chemical and transport origins. Since there is no overall mass creation in any of the reactions, for all $i \in [1; N]$:

$$ \sum_{\alpha \in A} \stackrel{\circ}{m}_{\rightarrow \alpha, i} = 0. \quad (9) $$

With $v_{\alpha,i}$ being the stoichiometric coefficient for species $\alpha \in A$ in reaction $i \in [1; N]$, positive for products and negative for reactants, and $M_\alpha$ its molar mass. The stoichiometry of the reaction imposes that for all $\alpha, \beta \in A$:

$$ \frac{\stackrel{\circ}{m}_{\rightarrow \alpha, i}}{v_{\alpha,i} M_{\alpha}} = \frac{\stackrel{\circ}{m}_{\rightarrow \beta, i}}{v_{\beta,i} M_{\beta}}. \quad (10) $$

Introducing the reaction rates $(\dot{\xi}_i)_{i \in [1:N]}$ as the above ratio for each reaction

$$ \stackrel{\circ}{\xi}_i = \frac{\stackrel{\circ}{m}_{\rightarrow \alpha, i}}{v_{\alpha,i} M_{\alpha}} \quad (11) $$

allows to rewrite the mass conservation (2):

$$ \frac{dm_\alpha}{dt} = -\nabla \cdot \mathbf{M}_\alpha + \sum_{i=1}^{N} v_{\alpha,i} M_\alpha \stackrel{\circ}{\xi}_i. \quad (12) $$

Assuming that there is no transport of species $\alpha$ except for the fluid ($\{l, v, a\}$), for all species $\alpha \in A_s$:

$$ \frac{dm_\alpha}{dt} = \sum_{i=1}^{N} v_{\alpha,i} M_\alpha \stackrel{\circ}{\xi}_i \quad (13) $$

indicating that the rates of reactions $(\dot{\xi}_i)$ are in fact time derivatives:

$$ \stackrel{\circ}{\xi}_i = \frac{d\xi_i}{dt} \quad (14) $$

with $\xi_i$ the extent of the $i$th reaction. We can then express the evolution of the mass of the constituents in terms of the advancement of the chemical processes, for $\alpha \in A_s$:

$$ \frac{dm_\alpha}{dt} = \sum_{i=1}^{N} v_{\alpha,i} M_\alpha \frac{d\xi_i}{dt}, \quad (15) $$

which allows to rewrite (8) as:

$$ \Phi = \sigma : \frac{d\epsilon}{dt} + p_g \frac{d\phi_g}{dt} + p_l \frac{d\phi_l}{dt} - \sum_{i=1}^{N} A_{s,i} \frac{d\xi_i}{dt} - S_s \frac{dT}{dt} - \frac{d\Psi_s}{dt} \geq 0, \quad (16) $$

where

$$ A_{s,i} = - \sum_{\alpha \in A_s} v_{\alpha,i} M_\alpha \mu_\alpha \quad (17) $$

is reminiscent of the expression of de Donder's chemical affinity of a reaction, restricted to the solid constituents.

### 2.1.3. Consideration of the unsaturated state

The current total pore volume $\phi_g d\Omega_0$ is filled with the liquid and gas phases, whose volumes are $\phi_l d\Omega_0$ and $\phi_g d\Omega_0$, respectively. Starting from an initial porosity $\phi_0$, deformations of the pore space induced by mechanical loading, pore pressure, and thermal deformations result in a porosity variation $\phi$, so that $\phi = \phi_0 + \phi$. Now, in the same way, the liquid and gas partial porosities can be split into two terms, one related to the reference, undeformed state, and the other related to the deformation of the partial pore volume due to the action of this phase's pressure. The so-called Lagrangian saturation $S_l$ of phase $i$ ($i \in \{l, g\}$) is defined as the fraction of the pore space that a phase $i$ would occupy in the undeformed state. The partial porosities can then be expressed as:

$$ \phi_i = S_i \phi_0 + \phi_l, \quad (18) $$

where $\phi_g + \phi_l = \phi$. Since $\phi_g + \phi_l = \phi$, we have $\phi = (S_g + S_l)\phi_0 + \phi$, so that ultimately $S_g + S_l = 1$.

The particularity of a chemically active porous medium such as cement paste lies in the redefinition of the initial state through the creation of additional solid volume in hydration reactions. Indeed, as hydration progresses, the solid volume fraction will increase and thus the pore volume will decrease, so that the initial porosity $\phi_0$ is in fact a function of $(\xi_i)_{i \in [1:N]}$. The change in solid molar volume in hydration reaction number $i$ is denoted by $\Delta\bar{V}_{s,i}$, and $v_s = 1 - \phi_0$ is the solid volume fraction. The initial porosity increment then reads:

$$ d\phi_0 = -dv_s = - \sum_{i=1}^{N} \Delta\bar{V}_{s,i} d\xi_i. \quad (19) $$

Using (18) and (19), we rewrite the porosity-related terms in inequality (16) so as to make the liquid saturation degree $S_l$ appear:

$$ 
\begin{aligned}
& p_g \frac{d\phi_g}{dt} + p_l \frac{d\phi_l}{dt} \\
&= p_g \frac{d}{dt} (S_g \phi_0 + \phi_g) + p_l \frac{d}{dt} (S_l \phi_0 + \phi_l) \\
&= p_g \frac{d\phi_g}{dt} + p_l \frac{d\phi_l}{dt} - \phi_0 (p_g - p_l) \frac{dS_l}{dt} - (S_l p_l + S_g p_g) \sum_{i=1}^{N} \Delta\bar{V}_{s,i} \frac{d\xi_i}{dt}
\end{aligned} 
\quad (20) $$

We can see here two different pressure quantities appearing: the capillary pressure $p_c$, and an equivalent pore pressure $p^*$, defined as:

$$ p_c = p_g - p_l, \quad (21) $$

$$ p^* = S_l p_l + S_g p_g. \quad (22) $$

In the saturated case, $p^*$ coincides with $p_l$.

Rewriting the dissipation, we get:

$$ 
\begin{aligned}
\Phi = \sigma : \frac{d\epsilon}{dt} + p_g \frac{d\phi_g}{dt} + p_l \frac{d\phi_l}{dt} - \phi_0 p_c \frac{dS_l}{dt} - p^* \sum_{i=1}^{N} \Delta\bar{V}_{s,i} \frac{d\xi_i}{dt} \\
- \sum_{i=1}^{N} A_{s,i} \frac{d\xi_i}{dt} - S_s \frac{dT}{dt} - \frac{d\Psi_s}{dt} \geq 0.
\end{aligned} 
\quad (23) $$

The free energy $\Psi_s$ is however not the free energy of the solid as it still includes the surface energy at the interface between the fluid and gas phases in the unsaturated case. Introducing $U$ as the interface
surface energy per unit of initial pore volume, the free energy of the solid $\psi_s$ is thus:

$$ \psi_s = \Psi_s - \phi_0 U. \quad (24) $$

This surface energy $U$ depends on the saturation degree $S_l$ and also on the pore structure. The pore structure can be altered through loading and through the chemical process; however, the influence of the former will be disregarded, as it is negligible. We also choose to neglect a potential direct influence of the temperature $T$ on the surface energy (see Section 5.3 for a more detailed discussion). Consequently, $U = U(S_l, \xi_l)$. Ultimately, we rewrite the dissipation in terms of the new free energy $\psi_s$:

$$ \Phi = \sigma : \frac{d\epsilon}{dt} + p_g \frac{d\phi_g}{dt} + p_l \frac{d\phi_l}{dt} - \phi_0 \left( p_c + \frac{\partial U}{\partial S_l} \right) \frac{dS_l}{dt} + \sum_{i=1}^{N} \frac{\partial}{\partial \xi_i} (\phi_0(p^* - U)) \frac{d\xi_i}{dt} - \sum_{i=1}^{N} A_{s,i} \frac{d\xi_i}{dt} - S_s \frac{dT}{dt} - \frac{d\psi_s}{dt} \geq 0. \quad (25) $$

As it would be more convenient to have the pressures as variables rather than the porosity variations, we use the free energy $G_s = \psi_s - p_g\phi_g - p_l\phi_l$ to rewrite the dissipation:

$$ \Phi = \sigma : \frac{d\epsilon}{dt} - \phi_g \frac{dp_g}{dt} - \phi_l \frac{dp_l}{dt} - \phi_0 \left( p_c + \frac{\partial U}{\partial S_l} \right) \frac{dS_l}{dt} + \sum_{i=1}^{N} \left[ \frac{\partial}{\partial \xi_i} (\phi_0(p^* - U)) - A_{s,i} \right] \frac{d\xi_i}{dt} - S_s \frac{dT}{dt} - \frac{dG_s}{dt} \geq 0. \quad (26) $$

From there, one would classically derive the constitutive thermo-poroelastic relations by considering no dissipation and by virtue of the mutual independence of the thermodynamic variables. However, as detailed in the next section, this cannot be applied to the case of a chemically solidifying medium, as some energy is locked inside the skeleton in the solidification process, even if the skeleton is otherwise perfectly elastic. This is due to the fact that the thermoporoelastic coefficients are dependent on the solidification process. Namely the stiffness tensor $C$, Biot's coefficient $b$ and Biot's modulus $N$ will be seen to depend directly on the chemical variables, while other coefficients (stress and porosity variations with respect to temperature as well as apparent heat capacity of the solid skeleton) depend on the former and on the current “initial” porosity $\phi_0$.

## 2.2. Chemo-mechanics of solidifying media

### 2.2.1. Recoverable and locked energies

Consider a simple case of a chemically active solid with energy $F$. The Clausius-Duhem inequality gives the following:

$$ \sigma : \frac{d\epsilon}{dt} + A_s \frac{d\epsilon}{dt} - \frac{dF}{dt} \geq 0. \quad (27) $$

If this solid is purely elastic and no other dissipation occurs, then the inequality becomes an equality. Furthermore, if the chemical reaction is a solidification reaction, a higher extent of reaction will lead to an increase in solid mechanical stiffness. Now consider the succession of a mechanical loading at a given stress state, followed by an increase in the reaction's extent, and then unloading to the previous stress state (Fig. 1). Since the additional solid is created in the current configuration by the chemical reaction, a reasonable assumption is to consider the new solid to be free from stresses and strains, thus having no impact on the stress-strain state of the solid. As the stiffness of the solid has increased between the loading and unloading phases, the unloading slope is steeper, causing the solid to not recover its initial configuration, but rather keep an irreversible strain $\epsilon^{ir}$. This means that a part of the elastic energy is chemically locked and cannot be recovered unless the chemical reaction can be reversed.

Accordingly, the solid's energy can be split in recoverable and chemically locked parts:

$$ F = R(\underline{\epsilon} - \underline{\epsilon}^{ir}, \xi) + L(\underline{\epsilon}^{ir}, \xi) \quad (28) $$

<figure><img src="image_4.png" /><figcaption>(a) Loading as a function of hydration.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Stress-strain response.</figcaption></figure>

<caption>Fig. 1. Illustration of irreversible strains.</caption>

where $R$ is the usual quadratic form:

$$ R = \frac{1}{2} (\underline{\epsilon} - \underline{\epsilon}^{ir}) : C(\underline{\xi}) : (\underline{\epsilon} - \underline{\epsilon}^{ir}). \quad (29) $$

Assuming that there are no dissipative mechanisms in the solid, the state law reads:

$$ \sigma = \frac{\partial R}{\partial (\underline{\epsilon} - \underline{\epsilon}^{ir})} = C(\underline{\xi}) : (\underline{\epsilon} - \underline{\epsilon}^{ir}), \quad (30) $$

which, in incremental form, gives:

$$ d\sigma = C(\underline{\xi}) : (d\underline{\epsilon} - d\underline{\epsilon}^{ir}) + \frac{dC}{d\underline{\xi}} : (\underline{\epsilon} - \underline{\epsilon}^{ir})d\underline{\xi}. \quad (31) $$

The hypothesis of new solid created free of stresses and strains in the current configuration implies that regardless of the chemical reaction's extent increment $d\underline{\xi}$, in a transformation where $d\underline{\epsilon} = 0$, there should be no increment of stress, which leads to the following expression for the evolution of the irreversible strains:

$$ d\underline{\epsilon}^{ir} = C^{-1}(\underline{\xi}) : \frac{dC}{d\underline{\xi}} : (\underline{\epsilon} - \underline{\epsilon}^{ir})d\underline{\xi}. \quad (32) $$

In fine, the incremental state equation becomes the usual elastic one:

$$ d\sigma = C(\underline{\xi}) : d\underline{\epsilon}, \quad (33) $$

with a stiffness tensor that is now dependent on the extent $\underline{\xi}$ of the chemical reaction.
### 2.2.2. Limitations

Caution should, however, be exercised when using the aforementioned assumption of new solid being created with a known stress state. Numerical integration using such a model should avoid resolution over large time steps, as the chronology of mechanical loading relative to that of the chemical reaction may be lost. A general condition that would ensure the validity of this assumption would be if the characteristic time $\tau_m$ of mechanical processes is greater than that $\tau_{ch}$ of the chemical processes: $\tau_m \gg \tau_{ch}$. At the other extreme, if the mechanical loading is close to instantaneous, one could consider the solid as non-reactive during loading. In practice, in the case of small monotonic perturbations, the error induced by the approximate chronology of events, considering sufficiently small time steps, is rather limited. Nevertheless, one should be wary of the practical application in the case of large and non-monotonic loading concomitant to a solidifying reaction.

It should also be noted that although $d\xi \ge 0$ has not been explicitly enforced here, it is advisable not to use this formulation in a case where both precipitation and dissolution could occur. This solidification theory is based on the fact that, at the microscopic scale, the stress is not uniform, since newly precipitated elements of solid are stress-free in an otherwise potentially stressed solid. This heterogeneous stress distribution, unknown at the macroscale, is then problematic in the case of dissolution, as the previous stress state of the dissolved solid is not known. The stress redistribution due to the dissolution of part of the solid skeleton is therefore hard to approach macroscopically. From the point of view of the chemical potentials of dissolution, solid undergoing the highest mechanical stress should dissolve first due to their higher potential. However, this information is not accessible at the macroscale. In contrast to solidification, which does not alter the stress state when other variables are kept constant, in the case of dissolution, the stress state at the REV scale may thus be altered by the chemical process, by virtue of removing a stressed solid. As accounting for heterogeneities in the stress state of the solid matrix requires modelling at a significantly smaller scale, combined precipitation and dissolution processes should not be considered within the scope of this theoretical framework.

### 2.3. Constitutive equations for the solid skeleton

For the sake of simplicity, we will detail here the constitutive equations taking into account only one hydration reaction with associated extent $\xi$. Hydrating cements are usually described through a single extent of reaction (or equivalently, a unique degree of hydration $\alpha = \xi/\xi_{\infty}$, where $\xi_{\infty}$ stands for the ultimate extent corresponding to complete hydration of the cement). Although multiple hydration reactions take place simultaneously, this amounts to considering a homogenised fictitious hydration reaction. This point of view will be adopted hereafter, although in the authors' experience, some cement-based materials, namely mixes with significant amounts of accelerating admixtures, display hydration kinetics, which would be better captured using more than a single extent of reaction.

Extending the above Section 2.2.1 to the multiphysics framework presented before, we split the energy $G_i$ in a recoverable and locked part. Due to the hydration dependency of the thermoporomechanical coefficients, not only the strain tensor $\epsilon$ but also the pressures and temperature will be divided into reversible and irreversible parts. Indeed, in a similar fashion to Fig. 1, let us replace $\epsilon$ with $p_i$, $p_g$ or $T$ and $\sigma$ with $\varphi_i$, $\varphi_g$ or $S_s$ respectively. In a transformation where other thermodynamic variables than those of interest are kept constant, and in the hypothesis of linear poroelasticity, there is a linear relationship between the variable and its associated thermodynamic force. However, similarly to E in Fig. 1, the coefficient in this linear relationship is a function of the hydration extent $\xi$.

In the case of pressures, the relevant property is the partial Biot's modulus $N_{il}$, where $i$ is the liquid ($l$) or the gas ($g$), while in the case of temperature, it is the apparent heat capacity of the solid

skeleton, which is proportional to the solid fraction (see Eq. (48) below). Fig. 2 shows how, similarly to Fig. 1, the hydration dependency of these parameters induces irreversible parts of the related external state variable. Loading and then unloading after a change in hydration extent will result in different loading and unloading slopes. Recovering the initial value of the associated thermodynamic force therefore requires a different value of the variable from its initial value prior to loading. This constitutes irreversible pressures and temperature. The recoverable energy is hence of the form:

$$R(\underline{\epsilon} - \underline{\epsilon}^{ir}, p_g - p_g^{ir}, p_l - p_l^{ir}, T - T^{ir}, S_l, \xi), \quad (34)$$

while the locked energy is so that:

$$\mathcal{L} = \mathcal{L}(\underline{\epsilon}^{ir}, p_g^{ir}, p_l^{ir}, T^{ir}, S_l, \xi). \quad (35)$$

Rewriting Clausius-Duhem's inequality:

$$\left(\underline{\sigma} - \frac{\partial R}{\partial(\underline{\epsilon} - \underline{\epsilon}^{ir})}\right) : \frac{d(\underline{\epsilon} - \underline{\epsilon}^{ir})}{dt} - \left(\varphi_g + \frac{\partial R}{\partial(p_g - p_g^{ir})}\right) \frac{d(p_g - p_g^{ir})}{dt} - \left(\varphi_l + \frac{\partial R}{\partial(p_l - p_l^{ir})}\right) \frac{d(p_l - p_l^{ir})}{dt} \dots, \quad (36)$$

and using the rate independence of the energy along with the mutual independence of the state variables, their rates can be chosen arbitrarily and consequently require the coefficient in front of each rate to be equal to zero in order to ensure the positivity of the dissipation at all times. This results in the following state equations:

$$\begin{align}
\underline{\sigma} &= \frac{\partial R}{\partial(\underline{\epsilon} - \underline{\epsilon}^{ir})}, \nonumber \\
\varphi_g &= -\frac{\partial R}{\partial(p_g - p_g^{ir})}, \nonumber \\
\varphi_l &= -\frac{\partial R}{\partial(p_l - p_l^{ir})}, \nonumber \\
S_s &= -\frac{\partial R}{\partial(T - T^{ir})}, \nonumber \\
\phi_0 \left(p_c + \frac{\partial U}{\partial S_l}\right) &= \frac{\partial R}{\partial S_l} + \frac{\partial \mathcal{L}}{\partial S_l}, \nonumber \\
\frac{\partial}{\partial \xi} \left(\phi_0(p^* - U)\right) - \mathcal{A}_s &= \frac{\partial R}{\partial \xi} + \frac{\partial \mathcal{L}}{\partial \xi}. \tag{37}
\end{align}$$

Note that we also have:

$$\begin{align}
\underline{\sigma} &= \frac{\partial \mathcal{L}}{\partial \underline{\epsilon}^{ir}}, \nonumber \\
\varphi_g &= -\frac{\partial \mathcal{L}}{\partial p_g^{ir}}, \nonumber \\
\varphi_l &= -\frac{\partial \mathcal{L}}{\partial p_l^{ir}}. \tag{38}
\end{align}$$

The fifth equation in (37) is the complementary evolution law which determines the water retention curve, that is the relationship between capillary pressure and saturation degree. As it stands, this relationship depends on $\partial R / \partial S_l$. However, in linear thermoporoelasticity, $R$ is a quadratic function of $\underline{\epsilon}$, $p_g$, $p_l$ and $T$ and here also of $\xi$. The coefficients of this quadratic function are related to the current material properties, and thus, in our case, depend also on $S_l$, $\xi$, and potentially on $T$. As a consequence, as regards saturation, $\partial R / \partial S_l$ is the same quadratic function as $R$ with its coefficients replaced by their derivatives with respect to $S_l$, and are then second-order terms which will be neglected with respect to $\partial U / \partial S_l$. A similar assumption regarding the locked energy $\mathcal{L}$ allows to also neglect it with respect to $\partial U / \partial S_l$. Under this assumption, the capillary pressure $p_c$ is linked to the surface energy $U$ through:

$$p_c = -\frac{\partial U}{\partial S_l}. \quad (39)$$
<figure><img src="image_6.png" /><figcaption>Fig. 2. Illustration of the origin of irreversible pressures (a) and temperature (b).</figcaption></figure>

We can then establish the state equations. The mechanical state equation is as follows:

$$ \begin{aligned} \underline{d\sigma}|_{S_i} = \mathcal{C}(\xi) : \underline{d(\xi - \varepsilon^{ir})} - b_g(\xi, S_i) \underline{1d(p_g - p_g^{ir})} & \\ - b_l(\xi, S_i) \underline{1d(p_l - p_l^{ir})} - 3\kappa(\xi) \underline{1d(T - T^{ir})} + \underline{\Sigma d\xi}, & \end{aligned} \quad (40) $$

where $b_l$ and $b_g$ designate the partial Biot coefficients for liquid and gas respectively, $\kappa = K\alpha_s$ with K the bulk modulus and $\alpha_s$ the solid volumetric expansion coefficient. The term $\underline{\Sigma}$ is given by:

$$ \begin{aligned} \underline{\Sigma} = \underline{\Sigma_0} + \frac{d\mathcal{C}}{d\xi} : (\underline{\varepsilon} - \underline{\varepsilon^{ir}}) - \frac{\partial b_g}{\partial \xi}(p_g - p_g^{ir}) \underline{1} & \\ - \frac{\partial b_l}{\partial \xi}(p_l - p_l^{ir}) \underline{1} - 3 \frac{\partial \kappa}{\partial \xi}(T - T^{ir}) \underline{1} & \end{aligned} \quad (41) $$

A solely chemical increment results in no stress increment as per our previous assumption, and similarly to (32), we write:

$$ \mathcal{C}(\xi) : \underline{d\varepsilon^{ir}} - b_g(\xi, S_i) \underline{1d}p_g^{ir} - b_l(\xi, S_i) \underline{1d}p_l^{ir} - 3\kappa(\xi) \underline{1d}T^{ir} = \underline{\Sigma} \underline{d\xi}, \quad (42) $$

resulting in the following state equation:

$$ \underline{d\sigma}|_{S_i} = \mathcal{C}(\xi) : \underline{d\varepsilon} - b_g(\xi, S_i) \underline{1d}p_g - b_l(\xi, S_i) \underline{1d}p_l - 3\kappa(\xi) \underline{1d}T. \quad (43) $$

Making a classical assumption of isodeformation of partial porosities [20], partial Biot coefficients are proportional to each phase's saturation degree:

$$ b_l(\xi, S_i) = S_i b(\xi). \quad (44) $$

In a similar fashion, as a purely chemical increment is supposed to result in no stress-induced (mechanical) change of the porosities (chemical variation of the porosity being fully encompassed in the variation of the initial porosity $\phi_0$), state equations for the liquid and gas are the following:

$$ \begin{aligned} \underline{d\varphi_l}|_{S_i} &= b_l(\xi, S_i) \underline{1} : \underline{d\varepsilon} + \frac{\underline{d}p_l}{N_{ll}(\xi, S_i)} + \frac{\underline{d}p_g}{N_{lg}(\xi, S_i)} - 3\alpha_{\varphi_l}(\xi, S_i) \underline{d}T, & (45) \\ \underline{d\varphi_g}|_{S_i} &= b_g(\xi, S_i) \underline{1} : \underline{d\varepsilon} + \frac{\underline{d}p_l}{N_{gl}(\xi, S_i)} + \frac{\underline{d}p_g}{N_{gg}(\xi, S_i)} - 3\alpha_{\varphi_g}(\xi, S_i) \underline{d}T. & \end{aligned} $$

Under the same previous assumption, parameters relative to the evolution of porosities can be simplified as such:

$$ \begin{aligned} \frac{1}{N_{ii}} + \frac{1}{N_{lg}} &= S_i \frac{b(\xi) - \phi_0(\xi)}{K_s} & (46) \\ \alpha_{\varphi_l} &= \alpha_s [b(\xi) - \phi_0(\xi)] S_i & \end{aligned} $$

Then, we can combine the two state equations for the liquid and gas partial porosities into:

$$ \begin{aligned} \underline{d\varphi}|_{S_i} &= \underline{d\varphi_l}|_{S_i} + \underline{d\varphi_g}|_{S_i} \\ &= \underline{b(\xi)} \underline{1} : \underline{d\varepsilon} + \frac{\underline{b(\xi) - \phi_0(\xi)}}{K_S} (S_l \underline{d}p_l + S_g \underline{d}p_g) & (47) \\ & \quad - 3\alpha_s (\underline{b(\xi) - \phi_0(\xi)}) \underline{d}T. \end{aligned} $$

Due to the symmetries of the constitutive relations, the first two coefficients of the thermal equations are known from the previous. Only remains an “apparent heat capacity” which depends on the current solid volume fraction:

$$ \begin{aligned} \underline{dS_s}|_{S_i} = 3\kappa(\xi) \underline{1} : \underline{d\varepsilon} - 3\alpha_s(b(\xi) - \phi_0(\xi))(S_l \underline{d}p_l + S_g \underline{d}p_g) & (48) \\ + \frac{(1 - \phi_0(\xi))C_{\alpha_s}}{T} \underline{d}T, & \end{aligned} $$

where $C_{\alpha_s}$ is the solid's volumetric heat capacity at constant stress.

The last state equation relates to the so-called solid chemical affinity, that is, to the chemical potential of solid constituents. However, because the activities of solids are always equal to unity, this solid chemical affinity is actually a combination of the standard chemical potentials and, as such, can be expected to depend only on pressure and temperature. This term is therefore unable to explain the kinetics of hydration of cementitious materials: indeed, this kinetics depends on ionic concentrations in the pore fluid from the dissolution of clinker phases and not on solid phases which constitute the matrix. The kinetics of hydration is then driven by the pore fluid, and the last state equation for the solid alone is of little usefulness. The hydration kinetics is then to be modelled separately and is detailed in the next section as part of the complementary evolution laws.

## 2.4. Strain decomposition

In the present model, strain is decomposed in three distinct parts: an elastic part, a visco-elastic part and a viscoplastic part:

$$ \underline{\varepsilon} = \underline{\varepsilon}^e + \underline{\varepsilon}^v + \underline{\varepsilon}^p. \quad (49) $$

### 2.4.1. Viscoelastic strains

Cement-based materials are known to be viscous, on a smaller time scale, especially at the early age, or on long time scales for creep of hardened concrete. The focus here lies more towards hydrating cement; hence, the model does not take into account long-term creep, which could be added as a supplementary viscous cell. We make use here of a Kelvin-Voigt element, although, as for elasticity, we have to take into account ageing, that is, parameters which depend on the hydration degree. Similarly to Carol and Bazant [26], we consider ageing through a modulus that is proportional to the solid volume fraction $1 - \phi_0$ with
<figure><img src="image_7.png" /><figcaption>Fig. 3. One dimensional representation of the considered viscoplastic behaviour with ageing.</figcaption></figure>

a constant relaxation time $\tau_v$, so that by defining a reference stiffness $E_v$, the stiffness and viscosity of the Kelvin-Voigt model are $(1 - \phi_0)E_v$ and $\tau_v(1 - \phi_0)E_v$, respectively. We consider the viscous Poisson ratio $v_v$ to be constant for the sake of simplicity. The differential equation governing the evolution of the viscoelastic strain is given then by:

$$ \frac{d\sigma'}{dt} = \tau_v C_v : \frac{d^2\varepsilon^v}{dt^2} + \left(1 - \phi_0 - \tau_v \frac{d\phi_0}{dt}\right) C_v : \frac{d\varepsilon^v}{dt}, \quad (50) $$

where $\sigma' = \sigma + \beta_v p_1 \frac{1}{\bar{\varepsilon}}$ is the effective stress driving viscoelastic strains, $\beta_v$ the related effective stress coefficient and $C_v$ the viscoelastic stiffness tensor. $\beta_v = 1$ corresponds to a viscoelastically incompressible solid. The evolution of porosity with viscoelastic strains can be expressed as:

$$ \frac{\partial \varphi}{\partial \varepsilon^v} = \beta_v \frac{1}{\bar{\varepsilon}}. \quad (51) $$

### 2.4.2. Viscoplastic strains

The goal of this model is to provide a continuous description of the behaviour of cement-based materials from the early age to the hardened state. This implies having unique yield criterion and hardening laws that continuously vary with cement hydration capable of describing both the pasty behaviour of the early age as well as the rock-like behaviour of the hardened state. A Cam-Clay yield criterion is used, as it was previously found to be suitable for describing the early age behaviour of cement paste [1]:

$$ f(\sigma', \tau) = \tau^2 + M^2(\sigma' + \sigma_l)(\sigma' - \sigma_c), \quad (52) $$

where $\sigma' = \text{tr}(\underline{\sigma} + \beta_p p^* \frac{1}{\bar{\varepsilon}})$ is a mean effective stress with an associated effective stress coefficient $\beta_p$, $\tau = \sqrt{\frac{3}{2}\underline{\sigma}^D : \underline{\sigma}^D}$ where $\underline{\sigma}^D$ is the deviatoric stress tensor is the deviatoric stress, and $\bar{M}$, $\sigma_c$ and $\bar{\sigma}_t$ are material parameters.

Usual isotropic hardening plasticity can be represented in one dimension by a frictional slider in parallel with a spring. To introduce time dependency and, in particular, no instantaneous plastic strains, we elect to add a dashpot in parallel with the slider and spring. Fig. 3 illustrates the viscoplastic element, with ageing of the component with the hydration extent $\xi$. In addition to an isotropic hardening coefficient $h_\xi$, we define a plastic relaxation time $\tau_p$, so that the spring has a stiffness $h_\xi \sigma_c$ and the dashpot a viscosity $\tau_p h_\xi \sigma_c$. The related hardening rule can then be expressed as:

$$ d\sigma_c|_{\xi} = \sigma_c h_\xi (1 - \phi_0(\xi))(\text{tr}(d\bar{\varepsilon}^p) + \tau_p \text{tr}(d\bar{\varepsilon}^p)). \quad (53) $$

The hydration progress results in an expansion of the elastic domain. Considering that in an otherwise unstressed state, the increase in yield stress is directly related to the additional volume of solid constituents $v_s$, the chemical hardening can be expressed as:

$$ \frac{\partial \sigma_c}{\partial \xi} = -h_{\xi} \frac{d\phi_0}{d\xi} = h_{\xi} \Delta \bar{V}_s, \quad (54) $$

<figure><img src="image_8.png" /><figcaption>Fig. 4. Evolution of elastic domain with hydration.</figcaption></figure>

where $h_\xi$ is the yield stress increase provided by the creation of a unit volume fraction of solid hydrates.

The tensile pre-consolidation pressure $\sigma_t$ should also follow a similar hardening law, and so we will consider that it is at all times proportional to $\sigma_c$ through a coefficient $f_t$:

$$ \sigma_t = f_t \sigma_c. \quad (55) $$

The transition from a quasi-fluid state to a solid state also calls for a second hardening law. Indeed, the ratio between shear and compressive strength is bound to be much lower for a paste than for hardened concrete [27]. This can be described through a hardening law for the friction ratio $M$, where we introduce a parameter $M_\xi$ such that:

$$ dM = M_{\xi} \Delta \bar{V}_{\xi} d\xi. \quad (56) $$

This results in the evolution of the yield surface depicted in Fig. 4.

Plastic strains also induce plastic deformation of the pore space:

$$ \frac{\partial \varphi}{\partial \varepsilon^p} = \beta_p \frac{1}{\bar{\varepsilon}}. \quad (57) $$

which combined with the assumption of a plastically incompressible solid matrix leads to $\beta_p = 1$ as plastic deformation can only be the result of deformation of the pores [19].

## 2.5. Complementary evolution laws

### 2.5.1. Water retention curve

Eq. (39) implies the existence of a relationship between $p_c$, $S_l$ and $\xi$, since $U = U(S_l, \xi)$. $p_c$ is classically linked to the saturation degree $S_l$ through the water retention curve, for which a widely accepted model in soil mechanics is the one proposed by Van Genuchten [28]. In this model, $S_l$ can be related to capillary pressure $p_c$ (21) through:

$$ S_l = \left(1 + \left(\frac{p_c}{p_r}\right)^n\right)^{\frac{1-n}{n}}, \quad (58) $$

where $p_r$ is a reference pressure and $n > 1$ a real exponent. However, this description lacks a dependency to the chemical variable $\xi$: as the pore structure of cement-based materials evolves with hydration, the water retention curve should shift. We deal with this evolution by considering a reference pressure, characteristic of the pore structure, which depends on hydration. As in [14], a parameter $k_{sl}$ is introduced so that:

$$ p_r(\xi) = p_{r,0} \exp \left( k_{sl} \frac{\xi}{\xi_{\infty}} \right). \quad (59) $$

The resulting water retention curves are illustrated in Fig. 5 for $n = 1.48$, $p_{r,0} = 10$ kPa and $k_{sl} = 7$.
<figure><img src="image_9.png" /><figcaption>Fig. 5. Example evolution of the water retention curve with hydration degree $\alpha = \xi/\xi_{\infty}$.</figcaption></figure>

### 2.5.2. Hydration kinetics

The law governing the hydration reaction links the reaction rate $\dot{\xi}$ with the chemical affinity $A$. As the chemical affinity $A$ is not given by the state equations for the solid skeleton, inspired by [29] we adopt the hypothesis of a cubic function for the chemical contribution to the free energy, resulting in a quadratic dependence of $A$ to $\xi$ and write:

$$ \frac{\partial A}{\partial \xi} = k_{\xi} \left( \xi_{\infty} - \frac{A_0}{k_{\xi} \xi_{\infty}} - 2\xi \right). \qquad (60) $$

The chemical affinity of the reaction also depends on the relative humidity in the pore space, as it has been reported both from a theoretical and experimental standpoint that hydration reactions stop under a certain relative humidity (typically around 80% for alite) [30–32]. For this model, we will use the empirical formula proposed by Wyrzykowski et al. [32] for the reduced reaction rate under unsaturated conditions and apply it to the chemical affinity:

$$ \frac{A(R_H)}{A(R_H = 1)} = \min \left( 1, 3.7 (R_H)^{25} \right). \qquad (61) $$

The internal relative humidity can be related to the capillary pressure $p_c$ through Kelvin's equation [33]:

$$ R_H = \exp \left( - \frac{p_c M_w}{RT \rho_w} \right) \qquad (62) $$

Neglecting effects other than relative humidity on the chemical affinity, its expression becomes:

$$ A(\xi, R_H) = \min \left( 1, 3.7 (R_H)^{25} \right) k_{\xi} \left( \frac{A_0}{k_{\xi} \xi_{\infty}} + \xi \right) (\xi_{\infty} - \xi). \qquad (63) $$

Hydration kinetics of cement is known to depend both on temperature and pressure, with higher pressure and/or temperature leading to a higher hydration speed, as exhibited experimentally through isothermal calorimetry [9]. We will assume an Arrhenius-style law with dependence on both the temperature and pressure:

$$ \frac{d\xi}{dt} = \frac{A}{\eta_{\xi}} \exp \left( - \frac{E_A + \Delta V p_l}{RT} \right), \qquad (64) $$

with $E_A$ the apparent activation energy at atmospheric pressure, $\Delta V$ an activation volume, and $\eta_\xi$ the viscosity due to the diffusion of water through the layer of already formed hydrates. $\eta_\xi$ must then be increasing with $\xi$, as the layer of already formed hydrate will get thicker as hydration advances, and we will adopt the following form [29]:

$$ \eta_{\xi}(\xi) = \eta_{\xi,0} \exp \left( \bar{\eta} \frac{\xi}{\xi_{\infty}} \right), \qquad (65) $$

where $\eta_{\xi,0}$ and $\bar{\eta}$ are material parameters.

### 2.6. Evolution of volume fractions and poroelastic parameters

The poroelastic parameters evolve with hydration as the volume fraction of the solid increases. They can be measured experimentally, e.g. [9], or estimated using micromechanical models as in [11,13]. We adopt here a two-step homogenisation approach similar to that of Agofack et al. [1]: a first homogenisation at the CSH level followed by a second one at the paste level. Both upscalings are done using a self-consistent homogenisation scheme [34].

For the first homogenisation step, a classical $C_{1,7}SH_4$ composition is considered and the porosity of CSH is calculated taking into account a gel composition of $C_{1,7}SH_{1,8}$ as shown by Allen et al. [35]. This is rendered in the model through a coefficient $v_w^{CSH}$ which describes the number of moles of gel water per mole of $C_{1,7}SH_{1,8}$ gel, equal to 2.2 to retain the classical CSH composition. This leads to a constant CSH porosity of 35%. This assumption is further discussed in Section 5.2.

The second homogenisation step is a function of the hydration degree $\alpha = \xi/\xi_{\infty}$. The considered phases are water, clinker components ($C_3S$, $C_2S$, $C_3A$ and $C_4AF$) and hydration products (CSH, Portlandite, and other products classified as aluminates). Chemical shrinkage represents the remaining volume compared to the initial one. Table 1 presents the phase parameters necessary to determine the volumetric assembly. If $v^l$ is the volume fraction of species $i$, then its evolution is as follows:

$$ \begin{aligned} v^l &= v_0^l - V_m^l n^l \xi, \\ v^{cap} &= v_0^{cap} + V_m^H \sum_i (v_w^j - v_w^{CSH} v_{CSH}^j) n^l \xi, \\ v^{CSH} &= (V_m^{C_{1,7}SH_{1,8}} + v_w^{CSH} V_m^H) \sum_i v_{CSH}^j n^l \xi, \end{aligned} \qquad (66) \\ v^{CH} = V_m^{CH} \sum_i v_{CH}^i n^l \xi, \quad \text{for } i \in \{C_3S, C_2S, C_3A, C_4AF\}. $$

The volume fraction of chemical shrinkage can be calculated from the chemical shrinkage of each clinker phase:

$$ v^{shr} = \sum_{i} \Delta \overline{V}^{i} n^{i} \xi, \quad i \in \{C_{3}S, C_{2}S, C_{3}A, C_{4}AF\}. \qquad (67) $$

The volume fraction of aluminates is then finally given as the remaining volume fraction:

$$ v^{Al} = 1 - v^{shr} - \sum_{i} v^{i}, \quad i \in \{C_{3}S, C_{2}S, C_{3}A, C_{4}AF, H, CSH, CH\}. \qquad (68) $$

### 2.7. Field equations and transport properties

In the absence of body forces other than self-weight induced by gravity's acceleration $\gamma$, the mechanical equilibrium reads:

$$ \underline{\underline{\nabla}} \cdot \underline{\underline{\sigma}} = \rho \gamma. \qquad (69) $$

The field equation for the fluid is obtained by writing the conservation of the fluid mass, and hence:

$$ \frac{d m_w}{d t} = - \underline{\underline{\nabla}} \cdot \underline{\underline{M}}_w + v_w \underline{\underline{M}}_w \frac{d \xi}{d t}, \qquad (70) $$

where

$$ \underline{\underline{M}}_w = - \rho_w (p_l, T) \frac{k_l(S_l, e)}{\eta_w} \underline{\underline{\nabla}} p_l, \qquad (71) $$

with $\eta_w$ the dynamic viscosity of water ($\approx 1$ mPa·s at 20°C). The liquid permeability $k_l$ is not considered constant, but is a function of the saturation degree $S_l$ and the void ratio $e = \phi/(1 - \phi)$. Permeability decreases with decreasing saturation degree, and here we use the widely adopted model from Mualem [40]:

$$ k_r(S_l) = \frac{k_l(S_l)}{k_l(S_l = 1)} = \sqrt{S_l} \left( 1 - \left( 1 - S_l^{\frac{n}{n-1}} \right)^{\frac{n-1}{n}} \right)^2, \qquad (72) $$
<caption>**Table 1**</caption>

Physico-chemical parameters of the clinker phases used in the model. Values are calculated from Tennis and Jennings [36], Bentz [37], and Barnes and Bensted [38].

<table><caption>**Table 2**</caption><thead><tr><th rowspan="2">Phase parameters</th><th colspan="4">Clinker Phases</th></tr><tr><th>C<sub>3S</sub></th><th>C<sub>2S</sub></th><th>C<sub>3A</sub></th><th>C<sub>4AF</sub></th></tr></thead><tbody><tr><td>Molar volume [cm<sup>3</sup>·mol<sup>-1</sup>]</td><td>72.48</td><td>52.56</td><td>89.17</td><td>130.28</td></tr><tr><td>Water stoichiometry v'<sub>w</sub></td><td>5.3</td><td>4.3</td><td>10</td><td>14</td></tr><tr><td>CSH stoichiometry v'<sub>CSH</sub></td><td>1</td><td>1</td><td>0</td><td>0</td></tr><tr><td>CH stoichiometry v'<sub>CH</sub></td><td>1.3</td><td>0.3</td><td>0</td><td>-1</td></tr><tr><td>Chemical shrinkage ΔV'<sub>s</sub> [cm<sup>3</sup>·mol<sup>-1</sup>]</td><td>-11.74</td><td>-8.19</td><td>-30.72</td><td>-39.56</td></tr><tr><td>Solid volume balance ΔV'<sub>s</sub> [cm<sup>3</sup>·mol<sup>-1</sup>]</td><td>83.75</td><td>69.28</td><td>149.43</td><td>212.66</td></tr><tr><td>Latent heat of hydration L'<sub>t</sub> [kJ·mol<sup>-1</sup>]</td><td>114.2</td><td>431.0</td><td>237.8</td><td>204.1</td></tr></tbody><thead><tr><th rowspan="2">Phase parameters</th><th colspan="4">Water and Products</th></tr><tr><th>H</th><th>C<sub>1.7SH1.8</sub></th><th>CH</th><th>Al</th></tr></thead><tbody><tr><td>Molar volume [cm<sup>3</sup>·mol<sup>-1</sup>]</td><td>18.02</td><td>72.25</td><td>33.08</td><td>-</td></tr></tbody></table>

Mechanical properties of the different phases in cement paste (values from [39], and inverse analysis for C<sub>1.7SH1.8</sub> from pastes hydrated at 7 °C [14]).

<table><thead><tr><th></th><th>C<sub>3S</sub></th><th>C<sub>2S</sub></th><th>C<sub>3A</sub></th><th>C<sub>4AF</sub></th><th>C<sub>1.7SH1.8</sub></th><th>CH & Al</th></tr></thead><tbody><tr><td>K [GPa]</td><td>112.5</td><td>116.7</td><td>120.8</td><td>104.2</td><td>31.3</td><td>32.5</td></tr><tr><td>G [GPa]</td><td>51.9</td><td>53.8</td><td>55.8</td><td>48.1</td><td>21.4</td><td>14.6</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Relative liquid permeability as a function of the liquid saturation degree $S_l$ (n = 1.48).</figcaption></figure>

where exponent n refers to the same exponent as the one from the water retention curve (58). Corresponding evolution of the relative permeability as a function of the liquid saturation degree for n = 1.48 is presented in Fig. 6.

The decrease of the permeability with decreasing porosity – or equivalently, decreasing void ratio e – has also been the focus of several models, the most used being that of Kozeny and Carman [41], and that of Taylor [42]. Following the work of Picandet et al. [43] on cement pastes, we will adopt Taylor’s model, linking the logarithm of the permeability to the void ratio:

$$ \frac{\partial}{\partial e}(\log_{10} k_l) = \frac{1}{C_{k_l}} \qquad (73) $$

The thermal field equation is obtained by writing the heat conservation and takes into account the heat of hydration as an additional source term:

$$ T \left( \frac{\partial S}{\partial t} - \underline{\nabla} \cdot (s_w \underline{M}_w) \right) = - \underline{\nabla} \cdot q + L \frac{d\underline{\xi}}{dt}. \qquad (74) $$

where q is the heat flux and L the latent heat of hydration. Combining this equation with Fourier’s law:

$$ q = -k_T \underline{\nabla} T, \qquad (75) $$

gives the following form for the thermal field equation:

$$ T \left( \frac{\partial S}{\partial t} - \underline{\nabla} \cdot (s_w \underline{M}_w) \right) = k_T \Delta T + L \frac{d\underline{\xi}}{dt}. \qquad (76) $$

## 2.8. Continuity with the nonpercolated system

Before percolation, the solid is diffuse and a contiguous matrix is not yet formed. The fresh cementitious material is fluid, and as such the stress state can be considered as hydrostatic, in equilibrium with the fluid pressure:

$$ d\underline{\sigma} = -d\underline{p}_l \underline{1} \qquad (77) $$

It is interesting to note that this corresponds to the previous poromechanical Eq. (43) in the limit case where $C \to 0$ (which implies $b = 1$) and $S_l = 1$. Degenerating the remaining state equations yields:

$$ d\phi = \text{tr}(d\underline{\epsilon}) + \frac{1 - \phi_0}{K_s} d\underline{p}_l - 3\alpha_s(1 - \phi_0)dT \qquad (78) $$

since $b = \beta^p = \beta^v = 1$, and

$$ dS_s = -3\alpha_s(1 - \phi_0)d\underline{p}_l + \frac{(1 - \phi_0)C_{\sigma_s}}{T} dT. \qquad (79) $$

In absence of a real solid skeleton, $\phi$ will be understood as the volume fraction of fluid. This set of equations behaves as expected, namely in the simple case of undrained, isothermal and isobare early hydration (as, for instance, the standard ASTM C 1608 chemical shrinkage test). Indeed, recalling Eqs. (19) and (71), the water mass balance in a closed system imposes $dm_w = v_w M_w d\underline{\xi}$ such that:

$$ \begin{aligned} \frac{dm_w}{\rho_w} &= d\phi \\ &= d\phi_0 + d\phi \\ &= -\Delta\bar{V}_s d\underline{\xi} + \text{tr}(d\underline{\epsilon}), \qquad (80) \\ &= v_w \frac{M_w}{\rho_w} d\underline{\xi} \\ &= v_w \bar{V}_w d\underline{\xi} \end{aligned} $$

Hence

$$ \begin{aligned} \text{tr}(d\underline{\epsilon}) &= \left( \Delta\bar{V}_s + v_w \bar{V}_w \right) d\underline{\xi} \\ &= \Delta\bar{V} d\underline{\xi}. \end{aligned} \qquad (81) $$

The total volumetric deformation is therefore equal to the chemical shrinkage $\Delta\bar{V} d\underline{\xi}$, as expected in such a test.

As homogenisation of the poromechanical properties results in null stiffness and b = 1 under a certain threshold (corresponding to a critical capillary porosity of 0.5 if spherical inclusions are considered), the
<figure><img src="image_11.png" /><figcaption>Fig. 7. Phase assembly of the considered class G cement. (w/c = 0.44).</figcaption></figure>

<table><caption>**Table 3**
Composition of the class G cement.</caption><thead><tr><th>Cement composition</th><th>C<sub>3</sub>S</th><th>C<sub>2</sub>S</th><th>C<sub>3</sub>A</th><th>C<sub>4</sub>AF</th></tr></thead><tbody><tr><td>Solid mass fraction [%]</td><td>58.9</td><td>15.5</td><td>2.3</td><td>17.4</td></tr><tr><td>ξ'<sub>∞</sub> [mol·m<sup>-3</sup>] (w/c = 0.44)</td><td>3353.6</td><td>1168.8</td><td>110.7</td><td>465.5</td></tr><tr><td>Molar fraction n'</td><td>0.658</td><td>0.229</td><td>0.022</td><td>0.091</td></tr></tbody></table>

poromechanical equations can be used continuously to represent the behaviour of both unpercolated and percolated systems.

The nature of the deformation at the very early age (elastic, viscous or plastic) is irrelevant to the equations detailed above, and as such, they can be attributed to any of the mechanisms and still yield the previous results. However, since the deformations at that stage correspond to a rearrangement of solid grains in suspension following the creation of “voids” from the volume imbalance of the hydration reaction, it is clearly irreversible and will be considered hereafter as purely plastic.

## 3. Model parameter identification

### 3.1. Physico-chemical parameters

The physico-chemical parameters of the model, such as porosity variation, water consumption and heat release due to hydration can be estimated from values for the different clinker phases (C<sub>3</sub>S, C<sub>2</sub>S, C<sub>3</sub>A and C<sub>4</sub>AF) detailed in Table 1. Using the cement composition and mix proportions shown in Table 3, it is possible to aggregate all phases under a unique hydration extent with averaged properties. Let ξ = Σ<sub>i</sub> ξ<sup>i</sup>, i ∈ {C<sub>3</sub>S, C<sub>2</sub>S, C<sub>3</sub>A, C<sub>4</sub>AF} be this aggregated hydration extent, then

$$ \xi_{\infty} = \sum_{i} \xi_{\infty}^{i} = 5224 \text{ mol·m}^{-3}, \qquad (82) $$

and the averaged properties can be calculated as follows:

$$ v_w = \sum_{i} n' v_{w,i}, \quad \Delta \bar{V}_s = \sum_{i} n' \Delta \bar{V}_{s,i}, \quad L = \sum_{i} n' L_i. \qquad (83) $$

### 3.2. Phase assembly and poroelastic parameters

Fig. 7 shows the phase assembly calculated according to the relations given in Section 2.6, using values from Table 1 and the cement composition detailed in Table 3, taking into account a water-to-cement ratio w/c = 0.44. Note that this water-to-cement ratio, often considered optimal for full cement hydration, is found to be close to the minimum amount of water needed for full hydration in the model, with 5% volume fraction of remaining capillary water.

<figure><img src="image_12.png" /><figcaption>Fig. 8. Evolution of poroelastic parameters with hydration from homogenisation.</figcaption></figure>

Poroelastic parameters as functions of the hydration degree α = ξ/ξ<sub>∞</sub> are then calculated from the phase assembly and stiffness of the phases given in Table 2 using a micromechanical model and homogenisation method [1], as explained in Section 2.6. The resulting bulk and shear moduli as well as Biot's coefficient can be seen in Fig. 8. The evolution is close to linear after a threshold at α = 0.2 which can be identified as the percolation threshold. The determination of a percolation threshold from homogenisation means is strongly dependent on the aspect ratio considered for hydrates, as demonstrated by Sanahuja et al. [44]. However, in the absence of an experimental setting measurement such as a Vicat needle, the shape of the CSH was not taken into account in this homogenisation scheme. Hence, it likely constitutes an overestimation of the percolation threshold.

### 3.3. Hydration kinetics

Hydration kinetics is fitted on calorimetric data at a temperature of 7°C. The cumulative heat release from the calorimetry measurement can be directly related to the extent of hydration through the latent heat of hydration L. Fig. 9 shows the experimental data along with the model fit. The model predicts a hydration degree α ≈ 0.61 at 150 h. Relevant model parameters can be found in the recapitulative Table 4.

## 4. Simulation of undrained oedometric experiments

We propose here to simulate the experimental results of Agofack et al. [9]. These experiments consist of two main load cases: a constant
<table><caption>**Table 4**
Summary of model parameters.</caption><thead><tr><th>Symbol</th><th>Value</th><th>Unit</th><th>Description</th><th>Determination</th></tr></thead><tbody><tr><td colspan="5"><strong>Thermoporoelastic parameters</strong></td></tr><tr><td>E<sub>v</sub></td><td>3·10<sup>9</sup></td><td>Pa</td><td>Base viscous Young's modulus</td><td>Fitting to experiment (Section 4.1, 4.2)</td></tr><tr><td>ν<sub>v</sub></td><td>0.2</td><td>[-]</td><td>Viscous Poisson's ratio</td><td>Assumption</td></tr><tr><td>τ<sub>v</sub></td><td>1500</td><td>s</td><td>Viscoelastic relaxation time</td><td>Fitting to experiment (Section 4.1, 4.2)</td></tr><tr><td>α<sub>s</sub></td><td>1·10<sup>-5</sup></td><td>K<sup>-1</sup></td><td>Solid thermal expansion</td><td>From clinker thermal expansion [45]</td></tr><tr><td>C<sub>ox</sub></td><td>2.4·10<sup>6</sup></td><td>J·K<sup>-1</sup>·m<sup>-3</sup></td><td>Solid heat capacity at constant stress</td><td>For clinker with ρ<sub>s</sub> = 3200 kg·m<sup>-3</sup> [46]</td></tr><tr><td>L</td><td>1.04·10<sup>5</sup></td><td>J·mol<sup>-1</sup></td><td>Latent heat of hydration</td><td>Cement composition & Table 1</td></tr><tr><td>φ<sub>0</sub>(0)</td><td>0.586</td><td>[-]</td><td>Initial porosity</td><td>Phase assembly</td></tr><tr><td>ΔV&#773;<sub>s</sub></td><td>9.36·10<sup>-5</sup></td><td>m<sup>3</sup>·mol<sup>-1</sup></td><td>Solid molar volume variation</td><td>Cement composition & Table 1</td></tr><tr><td colspan="5"><strong>Hydration</strong></td></tr><tr><td>ζ<sub>∞</sub></td><td>5224</td><td>mol·m<sup>-3</sup></td><td>Ultimate hydration extent</td><td>Cement composition & Table 1</td></tr><tr><td>E<sub>Λ</sub></td><td>3.458·10<sup>4</sup></td><td>J·mol<sup>-1</sup></td><td>Hydration activation energy</td><td>Ref. [47]</td></tr><tr><td>ΔV</td><td>1.84·10<sup>-5</sup></td><td>m<sup>3</sup>·mol<sup>-1</sup></td><td>Activation volume</td><td>Determined from calorimetries in [48]</td></tr><tr><td>A<sub>0</sub>/k<sub>z</sub></td><td>0.0587</td><td>mol<sup>2</sup>·m<sup>-6</sup></td><td></td><td>Fitting on calorimetric curve</td></tr><tr><td>k<sub>z</sub>/η<sub>z,0</sub></td><td>1146</td><td>m<sup>3</sup>·mol<sup>-1</sup>·s<sup>-1</sup></td><td></td><td>Fitting on calorimetric curve</td></tr><tr><td>η&#773;</td><td>7.35</td><td>[-]</td><td>Water diffusion through hydrate coefficient</td><td>Fitting on calorimetric curve</td></tr><tr><td>ν<sub>iv</sub></td><td>5.89</td><td>[-]</td><td>Water consumption stoichiometric coefficient</td><td>Cement composition & Table 1</td></tr><tr><td colspan="5"><strong>Retention curve</strong></td></tr><tr><td>p<sub>r</sub></td><td>2·10<sup>4</sup></td><td>Pa</td><td>Reference pressure</td><td>Fitting on MIP data</td></tr><tr><td>n</td><td>1.5</td><td>[-]</td><td>Exponent of VG model</td><td>Fitting on MIP data</td></tr><tr><td>k<sub>id</sub></td><td>9.24</td><td>[-]</td><td>Viscosity due to water diffusion through hydrates</td><td>Fitting on MIP data</td></tr><tr><td colspan="5"><strong>Plasticity</strong></td></tr><tr><td>h<sub>z</sub></td><td>1·10<sup>8</sup></td><td>Pa</td><td>Chemical hardening modulus</td><td>Fitting to experiment (Section 4.1, 4.2)</td></tr><tr><td>h<sub>ε</sub></td><td>80</td><td>[-]</td><td>Isotropic hardening coefficient</td><td>Fitting to experiment (Section 4.1, 4.2)</td></tr><tr><td>M<sub>ε</sub></td><td>3.5</td><td>[-]</td><td>Chemical friction hardening coefficient</td><td>Fitting to experiment (Section 4.1, 4.2)</td></tr><tr><td>M<sub>0</sub></td><td>0.1</td><td>[-]</td><td>Initial friction ratio</td><td>Assumption</td></tr><tr><td>f<sub>i</sub></td><td>0.1</td><td>[-]</td><td>Isotropic tensile to compressive strength ratio</td><td>Assumption</td></tr><tr><td>τ<sub>p</sub></td><td>600</td><td>s</td><td>Viscoplastic relaxation time</td><td>Fitting to experiment (Section 4.1, 4.2)</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Cumulative heat release during hydration: calorimetric measurement and model fit.</figcaption></figure>

load sustained for up to 144 h while the cement paste is hydrating, followed by a cyclic loading and unloading between 5 and 65 MPa. In addition, some samples are first subjected to the loading cycle prior to the 144 h mark (namely at 15, 18, 20, 40, 60 and 103 h), and then loaded again at 144 h to show the effect of the prior loading cycles on the response at 144 h. We propose here to simulate all of these experiments. All simulations are performed through an implementation of the previously described model in MFront [49].

Results from the first two Sections 4.1 and 4.2 have been used in order to calibrate some of the parameters from the model (see Table 4).

<figure><img src="image_14.png" /><figcaption>**Fig. 10.** First undrained oedometric load case with monotonous loading with an axial stress P followed by load cycles between 5 and 65 MPa [9].</figcaption></figure>

Simulation results in both of these sections are therefore to be understood as results of a fitting procedure. However, in Section 4.3, results on the asynchronous loading do not result from specific fitting on these experiments and can therefore be seen as model predictions.

### 4.1. Shrinkage under constant pressure

The samples are left to hydrate for 144 h in an oedometric cell (no radial displacement) under constant axial load and temperature conditions (Fig. 10, “Monotonous” phase). We present here results at 7°C and axial loads of 3, 15, 25 and 45 MPa. Fig. 11 shows the axial strain ε<sub>zz</sub> as a function of time. Simulation results (dotted line) show
<figure><img src="image_15.png" /><figcaption>Fig. 11. Axial strain evolution in time during hydration of the cement paste under different axial loads.</figcaption></figure>

good agreement with the experimental data (solid line), up to an initial
stain offset of 0.3%, which is considered to correspond to the very early
shrinkage associated with the dissolution of the clinker phases. Indeed,
shrinkage during the “dormant” period is known to be due in part to
the dissolution process [50]. As this mechanism is not included in the
model, we look to replicate only the deformation occurring after the
end of the dormant period. The stabilisation of axial deformation after a
certain time for pastes hydrated at 3, 10 and 15 MPa is well reproduced.

It is found to be caused by the desaturation of the cement paste, as can be seen in Fig. 12, which shows the evolution of the pore pressure and radial stress during hydration. The pore pressure steadily decreases while the cement paste is hydrating, and eventually reaches desaturation. This was observed experimentally in setups that mimic borehole conditions by Li et al. [51] and Meng et al. [52]. The time at which desaturation occurs is also naturally observed to increase with the applied overburden pressure (OBP) in the simulated well. However, since the aforementioned experiments were performed at higher temperatures and with different cement pastes (class A OWC and class G with a water-to-cement ratio of 0.4, lower than in this study), this desaturation occurs at an earlier time (typically, in a range of 5 to 50 h). Radial stress also decreases with hydration with the same tendency as pore pressure, which could be expected. Hydration in oedometric conditions under axial pressures below 3 MPa could result in a loss of contact between the paste and the oedometer ring, as the radial stress reaches a zero value.

Fig. 13 shows the decomposition of the strain in elastic, viscous, and plastic parts during hydration. It is interesting to note the importance of both plastic and viscous strains, which represent a majority of the strain in the case of 25 MPa axial load. Indeed, the elastic moduli of the cement paste are quickly expected to be in the range of several GPa, so that elasticity could not explain axial strains of several percent under the considered loads.

4.2. *Loading cycle at* 144 h

Following hydration under constant axial load, a loading-unloading
cycle between 5 MPa and 65 MPa is performed on each sample (Fig.
10, "Cyclic" phase). Fig. 14 shows the stress-strain curves of this
load cycle after hydration under different axial loads, comparing the
experimental results of [9] and the simulation results of the present
model. Simulation results show good agreement with the experimental
data, namely in reproducing the trend of decreasing axial strain with

increasing load during hydration. The shape of the response is generally
well reproduced by the model, with the exception of the load cycle
after hydration under 3 MPa, which does not exhibit the change in slope
observed at around 2% axial strain. This slope change can be identified
in samples hydrated under 10 and 15 MPa as the point of resaturation
of the paste. The shape difference for 3 MPa can therefore be linked to
the discrepancy already apparent in Fig. 11 where the model under-
estimates the axial strain, thus overestimating the desaturation of the
sample. This leads to a higher strain needed for resaturation, which
is consequently not observed in the simulation. Using some different
refitted parameters for this experiment specifically can greatly improve
the fit.

4.3. Asynchronous loading

Other experiments performed by Agofack et al. [9] include loading
the paste with the same previous load cycle at one point in time prior to
the 144 h mark and investigate its effect on the paste’s response to the
load cycle at 144 h. Fig. 15 illustrates these load cases. We replicated
this load case numerically, with a sustained axial load of 25 MPa and
load cycles performed at 20, 40, 60 and 103 h.

Fig. 16 shows the responses to the virgin loading cycles. Experi-
mental results are qualitatively well reproduced, with a mainly linear
response at 20 h (Fig. 16(a)), shifting towards a clear elasto-plastic
response at 103 h with approximately 0.2% irreversible deformations.
Simulations seem to overestimate the elasto-plastic slope at 40 and
60 h, which remains close, in experiments, to the response prior to
percolation. This suggests that plastic hardening of the skeleton is
negligible at this stage of hydration and that deformation is driven
purely by the compression of capillary water. More complex hardening
rules would thus be needed in order to achieve a better rendition of the
response at that time.

Fig. 17 displays the responses of previously-loaded samples to a load-cycle at 144 h, compared to a virgin loading at 144 h (labelled as “Reference”). Experimental results are again well reproduced, with no significant impact of a load-cycle at 20 h on the response at 144 h, small impact for a load-cycle at 40 h, and more significant impacts for load-cycles at 60 and 103 h. The latter exhibits an almost purely viscoelastic cycle. Indeed, prior load cycles have the effect increasing the pore pressure by plastic deformation of the pore space. This induces a decrease in effective stress, which is the driver of plasticity onset. Hence, less plastic strains are observed in the load-cycle at 144 h in the case of prior loading.

4.4. *Fitting parameters sensitivity*

To conclude this section, we present a short analysis of the sensitiv-
ity of the presented results with respect to the main fitting parameters,
namely the viscosity and plasticity-related parameters. To this end, the
effect of the viscous modulus *E*<sub>v</sub>, the viscous relaxation time τ<sub>v</sub> and
the isotropic hardening coefficient *h*<sub>e</sub> are presented by comparing the
results of simulations using the parameters listed in Table 4 and results
of simulations after changing one of the fitting parameters. For each
parameter of base value *p*, we consider two alternative parameters with
values *p*/2 and 2*p*.

Fig. 18 shows the influence of the viscous modulus *E*<sub>v</sub>. It has an
important influence on both the axial strain in the monotonous loading
phase, with a decrease in total strain when *E*<sub>v</sub> is increased. Later, due to
the previous effect on the monotonous loading phase, the strain level in
the load cycle at 144 h is also affected. Indeed, as lower deformations
occur in the monotonous phase, larger deformations occur during the
load cycle.

As seen in Fig. 19, the viscous relaxation time, has little influence
on the monotonous loading part, which is expected as its largest value
is still small compared to 144 h. However, it has a pronounced effect
<figure><img src="image_16.png" /><figcaption>Fig. 12. Evolution of pore pressure $p_l$ and radial stress $\sigma_{rr}$ during hydration of the cement paste under different axial loads.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 13. Decomposition of strain in elastic, viscous and plastic parts during hydration of the cement paste under different axial loads.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 15. Asynchronous loading in an oedometric cell as performed by Agofack et al. [9]. The red line shows the load history of the sample from the previous load case, while green lines show samples undergoing a virgin loading at different times before a synchronous reloading at 144 h.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 14. Load cycles at 144 h after hydration under oedometric conditions at different stress levels.</figcaption></figure>

on the load cycle response, with a lower relaxation time leading to a lower apparent unloading modulus and overall larger strains.

Finally, Fig. 20 shows the influence of the isotropic hardening coefficient $h_e$. A marginal influence can be seen in the unsaturated portion of the monotonous loading phase, that is after 120 h. The effect on the load cycle response is more obvious but is still limited compared to that of the previous viscous parameters, with a maximum strain staying within ±5% of the base response, compared with ±15% approximately for the other two parameters.

To conclude, it appears that the response of the material in the considered oedometric tests from Agofack et al. [9] is highly influenced by the viscous behaviour of the material, especially its viscous stiffness.

## 5. Current limitations and future outlooks

The aim of this work was to show how couplings between different physical phenomena amount to a high level of complexity capable of describing the behaviour of cement paste at the early age in intricate experimental settings. Each submodel describing the relationship between some of the variables have thus been kept simple in order to focus on
<figure><img src="image_20.png" /><figcaption>Fig. 16. Responses to a first load cycle at different times during hydration under 25 MPa.</figcaption></figure>

their interplay rather than on the unit mechanisms themselves, while also limiting the number of parameters to determine. This naturally results in limitations of the model in its current form, when compared with other, more advanced models available in the literature catering to specific parts. Herein, we review several points of interest on which existing literature can provide avenues for improvement.

## 5.1. Creep behaviour

Creep of cement-based materials has been extensively researched for decades, both experimentally and numerically. The simple approach to viscoelasticity presented here, taken from the works of Carol and Bažant [26], does not fully describe the creep behaviour of cement-based materials, especially long-term creep, as the focus of the present paper lied in the early-age, specifically. More advanced creep models, such as those based on the micro-prestress theory developed by Bažant et al. [53], could be introduced to improve the modelling of the creep behaviour, as was previously done by Sercombe et al. [23] in a similar chemo-mechanical model for shotcrete.

## 5.2. Micromechanical upscaling

Evaluation of the evolution with hydration of the different poroelastic parameters was done here in a crude way, with several assumptions known to be limiting. First, the spherical shape of hydrates is a strong assumption. It is known to influence namely the solid percolation threshold, determined here to be close to $\alpha = 0.2$, which appears unexpectedly late. Considering acicular particular tends to decrease the hydration degree at which solid percolation occurs to more realistic values [44]. The traditional assumption of a constant CSH porosity (here equal to 35%) is another limitation. Densification of the CSH is known to occur as hydration advances, due to the its precipitation

in increasingly confined spaces in the gel porosity, as evidenced for instance by Königsberger et al. [54]. This also has an influence on the evolution of the poromechanical properties of cement-based materials with hydration.

A comprehensive model for determining the poromechanical properties of cement pastes as a function of hydration degree and water-to-cement ratio has been developed by Königsberger et al. [55], taking into account the transition from disc-shaped hydrates at low hydration degrees towards more spherical hydrates at higher hydration degrees, as well as the progressive densification of the CSH. The values for the Young's Modulus, Poisson's ratio, Biot's coefficient and modulus given in [55] could be used directly as inputs for this model, given a certain water-to-cement ratio.

## 5.3. Temperature dependency of material behaviour

Due to the geothermal gradient, along a well, the cement paste is hydrated under temperatures increasing with depth. It is well-known that the hydration temperature has an important influence on the chemical composition, microstructure and mechanical properties of the cement paste. For hydration temperatures below 100 °C, higher temperature results in lower C/S and H/S ratios in CSH composition, higher CSH density, higher capillary porosity and lower mechanical properties [56]. For temperatures exceeding 100°C, Portland cement hydration is affected, due to a change in proportions of CSH in favour of more crystalline phases such as α-C₂SH and hillebrandite, leading to strength retrogression phenomenon [57]. These temperature dependencies of the cement chemistry and their influence on the mechanical properties are not taken into account here. However, they could be integrated by introducing some temperature dependency of parameters related to microstructure and cement chemistry, and in what concerns
<figure><img src="image_21.png" /><figcaption>Fig. 17. Responses at 144 h following previous loading compared with no prior loading (Reference).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 18. Influence of the viscous modulus $E_v$ on the simulated results for an axial stress of 25 MPa and a load cycle at 144 h.</figcaption></figure>

the strength retrogression, through an additional chemical process represented by another $\xi$ variable altering the poromechanical properties of the material.

The temperature dependency of the pore structure, and notably the pore size distribution, influences also the water-retention curve of the material and its evolution with time. Whether there exists a direct influence of temperature on the water retention curve, or if this apparent temperature dependence is solely due to the different nature of the chemistry of the material at higher temperatures is also not yet determined. Further research down this alley could help confirm or infirm the hypothesis made in Section 2.

## 6. Conclusion and perspectives

The behaviour of cement-based materials, from their very early age up to their hardened state, depends on many couplings between mechanical, hydraulic, thermal, and chemical processes. In addition, this behaviour is inseparable from the progress of the hydration reaction of cement. We presented a theoretical model for the thermoporomechanical behaviour of cement-based materials from the early-age to the hardened state, along with its formal derivation from thermodynamic principles, with a maximum of generality, making it of use for a broad range of cement-based materials and possibly for other types of chemically hardening materials such as geopolymers. Taking into account many aspects of fresh cement-based materials, such as early strength
<figure><img src="image_23.png" /><figcaption>Fig. 19. Influence of the viscous relaxation time $\tau_v$ on the simulated results for an axial stress of 25 MPa and a load cycle at 144 h.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 20. Influence of the isotropic hardening coefficient $h_c$ on the simulated results for an axial stress of 25 MPa and a load cycle at 144 h.</figcaption></figure>

growth, heat release, water content and transport evolutions, this also makes this model a good candidate for simulations of 3D printed structures. Specifically, reported effects of air curing on hardened properties may be explained through the relative humidity dependence of the hydration kinetics.

This model aims at reducing the number of fitting parameters in favour of material properties that can be identified from theoretical or elementary values, or calculated by knowledge of the chemical composition of the cement. In particular, the phase assemblage at a given hydration extent, calculated from the cement composition, allows for the estimation of poroelastic properties through homogenisation.

Using the experiments solely to fit the viscous and plastic parameters, the model is able to reproduce experiments on ageing cement paste under oedometric conditions of Agofack et al. [9] over a wide range of pressures. The model adequately simulates the early-age shrinkage of cement paste under high pressures and finds that it is a combination of elastic, viscous, and plastic deformations. In particular, experiments on the effect of mechanical loading at the early age on the behaviour of a paste at a subsequent age, paramount in oil-well cementing as it can result in the loss of integrity of the cement sheath, are well captured by the model.

Further experimental investigations, such as triaxial experiments on fresh cement paste with diverse stress paths, could allow to refine the yield criterion and hardening laws and improve the generality of the model.

### CRediT authorship contribution statement

**Maxime Pierre:** Writing – original draft, Validation, Methodology, Investigation, Conceptualization. **Marcos Samudio:** Writing – review & editing, Methodology, Formal analysis, Conceptualization. **Siavash Ghabezloo:** Writing – review & editing, Validation, Supervision, Project administration, Methodology, Conceptualization. **Patrick Dangla:** Writing – review & editing, Validation, Supervision, Methodology.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### References

[1] N. Agofack, S. Ghabezloo, J. Sulem, Chemo-poro-elastoplastic modelling of an oilwell cement paste: Macroscopic shrinkage and stress-strain behaviour, Cem. Concr. Res. 132 (2020) 106046.

[2] P. Parcevaux, P. Sault, Cement shrinkage and elasticity: a new approach for a good zonal isolation, in: SPE Annual Technical Conference and Exhibition, SPE, 1984, p. 13176.
[3] H. Justnes, D. Van Loo, B. Reyniers, P. Skalle, J. Sveen, E. Sellevold, Chemical shrinkage of oil well cement slurries, Adv. Cem. Res. 7 (26) (1995) 85-90.

[4] K. Backe, P. Skalle, O. Lile, S. Lyomov, H. Justnes, J. Sveen, Shrinkage of oil well cement slurries, J. Can. Pet. Technol. 37 (09) (1998).

[5] S. Ghabezloo, J. Sulem, S. Guédon, F. Martineau, J. Saint-Marc, Poromechanical behaviour of hardened cement paste under isotropic loading, Cem. Concr. Res. 38 (12) (2008) 1424-1437.

[6] S. Ghabezloo, J. Sulem, J. Saint-Marc, The effect of undrained heating on a fluid-saturated hardened cement paste, Cem. Concr. Res. 39 (1) (2009) 54-64.

[7] M.-H. Vu, J. Sulem, S. Ghabezloo, J.-B. Laudet, A. Garnier, S. Guédon, Time-dependent behaviour of hardened cement paste under isotropic loading, Cem. Concr. Res. 42 (6) (2012) 789-797.

[8] V.N. Lima, F. de Andrade Silva, H.J. Skadsem, K. Beltrán-Jiménez, J.K. Sunde, Effects of confinement pressure on the mechanical behaviour of an oil well cement paste, J. Pet. Sci. Eng. 208 (2022) 109769.

[9] N. Agofack, S. Ghabezloo, J. Sulem, A. Garnier, C. Urbanczyk, Experimental investigation of the early-age mechanical behaviour of oil-well cement paste, Cem. Concr. Res. 117 (2019) 91-102.

[10] N. Agofack, S. Ghabezloo, J. Sulem, A. Garnier, C. Urbanczyk, Onset of Creation of Residual Grain during the hydration of oil-well cement paste, Cem. Concr. Res. 116 (2019) 27-37.

[11] S. Ghabezloo, Association of macroscopic laboratory testing and micromechanics modelling for the evaluation of the poroelastic parameters of a hardened cement paste, Cem. Concr. Res. 40 (8) (2010) 1197-1210.

[12] S. Ghabezloo, Effect of the variations of clinker composition on the poroelastic properties of hardened class G cement paste, Cem. Concr. Res. 41 (8) (2011) 920-922.

[13] S. Ghabezloo, Micromechanics analysis of thermal expansion and thermal pressurization of a hardened cement paste, Cem. Concr. Res. 41 (5) (2011) 520-532.

[14] M. Samudio, Modelling of an Oil Well Cement Paste from Early Age to Hardened State: Hydration Kinetics and Poromechanical Behaviour (Ph.D. thesis), Université Paris-Est, 2017.

[15] T. Sasaki, K. Soga, M. Abuhaikal, Water absorption and shrinkage behaviour of early-age cement in wellbore annulus, J. Pet. Sci. Eng. 169 (2018) 205-219.

[16] M. Meng, L. Frash, J.W. Carey, Z. Niu, W. Zhang, N. Guy, Z. Lei, W. Li, N. Welch, Predicting cement-sheath integrity with consideration of initial state of stress and thermoporoelastic effects, SPE J. 26 (06) (2021) 3505-3528.

[17] T.A. Petersen, F.-J. Ulm, Radial fracture in a three-phase composite: Application to wellbore cement liners at early ages, Eng. Fract. Mech. 154 (2016) 272-287.

[18] T. Petersen, F.-J. Ulm, Early-age stress and pressure developments in a wellbore cement liner: Application to eccentric geometries, J. Appl. Mech. 83 (9) (2016) 091012.

[19] O. Coussy, Poromechanics, John Wiley & Sons, 2004.

[20] O. Coussy, Mechanics and Physics of Porous Solids, John Wiley & Sons, 2011.

[21] F.-J. Ulm, O. Coussy, Modeling of thermochemomechanical couplings of concrete at early ages, J. Eng. Mech. 121 (7) (1995) 785-794.

[22] F.-J. Ulm, O. Coussy, Strength growth as chemo-plastic hardening in early age concrete, J. Eng. Mech. 122 (12) (1996) 1123-1132.

[23] J. Sercombe, C. Hellmich, F.-J. Ulm, H. Mang, Modeling of early-age creep of shotcrete. I: model and model parameters, J. Eng. Mech. 126 (3) (2000) 284-291.

[24] C. Hellmich, H.A. Mang, F.-J. Ulm, Hybrid method for quantification of stress states in shotcrete tunnel shells: combination of 3D in situ displacement measurements and thermochemoplastic material law, Comput. Struct. 79 (22-25) (2001) 2103-2115.

[25] L. Ma, Q. Zhang, Z. Jia, C. Liu, Z. Deng, Y. Zhang, Effect of drying environment on mechanical properties, internal RH and pore structure of 3D printed concrete, Constr. Build. Mater. 315 (2022) 125731.

[26] I. Carol, Z.P. Bažant, Viscoelasticity with aging caused by solidification of nonaging constituent, J. Eng. Mech. 119 (11) (1993) 2252-2269.

[27] L.K. Mettler, F.K. Wittel, R.J. Flatt, H.J. Herrmann, Evolution of strength and failure of SCC during early hydration, Cem. Concr. Res. 89 (2016) 288-296.

[28] M.T. Van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Am. J. 44 (5) (1980) 892-898.

[29] M. Cervera, J. Oliver, T. Prato, Thermo-Chemo-Mechanical Model for Concrete. I: Hydration and Aging, J. Eng. Mech. 125 (1999) 1018-1027.

[30] B. Albert, B. Guy, D. Damidot, Water chemical potential: A key parameter to determine the thermodynamic stability of some hydrated cement phases in concrete? Cem. Concr. Res. 36 (5) (2006) 783-790.

[31] R.J. Flatt, G.W. Scherer, J.W. Bullard, Why alite stops hydrating below 80% relative humidity, Cem. Concr. Res. 41 (9) (2011) 987-992.

[32] M. Wyrzykowski, P. Lura, Effect of relative humidity decrease due to self-desiccation on the hydration kinetics of cement, Cem. Concr. Res. 85 (2016) 75-81.

[33] W. Thomson, 4. On the equilibrium of vapour at a curved surface of liquid, Proc. R. Soc. Edinb. 7 (1872) 63-68.

[34] A. Hershey, The Elasticity of an Isotropic Aggregate of Anisotropic Cubic Crystals, American Society of Mechanical Engineers, 1954.

[35] A.J. Allen, J.J. Thomas, H.M. Jennings, Composition and density of nanoscale calcium-silicate-hydrate in cement, Nat. Mater. 6 (4) (2007) 311-316.

[36] P.D. Tennis, H.M. Jennings, A model for two types of calcium silicate hydrate in the microstructure of Portland cement pastes, Cem. Concr. Res. 30 (6) (2000) 855-863.

[37] D.P. Bentz, Three-dimensional computer simulation of Portland cement hydration and microstructure development, J. Am. Ceram. Soc. 80 (1) (1997) 3-21.

[38] P. Barnes, J. Bensted, Structure and Performance of Cements, CRC Press, 2002.

[39] K. Velez, S. Maximilien, D. Damidot, G. Fantozzi, F. Sorrentino, Determination by nanoindentation of elastic modulus and hardness of pure constituents of Portland cement clinker, Cem. Concr. Res. 31 (4) (2001) 555-561.

[40] Y. Mualem, A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res. 12 (3) (1976) 513-522.

[41] P.C. Carman, Fluid flow through a granular bed, Trans. Inst. Chem. Eng. (London) 15 (1937) 150-156.

[42] D.W. Taylor, Fundamentals of Soil Mechanics, vol. 66, (2) LWW, 1948.

[43] V. Picandet, D. Rangeard, A. Perrot, T. Lecompte, Permeability measurement of fresh cement paste, Cem. Concr. Res. 41 (3) (2011) 330-338.

[44] J. Sanahuja, L. Dormieux, G. Chanvillard, Modelling elasticity of a hydrating cement paste, Cem. Concr. Res. 37 (10) (2007) 1427-1439.

[45] H. Kada, M. Lachemi, N. Petrov, O. Bonneau, P.C. Aïtcin, Determination of the coefficient of thermal expansion of high performance concrete from initial setting, Mater. Struct. 35 (2002) 35-41.

[46] T. Honorio, B. Bary, F. Benboudjema, Thermal properties of cement-based materials: Multiscale estimations at early-age, Cem. Concr. Compos. 87 (2018) 205-219.

[47] L. D'aloia, G. Chanvillard, Determining the “apparent” activation energy of concrete: $E_v$ — numerical simulations of the heat of hydration of cement, Cem. Concr. Res. 32 (8) (2002) 1277-1289.

[48] X. Pang, W.C. Jimenez, B.J. Iverson, Hydration kinetics modeling of the effect of curing temperature and pressure on the heat evolution of oil well cement, Cem. Concr. Res. 54 (2013) 69-76.

[49] T. Helfer, B. Michel, J.-M. Proix, M. Salvo, J. Sercombe, M. Casella, Introducing the open-source mfront code generator: Application to mechanical behaviours and material knowledge management within the PLEIADES fuel element modelling platform, Comput. Math. Appl. 70 (5) (2015) 994-1023.

[50] F. Beltzung, F. Wittmann, Early chemical shrinkage due to dissolution and hydration of cement, Mater. Struct. 34 (2001) 279-283.

[51] Z. Li, J.M. Vandenbosche, A.T. Iannacchione, A. Vuotto, Characterization of oil well cement performance during early hydration under simulated borehole conditions, SPE J. 26 (06) (2021) 3488-3504.

[52] M. Meng, L. Frash, J.W. Carey, W. Li, N. Welch, Measurement of cement in-situ mechanical properties with consideration of poroelasticity, SPE J. 27 (05) (2022) 2655-2667.

[53] Z.P. Bažant, A.B. Hauggaard, S. Baweja, F.-J. Ulm, Microprestress-solidification theory for concrete creep. I: Aging and drying effects, J. Eng. Mech. 123 (11) (1997) 1188-1194.

[54] M. Königsberger, C. Hellmich, B. Pichler, Densification of CSH is mainly driven by available precipitation space, as quantified through an analytical cement hydration model based on NMR data, Cem. Concr. Res. 88 (2016) 170-183.

[55] M. Königsberger, B. Pichler, C. Hellmich, Multiscale poro-elasticity of densifying calcium-silicate hydrates in cement paste: An experimentally validated continuum micromechanics approach, Internat. J. Engrg. Sci. 147 (2020) 103196.

[56] S. Bahafid, S. Ghabezloo, M. Duc, P. Faure, J. Sulem, Effect of the hydration temperature on the microstructure of Class G cement: C-S-H composition and density, Cem. Concr. Res. 95 (2017) 270-281.

[57] K.J. Krakowiak, J.J. Thomas, S. Musso, S. James, A.-T. Akono, F.-J. Ulm, Nanochemo-mechanical signature of conventional oil-well cement systems: Effects of elevated temperature and curing time, Cem. Concr. Res. 67 (2015) 103-121.