Integrated analysis of pore size, shape, and position on mesoscale compressive strength of cementitious materials: Numerical modelling and statistical regression

The analysis used a mixture of clay, CEM I cement (Górażdże 42.5 R), and water. The tested soil was predominantly composed of silt (Si), accounting for 68.8% of the total particle distribution (particle size between 0.002 and 0.063 mm). Clay particles (Cl), defined as particles smaller than 0.002 mm, constituted 28.73% of the material, while sand (Sa) content, with particle sizes ranging from 0.063 to 2 mm, was minimal at only 2.47%. This composition indicates a fine-grained soil with a dominant silt fraction and a relatively high clay content. Each batch consisted of 1,000 g of water, 900 g of cement, and 1,350 g of clay. The resulting water-to-cement (w/c) ratio of 1.1 is typical for injection works in geotechnical engineering, such as jet grouting [24]. The mixing process was carried out by the PN-EN 196-1 standard [25]. The prepared material was placed into syringes with a diameter of 10 mm. After a curing period of 28 days, the syringes were cut open and sectioned to obtain test specimens with a height of 10 mm (±0.2 mm). The samples were examined using a GE Phoenix M300 computed tomography scanner and a DEBEN CT5000TEC loading stage. A high-resolution scan was performed with an accuracy of 10 µm, which enables detailed evaluation of the internal structure [26]. Tomographic data were used to assess porosity through segmentation algorithms [27], which are particularly effective in detecting air pores. Image analysis of the tomographic data was conducted using VGStudio 3.5 software. Porosity analysis was also performed in the same software. The detection was carried out according to the VGDefx. Image analysis parameters and filtering criteria were adjusted individually based on the absorption histograms. These data also served as the basis for reconstructing the actual geometry of the analysed sample [20]. A numerical model was generated in FLAC3D 9.0, consisting of elements with a size of 50 µm. The material behaviour was described using the ‘concrete’ model [28,29]. This concrete model is a plastic-damage formulation based on the concepts of fracture-energy-driven damage evolution and modulus degradation, as defined in continuum damage mechanics. The numerical model was calibrated using damage patterns observed during in situ compression tests inside the CT scanner. A detailed description of the calibration process is presented in a separate paper [30]. The calibration involved adjusting the stress–strain curve to closely match the behaviour observed in the physical sample. Subsequently, energy accumulation parameters were refined to ensure that the simulated damage pattern corresponded with the one recorded in the CT chamber. Digital Volume Correlation was also employed in the displacement analysis. Figure 1a–c shows the sample geometry, a cross-section of the failed sample, and the result of the calibrated numerical simulation. An isometric view of the analysed sample, along with the porosity analysis, is shown in Figure 1d. The pore colours indicate variations in their sizes. A single elongated pore has a diameter exceeding 5 mm. Subsequently, Figure 1e presents the distribution of the analysed pores according to their diameters. The vast majority of pores are within the range of d > 100 µm. All pores identified in the analysis have diameters greater than 50 µm.

Figure 1

Material used in the study: (a) test specimen; (b) vertical CT cross-section of the failed specimen; (c) result of the calibrated numerical analysis; (d) isometric view; (e) distribution of pores.

The calibrated results were validated against the failure mechanism of a second specimen from the same material batch. With the calibrated material parameters, the damage pattern observed in shear strain maps for the second numerical model showed strong agreement with the experimental observations. The material parameters used in the analysis are summarised in Table 2, and the stress–strain relationship ϭ(ε) from both the experiment and the simulation is presented below. In the numerical calculations, the model accounted for the following parameters: Young’s modulus (E), fracture energy in uniaxial compression/tension (G _c, G _t), uniaxial initial yield strength compression/tensile (f _c0, f _t0), or maximum uniaxial yield strength compression/tensile (f _cm, f _tm). Mesoscale parameters differ significantly from those obtained directly through experimental measurements. All of these parameters were calibrated as shown in [30].

The graph shown in Figure 2 presents the σ(ε) relationship for both the experimental results and the calibrated numerical analysis. The simulation results are depicted with a dotted line. In the range of small strains (up to 0.1%), discrepancies are observed, associated with the settlement of the actual samples during testing. A nonlinear stress increase occurs at the initial loading stage, suggesting the closure of initial pores and microcracks within the cemented soil structure. The material exhibits predominantly elastic behaviour during this phase, with gradually increasing stiffness.

Figure 2

Uniaxial compression. Stress–strain σ(ε) curves for the tested samples and the corresponding simulation results.

In the quasi-linear phase (0.01 < ε < 0.03), stress increases linearly, indicating the development of elastic deformations. The value of σ _max (approximately 14–15 MPa for sample 6 and 12–13 MPa for sample 5) corresponds to the maximum load-bearing capacity of the material, beyond which degradation processes commence. In this phase, the microstructure begins to experience intense cracking and localised failure. The post-critical phase, occurring for ε > 0.05, follows the peak stress. After reaching σ _max, the stress decreases, indicating damage and disintegration of the cemented soil matrix. The experimental data reveal a more abrupt stress drop, suggesting brittle fracture and loss of structural continuity. In contrast, the numerical simulation shows a more gradual stress decline, resulting from using a plastic-damage material model that does not fully capture the sudden failure mechanisms of the real material. Another factor influencing the difference between experimental and numerical behaviour is the necessity to halt the loading process for CT scanning. The strain increment was not continuous at a constant rate but was interrupted at four discrete points to perform tomographic scans (under constant load without unloading), which may have influenced the experimental response.

The experimental and numerical data exhibit a typical pattern for compression tests: an initial phase of stress increase, reaching a maximum, followed by a decline associated with material degradation. The simulation curves (dotted lines) are generally smoother than the experimental curves (solid lines), which may be attributed to numerical averaging within the model. The numerical model is based on a continuum approach. The simulations predict higher maximum stress values than experimentally, particularly for sample 6. In the case of sample 5 (thin dotted vs. thin solid line), the maximum values are closer, although the simulation still slightly overestimates the stress. After the peak stress, the real cemented soil exhibits a more abrupt degradation, suggesting that the numerical model may lack the complete representation of localised cracking phenomena.

Numerical model

The study involved performing numerical analyses on models matching the dimensions of the specimens used for model calibration. In FLAC3D, 20 models were generated – one solid reference model and nineteen models containing pores. The geometric parameters of these models are summarised in Table 3 The models are labelled with a number (pore size, shape) and a letter (position). Sample 0 served as the reference specimen, against which subsequent models could be compared. The analyses focused on the influence of pore size (with radius ranging from 1 to 2 mm) and sphericity (ranging from 0.79 to 1.00). Sphericity is a dimensionless parameter that quantifies how closely the shape of an object approaches that of a perfect sphere. The value of sphericity ranges from 0 to 1, where 1 represents a perfect sphere. Objects with lower sphericity are more elongated, flattened, or irregular. Sphericity is an essential parameter in material science, geotechnics, and concrete research because it affects packing density, mechanical behaviour, and the development of stress concentrations around inclusions or pores. Changes in sphericity were achieved by scaling the pore dimensions along two axes. Scaling one horizontal axis by a factor of 0.66 and the other by 1.5 (while keeping the third axis unchanged) reduced the sphericity to 0.91 while maintaining the initial pore volume. A similar method was used to obtain a sphericity of 0.79. The numerical models can be grouped into several categories. Simulations 1–3 include spherical pores of varying sizes. Models 4 and 5 represent the smallest pores with different sphericities. In pairs 6 and 7, 8 and 9, the pore volume and the degree of sphericity vary. Figure 3 shows the geometries used in the simulations. Figure 3a illustrates three different pore shapes, highlighting the impact of sphericity changes.

Figure 3

Numerical models: (a) pores with constant volume and varying sphericity, from left to right: 1.0, 0.91, 0.79 – isometric view; (b) three pores of identical diameter shown at positions A, B, and C; (c) model’s geometry view – pore no. 1, position C, FLAC3D.

The individual pores were positioned at specific locations within the sample. Position A corresponds to the placement where the pore’s geometrical centre is located at the centre of the specimen, which under ideal conditions corresponds to the zone of maximum shear stress concentration. The height of the pore centre in position A is 5 mm. Position B is located below the specimen centre, with the pore centre at a height of 4 mm. Position C places the pore centre at 3 mm. Figure 3b presents an example of pore positioning within the sample.

The numerical model was composed of roughly 10 million cubic elements, each with a side length of 50 µm. To simulate a rigid, non-deformable loading surface, the top layer was constrained with zero-velocity boundary conditions in all spatial directions. The bottom layer was subjected to a constant vertical velocity of 0.05 µm per cycle, while horizontal movements were fully restricted. Numerical stabilisation was performed after every 200 cycles, corresponding to a total displacement of 10 µm.

The effect of pore shape was studied for position C (pore centre height at 3 mm). This analysis location was selected because it is below the sample centre, where maximum shear stress concentrations occur under peak loading conditions. The analysis was carried out for three pore sizes and three sphericity values.

The graphs show differences in the simulation results for samples with the smallest pore diameter (Figure 4a) that are not initially visible. As the strain increases, differences between the reference sample and the samples containing pores become more apparent. All samples with pores exhibit lower stress levels compared to Sample 0. Within the group of pores with a radius of 1 mm, the differences caused by pore shape remain quite subtle – the results for sphericity 1.0 (Sample 1C) and 0.91 (Sample 4C) are very similar. In contrast, Sample 5 (sphericity 0.79) shows slightly lower stress values.

Figure 4

Stress–strain curves with varying sphericity (a) pore radius: 1 mm; (b) pore radius: 1.5 mm; and (c) pore radius: 2 mm.

Maximum stresses are reached in the 5–6.5% strain range before degradation effects dominate. The reference sample maintains the highest stress level, while the presence of pores reduces the load-carrying capacity. Comparing the samples with a 1 mm radius, it is observed that Sample 1C (perfectly spherical) generally reaches slightly higher stress levels than the samples with lower sphericity. Samples 4C and 5C exhibit very similar behaviour, although Sample 5 (sphericity, 0.79) occasionally achieves marginally lower stress values, which may indicate higher stress concentrations due to the irregular pore shape.

As shown in Figure 4b, all specimens – the reference and those with a 1.5-mm pore – show similar stress values at very low strains. As the strain increases, the reference specimen’s stress remains noticeably higher. For example, at a strain of 3%, Sample 0 reaches approximately 11.14 MPa, while the 1.5 mm specimens show lower stresses (around 10.99 MPa for Sample 2, with comparable values for Samples 6 C and 7 C). At very low strains, all specimens – the reference and those with a 1.5-mm pore – show similar stress values. At a strain of 5%, Sample 0 exhibits a stress of roughly 17.36, whereas the 1.5 mm specimens show markedly reduced stresses – approximately 15.37 MPa for Sample 2, 16.40 MPa for Sample 6C, and 16.58 MPa for Sample 7C. While one might initially expect the ideally spherical pore (Sample 2C, sphericity = 1) to produce the least stress reduction, the data indicate that for a pore of 1.5 mm, the sample with the lowest sphericity (Sample 7, sphericity = 0.79) maintains slightly higher stress values at higher strains compared to its more spherical counterpart. For the 1.5 mm pore size, the irregularity of the pore (lower sphericity) appears to modify the local stress distribution such that Sample 7C sustains higher stresses at elevated strains than Sample 2C. Nonetheless, all pore-containing specimens are still significantly below the performance of the pore-free reference.

The samples with a 2 mm pore located at position C are shown in Figure 4c. All specimens – the pore-free reference and those containing a 2 mm pore – exhibit similar stress values at low strains. Although the presence of the pore remains the dominant factor reducing mechanical performance, subtle effects related to pore shape are observable. The ideally spherical pore (Sample 3C) generally corresponds to a lower stress response at higher strain levels. Specimens with lower sphericity (Samples 8C and 9C) occasionally record slightly higher stresses, indicating a complex interaction between pore shape and stress distribution under large deformations.

Compared to the reference sample without pores, the specimens with a 1-mm pore radius exhibit lower stress values at a given strain level. This reduction results from the presence of the pore, which induces local stress concentrations. These findings demonstrate that, although the impact of changing sphericity at a constant 1 mm radius is not appreciable, there is a clear tendency for more irregular pore shapes to decrease the load-bearing capacity further. In practice, pore shape cannot be directly prescribed, particularly in injection-based cementing technologies; however, experimental studies on grouted granular media demonstrate that injection parameters and grout properties can significantly modify pore structure and effective porosity by filling existing voids and altering the pore size distribution. As a result, technological control primarily affects the void system indirectly through changes in porosity and microstructural homogenisation rather than through direct geometric control of individual pores [31]. While a pore of 1.5 mm in size compromises the material’s overall stiffness and strength compared to a pore-free specimen, the irregularity of the pore (as evidenced by a lower sphericity) may mitigate the reduction in stress-bearing capacity at higher strains. This behaviour contrasts with the observations for smaller pores (e.g., 1 mm), where a more spherical shape tends to preserve stiffness better under higher strain. These nuances highlight the importance of considering pore size and shape when evaluating cement grout’s mechanical performance. In some strain ranges, Sample 8C and Sample 9C exhibit slightly higher stress values than Sample 3C. This suggests that a deviation from perfect sphericity may alter the stress distribution for larger pores to partially compensate for the reduced load-bearing capacity. However, it should be emphasised that the variation in strength for the specimens with the largest pore did not exceed 4% of the maximum strength. It is important to note that the analysis was performed only for a single pore within a mesoscale specimen.

To investigate the effect of vertical pore placement on the mechanical behaviour of cemented soils, uniaxial compression simulations were performed on specimens with a single spherical pore of radius 1, 1.5, or 2 mm. For each configuration, the pore was located at three different heights: mid-height (5 mm from the base, denoted as A), slightly lower (4 mm, B), and near the bottom of the specimen (3 mm, C). Stress–strain responses were compared to identify the influence of pore position on global material strength. The stress–strain relationships are presented in Figure 5a.

Figure 5

Results for the effect of pore position: (a) stress–strain curves for different pore locations; (b) vertical stress distribution map for position A; (c) vertical stress distribution map for position B; (d) vertical stress distribution map for position C; (e) shear strain distribution map for position A; (f) shear strain distribution map for position B; and (g) shear strain distribution map for position C.

The results consistently indicate that specimens with pores closer to the sample’s base (position C) exhibit higher peak stresses and greater residual load-bearing capacity than those with centrally located defects (position A). This trend becomes more pronounced at larger strains, beyond the linear elastic level. The highest degradation in mechanical response was observed when the pore was placed near the cones of maximum shear stress (positions A and B).

This phenomenon can be linked to the well-established stress distribution in uniaxial compression tests, where the central region of the specimen is known to accumulate the highest shear stresses. The presence of a defect in this critical zone leads to earlier initiation of plastic deformation and microcracking, accelerating the overall material degradation. Conversely, pores near the boundaries experience lower local shear stress intensity, and their impact on global failure mechanisms is less severe. These findings underscore the importance of pore size and geometry and their positioning relative to internal stress fields, particularly in heterogeneous materials such as cemented soils.

Based on the results of numerical and geometrical analyses, it is shown that a pore whose volume does not intersect the shear stress concentration cone (as in position C) exerts a significantly smaller influence than a pore located at or intersecting the boundaries of these high-stress zones. Figure 5b–d presents vertical stress distribution maps for a 1 mm pore positioned at locations A–C, along a vertical cross-section through the centre of the specimen. The corresponding shear strain distributions are shown in Figure 5e–g. At a strain level of ε = 5%, the most remarkable asymmetry in vertical stress is observed for position B. For this configuration, the numerical analysis indicated the earliest onset of cracking.

Figure 6 illustrates the influence of total porosity on the peak compressive strength of cemented soil specimens. The peak strength generally decreases with increasing porosity, which aligns with theoretical expectations. The reference sample with 0% porosity achieved the highest strength (19.15 MPa). At a total porosity of just 1%, the strength dropped to a range of 16–18.3 MPa, representing a maximum reduction of more than 11%.

Figure 6

Relationship between uniaxial compressive strength and total porosity.

With further increases in porosity, the strength followed a downward trend – for example, at approximately 3% porosity, the average strength was around 14.3 MPa. The most porous specimens (3–4.5%) exhibited the lowest strength values, in the 13–15 MPa range, which is approximately 20–30% lower than the solid model. This inverse relationship is consistent with the well-known mechanism of strength reduction in cementitious materials due to a decrease in the load-bearing phase [4].

However, it is essential to note that this dataset’s relationship between strength and porosity is not perfectly monotonic due to variability among the specimens. Some individual samples may even suggest a trend opposite to that reported in the literature. For instance, a specimen with 2.97% porosity reached a peak strength of approximately 15.6 MPa, whereas specimens with 1.25% porosity exhibited lower peak stresses (∼14.44 MPa). Therefore, in mesoscale studies, it is essential to provide detailed information on the shape and position of pores, as these factors can significantly influence the mechanical response.

The deviations shown in Figure 6 result in significant data dispersion. Two trend lines were added to the plot: linear and exponential. Both were constrained to satisfy a boundary condition: at 0% porosity, the trend line was required to pass through the value of 19.15 MPa, corresponding to the peak compressive strength of the solid (pore-free) specimen. Based on the coefficient of determination (R ²) and following the verbal interpretation for linear correlations by Sobczyk [32], a strong correlation is indicated between porosity and maximum compressive strength. The mean squared error (MSE) was also computed for linear and exponential fits, yielding values of 1.88 and 1.70 MPa², respectively.

The numerical results demonstrate that variations in pore shape and pore position, in addition to total porosity, significantly influence the final mesoscale response of the specimens. This indicates that porosity alone is insufficient for a reliable prediction of compressive strength and that additional geometric descriptors of pores must be considered. To quantify the combined effect of these parameters in a statistically consistent and fully interpretable manner, a multiple linear regression model was adopted.

Compressive strength f c {f}_{{\rm{c}}} was treated as the dependent variable, while porosity n n , pore sphericity S , S, and pore position p p served as independent predictors. The regression analysis was performed on the complete dataset (n = 19). The resulting model achieved a high coefficient of determination R 2 = 0.80 {R}^{2}=0.80 and an adjusted coefficient R adj 2 = 0.77 {R}_{{\rm{adj}}}^{2}=0.77 , indicating very good agreement between the numerical simulations and the regression-based prediction despite the limited number of samples.

The final regression equation can be written as: f c = 23.46 − 0.803 n − 4.228 S − 0.660 p , {f}_{{\rm{c}}}=23.46-0.803\hspace{0.25em}n-4.228\hspace{0.25em}{\rm{S}}-0.660\hspace{0.25em}p, where f c {f}_{{\rm{c}}} is expressed in MPa, n n denotes porosity, S S the pore sphericity, and p p the pore position index.

All three predictors were found to be statistically significant at the 0.05 level. Porosity exhibits the strongest and most robust influence on compressive strength (β = −0.803, p = 5.0 × 10⁻⁵), confirming that increasing void content systematically reduces the load-bearing capacity of the cemented soil matrix. This observation is in full agreement with classical micromechanical interpretations, according to which increasing porosity leads to enhanced stress concentration and earlier initiation of damage.

Pore sphericity is also statistically significant (β = −4.228, p = 0.034), demonstrating that pore shape measurably affects the macroscopic strength response. Within the analysed geometric configuration, the results indicate that horizontally elongated pores lead to systematic changes in strength compared to nearly spherical voids. Although the influence of sphericity is weaker than that of porosity, it remains non-negligible and highlights the importance of pore geometry at the mesoscale.

The effect of pore position is likewise statistically significant (β = −0.660, p = 0.0035), confirming that the vertical location of a dominant defect within the specimen substantially influences compressive strength. Pores located closer to the lower part of the specimen were found to be less detrimental than those positioned higher. This observation is consistent with the spatial distribution of stress and deviatoric strain in the numerical models and supports the conclusions drawn in the previous subsection regarding the role of defect location. A summary of the regression parameters describing the mechanical influence of pore-related features on compressive strength is provided in Table 5.

Since the predictors are expressed in different physical units and have different numerical ranges, their relative importance was additionally assessed using standardised regression coefficients. The standardised values indicate that porosity has the strongest relative influence on compressive strength (β* = −0.69), followed by pore position (β* = −0.42) and pore sphericity (β* = −0.28). This ranking confirms that porosity is the primary mesoscale control on compressive strength, while pore position and shape provide additional, independent contributions to the observed variability.

Overall, the regression-based multi-parameter model demonstrates that compressive strength of cemented soil specimens is governed by the combined interaction of porosity, pore shape, and pore position. The high coefficient of determination and the statistical significance of all three predictors confirm that a simple and fully interpretable linear model is sufficient to capture the dominant physical trends present in the analysed dataset.

The comparison of prediction errors demonstrates (Table 6) a clear advantage of the multi-parameter regression approach over single-variable trend-line models. The single-parameter linear and exponential trend-line fits yield relatively high prediction errors (MSE = 1.88 MPa² and MSE = 1.70 MPa², respectively), indicating limited predictive capability when porosity alone is considered. In contrast, the multiple linear regression model, which simultaneously accounts for porosity, pore sphericity, and pore position, achieves a substantially lower prediction error (MSE = 0.34 MPa²). This represents more than a threefold reduction in error compared to the single-parameter models and confirms that compressive strength is governed by the combined interaction of several geometric features rather than porosity alone.