Multiscale poromechanics of wet cement paste.

Capillary effects, such as imbibition drying cycles, impact the mechanics of granular systems over time. A multiscale poromechanics framework was applied to cement paste, which is the most common building material, experiencing broad humidity variations over the lifetime of infrastructure. First, the liquid density distribution at intermediate to high relative humidity is obtained using a lattice gas density functional method together with a realistic nanogranular model of cement hydrates. The calculated adsorption/desorption isotherms and pore size distributions are discussed and compare well with nitrogen and water experiments. The standard method for pore size distribution determination from desorption data is evaluated. Second, the integration of the Korteweg liquid stress field around each cement hydrate particle provided the capillary forces at the nanoscale. The cement mesoscale structure was relaxed under the action of the capillary forces. Local irreversible deformations of the cement nanograins assembly were identified due to liquid-solid interactions. The spatial correlations of the nonaffine displacements extend to a few tens of nanometers. Third, the Love-Weber method provided the homogenized liquid stress at the micrometer scale. The homogenization length coincided with the spatial correlation length of nonaffine displacements. Our results on the solid response to capillary stress field suggest that the micrometer-scale texture is not affected by mild drying, while nanoscale irreversible deformations still occur. These results pave the way for understanding capillary phenomena-induced stresses in heterogeneous porous media ranging from construction materials to hydrogels and living systems.

Over the lifetime of cement paste, the degree of water saturation can span a wide range: the initial reaction and precipitation of cement hydrates happen when cement powder is mixed with water to form the cement paste. Later, during setting, remaining water (not used by the dissolution–precipitation reaction) gradually evaporates from the pore space (1) that is composed of gel pores (1–5 nm) and capillary pores (5–50 nm). The set paste in the construction environment is then exposed to weather conditions that correspond to different relative humidity (RH) values. These changes in saturation level of cement paste induce capillary stress and material deformation as demonstrated in drying shrinkage experiments (2–4), which potentially contribute to the degradation and failure of buildings and bridges.

Structural changes due to liquid intake/drainage in porous media are known as “sorption-induced deformation” (5). Irreversible deformations during drying shrinkage of cement paste have been observed. Under mild drying conditions, small angle neutron scattering experiments (3) show plastic rearrangements at very small length scales, but no significant structural changes are detected at larger scales. However, harsh drying conditions can lead to large strain irreversibility at the macroscale (2, 4). It is also known that pressurization of the pore fluid in liquid-saturated rocks is an important weakening mechanism that can lead to fracture (6).

Capillary condensation and induced mechanical strain in materials consisting of simple pore structures [such as Mobil Composition of Matter No. 41 (MCM)-type mesoporous silica and carbon nanotubes] have been successfully described using independent pores of cylindrical geometry (5, 7–9). However, the multiscale pore structure of a heterogeneous material, such as cement paste, invalidates these theories mainly for two reasons. (i) The pores are connected to form a complex percolating topology, which gives rise to more complicated hysteretic behaviors in the sorption isotherms and liquid distributions. (ii) The liquid distribution inside the pores renders a heterogeneous force boundary condition, which along with a highly heterogeneous solid structure, challenges the applicability of classical continuum tools (10). Thus, to study the mechanics of undersaturated cement paste, a multiscale approach must be undertaken.

Homogenization methods play a central role in bridging different scales to understand how large-scale properties emerge from small-scale interactions. When the material is heterogeneous and partially saturated, it is necessary to reexamine the emergence of a continuum description. One important issue is how to determine a proper physical homogenization scale (11). For some cases of dense granular flows, the continuum postulate has been shown to be valid (11, 12). However, this has not been tested for poromechanics in the presence of capillary forces.

In this paper, we first simulate the liquid distribution in a heterogeneous mesoscale model of cement paste at partial saturations using lattice gas density functional theory (DFT). Fig. 1 shows the mesoscale model of cement paste composed of cement hydrates nanograins (5 nm) interacting with effective potentials from atomistic simulations (13–15). Our simulated hysteric adsorption/desorption isotherms of water and nitrogen compare well with experiments. Our subsequent calculations reveal the nanoscale details of capillary forces by constructing the Korteweg stress tensor field and integrating over Voronoi cells centered on cement hydrate nanograins (16). Using the capillary forces, we simulate structural relaxation and a subsequent sorption cycle. Analyses of the statistics of solid and liquid stress distributions point to the same homogenization length. We rationalize it as the length scale of structural heterogeneity. This example examines the continuum postulate for unsaturated porous media. Our results provide insights into drying/wetting of cement paste.

Fig. 1.

Multiscale cement paste model. Mesoscale structure of a realistic cement paste model spanning a cubic simulation box of 0.4 m. Gray spheres represent nanograins of cement hydrates of polydisperse diameters nm. Upper Inset is a zoomed-in illustration of a few cement hydrates showing the gel pores (1–5 nm) among the close-packed nanograins. The molecular structure of cement hydrates is shown in Lower Inset—yellow and red sticks are silicon and oxygen atoms, respectively, in silica tetrahedra; the blue and white spheres are oxygen and hydrogen atoms, respectively, of water molecules; and the green and gray spheres are inter- and intralayer calcium ions, respectively. The effective pairwise interaction potentials of the nanograins have been upscaled from atomistic simulations (13, 50). The mesoscale cement paste model is based in out-of-equilibrium precipitation and aggregation of cement hydrates (15) (). It reproduces realistically the mesoscale structure and mechanics of cement paste, including the complex morphology of capillary (5–50 nm) and gel pores (33). Lower Inset is reproduced with permission from ref. 14.

Sorption Hysteresis and Pore Size Distributions

To obtain the liquid distributions in cement paste, we simulated the adsorption/desorption isotherms for two adsorbates nitrogen and water in mesoscale configurations of cement hydrates using lattice gas DFT (17, 18) as described in . The sorption isotherms calculated for both adsorbates are in good agreement with experimental data as shown in Fig. 2 for nitrogen (77 K) and water (300 K), respectively. The nitrogen adsorption/desorption isotherm displays only minimal hysteresis in agreement with experimental observations of ref. 19. The mesoscale configuration of cement hydrates used for the nitrogen and water lattice DFT simulation is called “hardened,” as it is obtained right after out-of-equilibrium precipitation of nanograins (15) and has internal eigenstresses (20) ( and ).

Fig. 2.

Isotherms and PSDs. Adsorption/desorption isotherms for (A) nitrogen (at 77 K) and (B) water (at 300 K) in cement paste. Symbols show the mass ratios of wet to dry sample (; blue squares and red circles for adsorption and desorption, respectively) from experiments in ref. 19 for A and ref. 21 for B. Black lines show simulated isotherms on cement model configurations with realistic 3D pore network of water to cement ratio of 0.52. (B) For water sorption, the degree of saturation (DOS) estimated from the simulations is shown on the right y axis. C and D show PSDs calculated by the MC method on the three-cement model configuration and from applying the Kelvin equation to the isotherms of nitrogen and water. The y axis value is normalized probability density, and the x axis value is pore radius.

The water adsorption/desorption isotherm features a significant hysteresis. We focus our analysis on pure capillary effects at , and the results compare well with experimental data in ref. 21. Under harsh drying conditions, other chemophysical processes may take place, and these are beyond our scope (22, 23). Our approach allows us to explain the water hysteresis loop from metastable ink bottleneck states on desorption due to the pore constrictions in cement paste. This contrasts with previous results based on equilibrium thermodynamics for both adsorption and desorption branches and pore connectivity assumption (23, 24).

To investigate the poromechanical effect of internal eigenstresses and water capillary forces, we simulated a second cycle of adsorption/desorption on (i) the “aged” [the hardened configuration relaxed to zero average eigenstress by means of molecular dynamics (MD) simulation in the NPT (constant particle number, pressure, and temperature) ensemble] and (ii) the “capillary aged” (the hardened configuration relaxed to zero average eigenstress by MD simulation in the NPT ensemble under the action of capillary forces at , where capillary forces are the largest) (see below). All curves for the adsorption branch of water are superposable (Fig. 2). The desorption branches differ between the hardened and the two relaxed configurations, showing a densification tendency altering pore constrictions due to structural relaxation. We note that the distribution of pore sizes between the hardened and relaxed configurations shows no significant differences (Fig. 2).

Pore size distribution (PSD) is an important piece of information in the characterization of porous media, which experimentally is often derived from adsorption/desorption isotherms (25–27) using the Kelvin equation and surface adsorption models (28–30). However, isolated and geometrically idealized pore shapes, such as cylinders or spheres, are usually assumed in these calculations. Another way to calculate PSD is using Monte Carlo (MC) sampling of the pore space of simulated samples (as in this work) or experimentally acquired 3D tomography images (31). To access the pore connectivity down to the nanometer scale, techniques, such as NMR or X-ray tomography (27, 32), are used. Here, we apply the MC algorithm to the mesoscale cement paste models that have a realistic PSD in agreement with experimental data (33).

The extracted PSDs are compared in Fig. 2 for nitrogen and water, respectively. The MC results capture a population of gel pores ( nm) and large pores that extends to 25 nm. The Kelvin equation results for nitrogen show similar features. However, for water, the Kelvin equation results feature bimodal distributions with both peaks smaller than 5 nm, similar to the curves derived from experimental data in ref. 21. The lack of large pores derived from the Kelvin equation using water sorption isotherm data is due to the nonlinear relationship between the Kelvin radius and the so that pores larger than 5 nm will correspond to > 95%. Our results indicate that the PSDs derived from nitrogen sorption isotherms are consistent with the real morphology of the pore space, while those derived from water sorption isotherms lack the large capillary pores. However, such incomplete PSDs from water sorption experiments are often reported in the literature.

Capillary Forces at Nanograins Level

Having the water distributions from the lattice gas DFT simulations, we proceed to calculate the capillary forces on the nanograins following our previous work described in ref. 16. At low water saturation, capillary bridges between two or more nanograins can be described by the Kelvin equation. This traditional method requires increasing effort for more nanograins (34, 35) and fails when liquid clusters starts percolating. The advantage of our method (16) is that it enables us to calculate capillary stress at arbitrary liquid saturation levels. We first construct the Korteweg stress field (36, 37) based on the liquid density and its gradient (). Then, we integrate the stresses over the surface of Voronoi cells of each nanograins to obtain the capillary forces.

Fig. 3 shows the water capillary force per nanograin distribution computed at different values. Qualitatively, capillary forces increase in magnitude when decreases. The forces at (the closing point of the hysteris loop) exhibit a long-tailed distribution at the particle level, with a mean value around 0.7 nN, and the largest is above 1.8 nN (Fig. 3). The maximum of (3 nN) is still smaller than the magnitude of solid–solid particle interaction forces ( 5–10 nN). Fig. 3 shows the magnitude of the capillary force on each particle of a thin section of the 3D hardened cement paste.

Fig. 3.

Distribution and visualization of capillary forces. (A) Distribution of capillary forces at the particle level. (B) Color code shows the magnitude of on a thin section of the sample at . The large forces are concentrated on the surface of pores as highlighted in C and D. The green segments in C and D compose the pore–solid interface. In C, only particles having nN are shown in the same color code as B, while D shows only force vectors in light blue. Most of these forces are pointing into the solid packing region, effectively densifying the gel texture. The simulation box size is 400 nm.

Fig. 3 depicts the particles experiencing the largest capillary forces at (magnitude and vector, respectively). The small gel pores are always filled until low , while large capillary pores remain empty even at high , thus creating the largest capillary forces between the gel pores (1–5 nm) and the capillary pores (5 nm). The overall capillary effect at is densification of the solid texture, since the capillary forces enhance attractions between nanograins surrounding small capillary pores.

Irreversible Rearrangements Due to Capillary Stresses

We further investigate the relaxation of the solid structure with eigenstresses (hardened sample) under wet conditions using MD simulations with capillary and pairwise particle interactions. We assume the capillary forces to be constant along the MD trajectories. This assumption is proven to be valid, as capillary forces do not induce large particle displacements as shown in Fig. 4, where the distributions of particle displacement magnitudes (, being the average particle diameter) are displayed. Under the action of capillary forces in all , only 0.1% of the particle population undergo displacement comparable with their own size (i.e., a few nanometers). This is due to the fact that the magnitude of the capillary forces is in general five times smaller than interparticle cohesion. As the value decreases, capillary forces increase, and the peak of distribution shifts to larger values. At larger displacements, the distributions at all exhibit an algebraic decay with characteristic exponent equal to −2.5.

Fig. 4.

Structural relaxation and nonaffine displacements due to capillary forces. (A) Distributions of particle displacement magnitudes ( is the average particles diameter) at different values. (B) Distributions of nanograins nonaffine displacements . The aged sample is the reference state. (C) Spatial correlation function of . Power law fits are shown as dashed lines. (D) PE after structural relaxation simulated with and without capillary forces (aged and capillary aged, respectively) at different values. The initial PE of the hardened paste before relaxation is also shown. For , capillary stress assists the relaxation process and lowers the PE of the system. PE values are in reduced units.

To investigate further the effect of capillarity in the structural relaxation of cement, we measure the drying-induced irreversible displacements in the capillary aged sample, taking as reference the aged sample. This quantity is akin to what in glass physics is coined nonaffine displacement, which in our case, is a result of capillary pressure (more details are in ).

Fig. 4 shows the distributions of nonaffine deformations that exhibit long-range correlations with algebraic decay with exponent of −1 at all . However, additional analysis of by means of the spatial correlation functionshows an algebraic behavior of exponent restricted to short distances (R 35 nm) (Fig. 4) by contrast to observations for dense colloidal glasses under shear deformations (38, 39). We tentatively attribute this behavior to the presence of capillary pores that localized these nonaffine deformations at their vicinity.

We note that high RH values () result in a lower potential energy (PE) of the system than the fully wet aging condition (no capillary effects) after relaxation as seen in Fig. 4. This range of values corresponds to capillary forces of smaller amplitude (Fig. 2). Capillary forces of larger amplitude at RH increase the PE.

Stress Homogenization

The capillary forces calculated on the cement hydrate nanograins show a heterogeneous distribution (Fig. 3). We adopt the Love–Weber homogenization (40, 41) scheme to coarse grain the nanometric capillary stress distribution (). For solid and capillary stresses, we tested different homogenization scales (Fig. 5) ranging from 10 to 100 nm in all of the cement paste models considered in this work. The distribution of solid stress in the hardened sample shows a maximum at nonzero pressure due to the eigenstress acquired during setting. The distributions of solid stresses in the aged and capillary aged samples are similar in amplitude and width (). The shape and width of the capillary stress distribution at for the hardened sample change with (skewness in Fig. 5 and ). However, Fig. 5, which presents the first four moments of these distributions, indicates that a characteristic length scale emerges around 35 nm. Above this scale, the first moments converge to the macroscopic mean, the second moments monotonically decreases, and the excess kurtosis is negative or close to zero, indicating shorter tails than the normal distribution.

Fig. 5.

Homogenization length scales for solid and liquid stresses. A–D show the statistics of solid volumetric stress () in hardened, aged, and capillary aged samples. E–H show the statistics of liquid pressure . Black lines denote the homogenization length = 35 nm, common for all samples/stresses. The simulation boxes are 400 nm3.

Other than the heterogeneity of the structure, it is interesting to note that is an order of magnitude smaller than the simulation box size ( nm), and it is also similar to the correlation length of in Fig. 4. The emerging homogenization length common to both capillary stress and solid stress rationalized as the scale of long-range spatial correlations inherent in the heterogeneous texture has been shown in ref. 33 via scattering data and chord length distributions.

Averaging over the entire simulation box results in the overall capillary pressure at the micrometer level as shown in Fig. 6. Cement paste, although a heterogeneous porous medium, displays a hysteretic capillary pressure similar to that found for simple geometry porous media (5, 9). Fig. 6 shows the homogenized capillary pressure calculated at = 35 nm at different values during adsorption (Fig. 6, Upper) and desorption (Fig. 6, Lower). In absolute value, the homogenized capillary stress is larger on desorption, and it is more intense in the densest matrix regions.

Fig. 6.

Homogenized stresses. (A) Hysteresis of average capillary pressure (in MPa). Negative sign indicates shrinkage. (B) Hysteresis of homogenized capillary pressure displayed at the mesoscale. Upper shows adsorption, and Lower shows desorption. The simulation boxes are 400 nm3. Homogenization scale 35 nm.

Conclusion

In this work, we investigated the role of capillary forces in the relaxation mechanism of a granular model of cement paste under partial saturation from the nanograin level (5 nm) to the mesoscale (400 nm). We achieved agreement with nitrogen and water sorption experiments and assessed the validity of the Kelvin equation in determining the PSD of cement paste. We showed that nitrogen desorption data give better PSD than water compared with the PSD obtained from 3D stereological analysis. We also considered one more drying/wetting cycle showing structural densification but without significant change of the PSDs.

The framework that we developed based on the liquid distribution and stress field at the nanograin level enables us to show that the largest capillary forces are concentrated at the interface between gel pores and larger capillary pores. We found that the largest capillary forces are obtained at the closure point of the hysteresis loop (RH = 31%) and correspond to menisci located at the interface between gel and capillary pores. These forces induce local irreversible rearrangements that have relatively short spatial length correlations due to the presence of the capillary pores spanning the samples and hence, do not exhibit the long-range character found in colloidal glasses during inhomogeneous shear (38, 39). Under additional mechanical loading and longer timescales or with multiple drying/wetting cycles, we speculate that they may contribute altogether to creep and fracture, similar to observations in amorphous materials (42, 43) and fluid saturated rocks (6).

The results suggest a homogenization length scale common to both liquid and solid stress fields correlated with the scale of structural heterogeneities. Therefore, this length sets the representative volume element size for upscaling poromechanics using either discrete or continuum descriptions, such as the lattice element method (LEM) method (44).

In soft porous materials, capillary forces are comparable with the solid cohesive forces and may induce large deformations, such as drying of colloidal suspensions (45, 46). The framework that we proposed can be readily adopted in these materials.

Materials and Methods

Sorption Simulations and Stress Calculation.

The sorption simulations are based on the lattice gas DFT model (17, 18) by minimizing the grand potential: , where denotes the normalized density of fluid on site i continuously varying from 0 to 1; or 1, indicating that site i is occupied by solid or vacant, respectively; and and are the fluid–fluid interaction and fluid–solid interaction energy, respectively. is determined by the bulk critical point , with the number of nearest neighbors. is imported from atomistic simulation data as the isosteric heat of adsorption (i.e., the differential heat of adsorption) in the limit of zero coverage. For water and nitrogen adsorption in cement, (22). To estimate the lattice spacing , we estimate the surface tension, which is energy per area, . For nitrogen at K, 8.94 mN/m, which gives 0.345 nm. For water at K, 72 mN/m, which gives 0.24 nm. Based on the estimates, we choose a fine-grid cell size of Å. Taking g/mL for solid-phase density of cement hydrates, we convert degree of saturation measured from the simulations to total amount of absorbed liquid in cement paste. The Korteweg stress tensor can be derived as () , where is the identity tensor and is an arbitrary constant tensor. is the asymptotic bulk value of the hydrostatic pressure.

Simulations of Structural Relaxation.

The initial mesoscale configuration of cement hydrates is the hardened configuration that has an average eigenstress −50 MPa. Structural relaxation MD simulations of the hardened sample were carried out using LAMMPS (47) on NPT conditions of zero pressure and room temperature. Capillary forces at different values were constantly applied during the entire simulation. The hardened configuration under the action of capillary forces at , where the distribution of capillary forces has the longest tail, is called capillary aged. All simulations were terminated after 500,000 MD steps with time step in reduced unit when the system had converged. From the particle positions, nonaffine displacements were computed. The nonaffine displacement is defined as , where runs over neighbors of the central particle indexed by 0 (48, 49). corresponds to the deformation gradient in continuum mechanics, and essentially captures the higher-order deformations that are usually not considered in linear theories. The reference configuration for the nonaffine displacements was the aged configuration. The aged configuration is final state of the hardened configuration after relaxation to zero average eigenstress by MD simulation in the same NPT conditions (zero pressure and room temperature) but without capillary forces (only interparticle interactions).