[orcid=0000-0002-7219-249X]

[orcid=0000-0002-6478-1833]

[orcid=000-0001-5734-6238]

[orcid=0000-0001-5557-144X] [1]

label1]addressline=Department of Continuum Mechanics and Structural Analysis,
University of Seville, Camino de los Descubrimientos s/n, 41092 Seville,
Spain.

label2]addressline=Department of Continuum Mechanics, University of Seville,
Avenida Reina Mercedes, 41012 Seville, Spain.

label3]addressline=Department of Structural Mechanics and Hydraulic
Engineering, University of Granada, C/ Dr. Severo Ochoa s/n, 18071 Granada,
Spain.

[cor1]Corresponding author - Department of Structural Mechanics and Hydraulic
Engineering, University of Granada, Av. Fuentenueva sn, 18002 Granada, Spain.
Phone: +34 958241000, Ext. 20668.

# Hierarchical meta-modelling for fast prediction of the elastic properties of
stone injected with CNT/cement mortar

Rubén Rodríguez-Romero<EMAIL_ADDRESS>Víctor Compán<EMAIL_ADDRESS>Andrés
Sáez<EMAIL_ADDRESS>Enrique García-Macías<EMAIL_ADDRESS>[  [  [

###### Abstract

This paper presents an innovative hierarchical meta-modelling approach for
predicting the effective elastic properties of stone injected with composite
cement grouts containing Carbon NanoTubes (CNTs) for structural rehabilitation
purposes. The study first addresses the homogenization of the nano- and micro-
scales through numerical representative volume elements. The first numerical
model focuses on the cement matrix doped with randomly oriented CNTs, while
the second model represents the porous stone with CNT/cement grout filling the
pores. However, the computational burden associated with these models poses a
significant limitation when analysing large-scale macrostructural elements. To
overcome this challenge, a hierarchical meta-modelling approach based on two
nested Kriging surrogate models is proposed. This approach offers a fast and
accurate alternative to bypass the time-consuming numerical homogenization
process. The first surrogate model establishes a connection between the
microstructural characteristics of the CNTs and cement and the effective
elastic properties of the composite grout. Subsequently, a second meta-model
expands over the first one by introducing two additional variables: the
elastic modulus of the stone and the porosity, allowing for estimation of the
overall elastic properties of CNT/cement/stone composites. The accuracy and
efficiency of the proposed hierarchical meta-model is demonstrated through
detailed parametric analyses. To showcase its potential for optimizing
rehabilitation interventions, the paper is closed with its application to a
benchmark case study of a stone column injected with CNT/cement grout.

###### keywords:

Cement grounts Kriging Numerical Homogenization Porous Composites Surrogate
modelling

Composite CNT/cement grouts are investigated for the first in the literature.

A novel nested surrogate modelling approach is proposed for multi-scale
composites.

Fast stiffness estimates of porous stone injected with CNT/cement grout are
obtained.

The proposed approach is used to assess the performance of grout injected
columns.

##  1 Introduction

Masonry structures represent an integral part of Europe’s architectural
heritage and constitute some of the most iconic secular and sacred monumental
buildings. In this light, important research efforts have been devoted to
developing structural rehabilitation techniques pursuing the principles of
minimal intervention and safeguard of the architectural/artistic values [1,
2]. Traditional retrofitting methods include, among others, the addition of
shear or infill walls, steel braces, confinements (ring beams or tie bars),
transversal anchorage systems, steel reinforced grouts, or grout injections
(refer to [3] for a comprehensive state-of-the-art review). In this context,
the recent strides in Nanotechnology and Materials Science have unveiled
unprecedented opportunities for the development of innovative micro- and nano-
engineered materials with vast possibilities across disciplines. Although
still largely unexplored, these advancements offer exciting possibilities for
creating new minimally invasive rehabilitation techniques endowed with
superior mechanical properties.

There is an increasingly frequent use of composite materials for strengthening
and repairing both modern and historic masonry constructions. Especially in
the form of externally bonded composite layers, fiber-reinforced polymers
(FRP) such as carbon fibre and glass fibre reinforced polymers (CFRP/GFRP),
and textile reinforced mortars (TRM) [4, 5, 6] have been increasingly used.
FRP systems have been rapidly adopted due to their ease of application,
versatility, high tensile strength and stiffness-to-weight ratio, as well as
the progressive reduction in their manufacturing and distribution costs.
Nonetheless, such systems suffer from durability issues related to debonding,
poor performance at high temperatures and sensitivity to ultra violet rays, as
well as improper transpiration of the masonry substrate causing compatibility
issues [7]. Another promising field of application of composite materials for
rehabilitation of built heritage regards grout injections. Reinforcing grout
injections consist of low-pressure (0.5–1 atm) injections of fine cement-based
grout through holes drilled into the masonry [8]. This technique represents a
widely accepted approach in rehabilitation [9] for several purposes: to repair
and halt the spread of cracks; to increase resistance to moisture penetration;
and to strengthen historic structures in seismic areas. These grouts are
commonly made of Portland cement (PC) or hydraulic lime as the main binder,
aggregates, and additives such as brick dust, epoxy, pozzolana, bentonite,
ultra-fine slag, and fly ash or clays [10, 11]. However, despite their widely
spread use in masonry rehabilitation, cement-based grouts suffer from some
limitations such as low tensile strength and toughness, being susceptible to
cracking with the subsequent penetration of water and chemicals accelerating
the degradation of the material [12]. Recent developments in Nanotechnology
have opened vast new possibilities for the development of high-performance
composite materials. Notably, nanomaterials such as nano calcium carbonate
(nano-CaCO3) [13], nano silica (nano-SiO2) [14], graphene and derivatives such
as graphene nano-platelets (GnPs) or graphene oxide (GO) [15], and carbon
nanotubes (CNTs) [16] have proved exceptional potential to develop high-
performance multi-functional cement-based composites (refer to [17] for a
comprehensive literature review). Research on cement mortars incorporating
CNTs has been particularly copious in the last years, reporting increases of
around 10%, 25%, and 30% the compressive strength, bending and Young’s moduli
of pristine mortar, respectively, by adding low CNT dosages ($\leq 1$ wt%)
[18, 19]. However, specific research on the development of nano-modified
grouts is not so abundant. One noteworthy contribution from the limited body
of existing literature is the one by Liu and co-authors [20], who investigated
the addition of CNTs for the development of cement composite grouts in
inclined fracture grouting with flowing water. Their results demonstrated
increases in the sealing efficiency of the composite grout from $\approx$60 to
70% for CNT contents increasing from 0.5 to 2%. Another pioneering work was
conducted by Lee et al. [21], who reported the development of self-heating
electrically conductive CNT-cement grouts for the detection of voids in
prestressed concrete through thermal imaging and electrical resistance
analyses.

The modelling of CNT-based composite materials, while critical for optimizing
the composite admixture, poses a formidable challenge given their multi-scale
nature [22]. In the first place, inclusions are defined at the nano-scale,
where the filler-matrix load transfer properties are characterized by inter-
atomic potentials. At this scale, the most common approaches are [23]: (i)
Molecular Dynamics (MD) simulations; (ii) Continuum modelling; and (iii)
atomistic-based continuum approaches. Although MD can analyse high-fidelity
representations of the filler-matrix atomic interaction, their application is
often limited by computational constraints when analysing systems involving
large numbers of atoms. Continuum mechanics-based models can alleviate such
limitations by simulating inclusions through continuum structural elements
(e.g. beams or shells), although their accuracy may be questionable when
adopting concepts from conventional continuum mechanics at the nanoscale [24].
Atomistic-based continuum techniques combine the advantages of both methods by
properly correlating the interatomic potentials with the stiffness of certain
continuum elements such as truss or shell elements [25]. Once the atomic
properties of the composites are determined, a micromechanics approach with an
appropriate upscaling rule is necessary to estimate the
microscopic/macroscopic properties. To this aim, particularly popular
approaches are mean-field homogenization (MFH) methods based on the Eshelby’s
fundamental solution [26]. The extension of these methods after the work by
Mori and Tanaka [27] has proved a computationally efficient solution for
simulating large-scale and complex microstructures in an analytical or semi-
analytical framework. These techniques can incorporate multiple
micromechanical features, including geometrical characteristics of the
inclusions and their orientation distribution, as well as filler waviness,
agglomeration, and interphase effects [28, 29, 30]. For instance, García-
Macías and Castro-Triguero [31] demonstrated, using an extended Mori-Tanaka
model, that the interplay between filler waviness and agglomeration leads to
coupled weakening effects on the overall stiffness of CNT-composites. Rao et
al. [32] developed a double-inclusion MFH model to predict the thermo-
viscoelastic properties of CNT-reinforced polymers, considering imperfect
interfaces. Their findings highlight the significant role of interface effects
in influencing the overall properties of CNT-based composites, emphasizing
that assuming perfect interfaces leads to overestimates of the stiffness of
the composites. These approaches can also be successfully extended to predict
other physical properties of composites, including the effective
piezoresistive [33] or piezoelectric [34] properties, among others. A
noticeable contribution was made by Hassanzadeh-Aghdam and co-authors [35] who
reported the development of an effective medium (EM) micromechanics-based
method to estimate the thermal conductivity of cementitious nanocomposites
doped with CNTs, considering filler waviness and interfacial region
contributions. However, the mechanical interaction between fillers is defined
through simplifying assumptions, with no consideration of the actual micro-
mechanical stress-strain state. To overcome this limitation, computational
homogenization methods offer a more sophisticated yet computationally
intensive approach. These methods consist of simulating a set of virtual tests
on a representative volume element (RVE) of the microstructure discretized
using numerical simulation approaches such as the finite element (FEM) or the
boundary element method (BEM) [36, 37]. In this context, it is noteworthy to
mention the contribution made by Charitos and co-authors [38] who evaluated
the effective elastic properties of CNT/polymer nanocomposites by direct 3D
FEM. This approach relies on the explicit discretization of the composite
phases, requiring advanced mesh strategies and demanding FEM calculations. As
a result, its applicability is limited to fillers with short to moderate
aspect ratios and low volume fractions. In this context, an embedded element
techniques (EETs) have been proposed as a computationally efficient
alternative to simulate more complex fiber-reinforced microstructures. These
approaches simulate fillers as 1-D line elements embedded in the matrix
through certain constraint conditions. In this line, Lu et al. [39] developed
an EET with a two-step searching & coupling technique to constraint the
movement of fillers inside the matrix. That approach allowed discretizing the
matrix phase with a regular non-conforming mesh, enabling the simulation of
complex 3D fiber arrangements with large aspect ratios ($\leq 300$) and volume
fractions (10-40 wt%) with limited computational costs. Other advanced
homogenization techniques comprise numerical simulation approaches accounting
for inclusion-matrix interaction and interfacial effects, Fast Fourier
Transform techniques, damage and fracture mechanics, and multi-scale
homogenization aided by Artificial Intelligence. Interested readers may
consult references [40, 41] for a recent state-of-the-art review.

In light of the literature review conducted above, it is clear that there is
great potential for the development of nano-modified cement grouts. In this
context, the present work presents a numerical investigation of the
possibilities of CNTs to develop high-performance cement-based grouts for the
rehabilitation of weak porous stone. This study builds upon recent research
activities conducted by the authors, involving the assessment of the
structural condition of heritage constructions erected in Andalusia (Spain)
between the 15 and 18${}^{\textrm{th}}$ centuries. Many of these monumental
buildings were made by the so-called San Cristóbal stone, including iconic
constructions such as the Monastery of San Jerónimo de Buenavista [42], the
cathedral of Seville [43], or the Santiago church [44]. The San Cristóbal
stone is characterized by a low elastic modulus (4-10 GPa), compressive
strength (1.5–2.6 MPa), and high porosity (22.4–27%) [45]. As a consequence,
many of these structures have developed severe cracking patterns,
necessitating multiple rehabilitation interventions in recent decades. In this
light, the idea of producing high-performance CNT-reinforced grouts suggests
significant potential for the development of efficient rehabilitation
techniques with minimal impact on architectural values. To this aim, a
hierarchical multi-scale FEM homogenization approach is developed to estimate
the effective elastic properties of porous San Cristóbal stone doped with
CNT/cement grouts. Additionally, to enable the mechanical simulation of large-
scale structural elements with minimal computational costs, the numerical
homogenization model is bypassed by two nested surrogate models, which
represents the main novelty of this work. The effectiveness of the developed
approach and the effect of the main properties of the constituent phases are
assessed through a series of detailed parametric analyses. Finally, the
potential of the proposed approach is illustrated through a demonstration
example evaluating the modal properties of a full-scale column rehabilitated
with grout injections.

##  2 General methodology: hierarchical meta-modelling

The proposed multi-scale modelling approach is sketched in Fig. 1. It spans
across three different scales, including the (a) micro-scale analysing the
interaction between CNTs and the cement matrix of the grout, (b) the meso-
scale studying the interaction between the matrix stone and the CNT/cement
grout occupying the volume of pores and, finally, the (c) macro-scale where
the behaviour of complete structural elements is investigated. The effective
elastic properties from the first two-scales are estimated through numerical
3D RVEs using FEM. Given the considerable computational burden of such models,
which represents a serious difficulty to estimate the response of macroscopic
elements, a new hierarchical meta-model is proposed in this work. This
involves a set of two nested Universal Kriging models $\mathcal{M}_{1}$ and
$\mathcal{M}_{2}$. Firstly, meta-model $\mathcal{M}_{1}$ relates the
constitutive properties of the CNTs and the cement matrix with the effective
elastic properties of the composite grout. At the meso-scale, the RVE is
constituted by spherical inclusions made of CNT/cement composite and a stone
matrix. Herein, the effective properties of the grouted pores are estimated
using $\mathcal{M}_{1}$. On this basis, a second meta-model $\mathcal{M}_{2}$
is constructed to relate the microstructural features of the CNTs, cement, and
the stone matrix with the overall elastic properties of the grouted stone. The
estimates of this second meta-model are finally used to investigate the
macroscopic response of large scale structural elements, and so investigate
the effectiveness of the proposed composite grout for rehabilitation purposes.

Figure 1: Flowchart of the proposed hierarchical meta-modelling: (a)
CNT/cement microstructure, (b) CNT/cement/stone meso-scale, and (c) macro-
scale.

##  3 Theoretical background

###  3.1 Numerical Mechanical homogenization

A classical approach to estimate the effective properties of heterogeneous
materials is to adopt a sub-scale modelling approach based on the concept of a
RVE. This concept lays on the assumption that any material point X in the
continuum scale is associated to a local RVE of volume $\Omega$ and boundary
$\Gamma$ (see Fig. 2). The concept of RVE is only valid in two situations: (i)
unit cell in a periodic microstructure: and (ii) volume containing a
sufficient number of inclusions to achieve statistical homogeneity and
ergodicity in the overall behaviour [46]. On this basis, in order for the sub-
scale-modelling to be energetically consistent, the deformation energy at the
macroscopic level must equate to the volume-averaged micro-scale stress work,
which is encapsulated by the Hill-Mandel condition [47]:

Figure 2: Continuum macro-structure associated to a two-phase RVE consisting
of a matrix phase with stiffness tensor $\textbf{C}_{m}$ and doping inclusions
with stiffness tensor $\textbf{C}_{i}$.

$\overline{\sigma}_{ij}\overline{\varepsilon}_{ij}=\frac{1}{\Omega}\int_{\Omega}\sigma_{ij}\,\varepsilon_{ij}\textrm{d}\Omega,$
(1)

where overbar represents the volume averages of the heterogeneous Cauchy
stress tensor $\bm{\sigma}$ and the (virtual) strain tensor $\bm{\varepsilon}$
given by:

$\overline{\sigma}_{ij}=\frac{1}{\Omega}\int_{\Omega}\sigma_{ij}\textrm{d}\Omega,\quad\overline{\varepsilon}_{ij}=\frac{1}{\Omega}\int_{\Omega}\varepsilon_{ij}\textrm{d}\Omega.$
(2)

In the absence of body forces, the Hill-Mandel principle can be rewritten by
applying the Gauss theorem as follows [48]:

$\int_{\Gamma}\left(t_{i}-\overline{\sigma}_{ij}n_{j}\right)\left(u_{i}-\overline{\varepsilon}_{ik}x_{k}\right)\textrm{d}\Gamma=0,$
(3)

where the terms t and u denote the external traction vector and the
displacement field on the RVE boundary with normal vector n, respectively. On
this basis, to fulfil the Hill-Mandel principle, the boundary conditions (BCs)
of the sub-scale problem must satisfy Eq. (3) in order for the homogenization
to be energetically consistent. The corresponding Boundary Value Problem (BVP)
is typically solved by applying the macro-strain tensor to the RVE. Then, the
macro-stress and macro-strain tensors are directly estimated from Eq. (2) or,
alternatively, through surface integration:

$\begin{split}\overline{\sigma}_{ij}&=\frac{1}{\Omega}\int_{\Gamma}t_{i}x_{j}\textrm{d}\Gamma,\\\
\overline{\varepsilon}_{ij}&=\frac{1}{2\Omega}\int_{\Gamma}\left(u_{i}n_{j}+u_{j}n_{i}\right)\textrm{d}\Gamma.\end{split}$
(4)

Once the proper BCs are applied to the RVE, the effective stiffness tensor
$\mathbb{C}$ can be estimated from the linear relation between the macrostress
tensor $\bm{\overline{\sigma}}$ and the macro-strain tensor
$\bm{\overline{\varepsilon}}$ as:

$\overline{\sigma}_{ij}=\mathbb{C}_{ijkl}\overline{\varepsilon}_{kl}.$ (5)

####  3.1.1 Nanoscale - CNT/cement composite

The first step of the proposed homogenization approach consists of the
estimation of the effective properties of the cement mortar doped with CNTs at
the nano-scale. Let us consider a cubic RVE of edge length $L$ constituted by
a cement matrix doped with randomly oriented CNTs (Fig. 3 (a,b)). The RVE is
designed to comply with four key assumptions: (i) it contains a sufficiently
representative number of CNTs to statistically represent the composite
material as a whole; (ii) the length $L_{cnt}$ and diameter $D_{cnt}$ of the
CNTs are constant; (iii) the fillers are randomly oriented; and (iv) the
mortar matrix can be modelled as a continuum medium. The microstructures are
generated through a set of custom-made scripts developed in the MATLAB
environment. Each CNT is described by its centre point $C$ as a random
variable within the volume of the RVE, and its orientation is defined by means
of two Euler angles $\theta$ and $\varphi$ randomly distributed in the
intervals $\theta\in\left[0,2\pi\right]$ and $\varphi\in\left[0,\pi/2\right]$,
respectively (see Fig. 3 (c)). Fillers are added sequentially to the RVE until
achieving the desired volume fraction, considering non-intersection conditions
between fillers.

Figure 3: Nano-scale: RVE of cement mortar doped with CNTs. (a) Mortar matrix;
(b) CNTs; and (c) Euler angles used to define the orientation of CNTs.

Both the matrix phase and the CNTs have been modelled assuming isotropic
constitutive properties, as a reasonable simplification commonly adopted in
the literature (refer e.g. to [49, 50]). Given the computational difficulties
involved in the modelling of large volume fractions of randomly oriented CNTs
(aspect ratio$\approx$1000) using solid elements, the modified RSA algorithm
proposed in [51, 52] has been adopted in this work. This approach allows
modelling fiber-like inclusions by means of computationally light beam
elements, and the matrix phase by a non-conforming (planar or three-
dimensional) regular mesh (Fig. 3 (a)). On this basis, the deformations of the
fibers and the matrix are constrained considering perfect interfacial bonding.
This alleviates the serious computational demands of complex conforming
meshing, resulting in important reductions in the computational time involved
in the mesh generation. The algorithm comprises two steps for each fibers’
nodes: finding the closest node in the matrix phase using a searching
algorithm, then coupling the translational degrees of freedom (DOFs) of the
identified nodes (Fig. 4). To compensate the augmented stiffness resulting
from the mass matrix covering the volume space that the CNTs should occupy,
the Young’s modulus of the fibers, $E_{f}$, is modified as
$\hat{E}_{f}=E_{f}-E_{m}$ [39], with $E_{m}$ being the elastic modulus of the
matrix phase.

Figure 4: Schematic representation of the ASC algorithm (a), and application
to the FEM of a RVE of cement mortar doped with CNTs (b).

The resulting RVEs have been meshed and simulated using the FEM code ANSYS19
R3. Specifically, CNTs have been discretized using 2-nodes (6 DOFs per node)
3D Timoshenko beam elements (BEAM188) with circular cross-section, while the
mortar matrix has been discretized by a regular array of 8-nodes (3 DOFs per
node) solid elements (SOLID185). Given that geometrical periodicity has not
been imposed because of the difficulties in achieving high volume fractions
when non-intersection conditions among the fibers are considered, the
effective elastic properties of the RVE have been estimated from the BVP in
Eq. (5) after adopting linear displacement BCs (Dirichlet). Such BCs are
defined to be compliant with the Hill-Mandel principle in Eq. (3) by imposing
a set of displacements $u_{i}$ at every point x on the boundary $\Gamma$ of
the RVE such that:

$u_{i}=\overline{\varepsilon}_{ij}x_{j},\quad\forall\textbf{x}\in\Gamma.$ (6)

The stiffness tensor is determined after applying six independent loading
conditions (3 pure dilations and 3 pure distortions) [53, 54] to the RVE. For
each loading case, one of the six components of the macroscopic strain tensor
$\overline{\varepsilon}_{ij}$ is imposed, while the others are set to zero.
Afterwards, the stress field in the heterogeneous material is estimated from
the static equilibrium analysis, and the average stress field components
$\overline{\sigma}_{ij}$ are calculated using Eq. (4). Next, the effective
stiffness tensor $\mathbb{C}$ can be directly derived from a generalised
Hooke’s law in Eq. (5). Note that, despite perfectly random orientation
distributions of CNTs are well known to lead to perfectly isotropic
conditions, the randomness in the construction of the RVE will always induce
certain residual anisotropy. This can be minimized through integration over
all possible orientations in the Euler space, often referred to as the
orientational average. To this aim, the tensor $\mathbb{C}$ can be arbitrarily
rotated by the transformation matrix $\bm{a}$ consisting of $\theta$,
$\varphi$, and $\psi$ Euler angles as
$\mathbb{C}^{*}_{ijkl}(\theta,\varphi,\psi)=a_{ip}a_{jq}a_{kr}a_{ls}\mathbb{C}_{pqrs}$.
On this basis, the orientational average of $\mathbb{C}^{*}$ can be computed
as [55]:

$\langle\mathbb{C}^{*}\rangle=\displaystyle\frac{1}{4\pi^{2}}\int_{0}^{2\pi}\int_{0}^{2\pi}\int_{0}^{\pi/2}\mathbb{C}^{*}(\theta,\varphi,\psi)\sin\varphi\,\textrm{d}\theta\,\textrm{d}\varphi\,\textrm{d}\psi.$
(7)

and the effective elastic modulus $\overline{E}$ and the Poisson’s ratio
$\overline{\nu}$ can be directly extracted from the compliance tensor
$\mathbb{S}=\langle\mathbb{C}^{*}\rangle^{-1}$ (in Mandel-Voigt contracted
notation) as:

$\overline{E}=\dfrac{1}{\mathbb{S}_{11}},\quad\overline{\nu}=-\overline{E}\,\mathbb{S}_{12}.$
(8)

####  3.1.2 Microscale - CNT/cement/stone composite

The second step in the proposed homogenization approach is aimed at estimating
the effective constitutive properties of porous stone injected with CNT/cement
mortar. It is assumed, for the sake of simplicity, that the pores within the
stone are uniformly distributed and possess approximately spherical geometries
with a constant diameter. Under these assumptions, as often adopted in the
literature as a first approximation for porous materials (refer e.g. to [56,
57]), the RVE of the microscale can be defined through a face-centred cubic
cell as shown in Fig. 5 (a).

Figure 5: RVE of CNT/cement/stone composite (a), FEM mesh (b) (porosity =
20%), and schematic representation of periodic BCs (c).

To compute the effective properties of this second RVE, periodic BCs have been
applied. To do so, consider the notation of the cell surfaces $A^{-}/A^{+}$,
$B^{-}/B^{+}$, and $C^{-}/C^{+}$ indicated in Fig. 5 (b). On this basis, the
displacements on a pair of opposite boundary cells (with their normal along
the $x_{j}$ axis) read:

$u_{i}^{K+}=\hat{\varepsilon}_{ij}x_{j}^{K+}+v_{i}^{K+},\quad
u_{i}^{K-}=\hat{\varepsilon}_{ij}x_{j}^{K-}+v_{i}^{K-},$ (9)

with indexes $K+$ and $K-$ indicating the displacements along the positive and
negative $x_{j}$ direction, respectively. The local fluctuations $v_{i}^{K-}$
and $v_{i}^{K+}$ should be identical on every two opposite faces to satisfy
the periodicity of the RVE. Thus, the local displacements can be dropped from
Eq. (9) by the difference between the two expressions, which leads to:

$u_{i}^{K+}-u_{i}^{K-}=\hat{\varepsilon}_{ij}\left(x_{j}^{K+}-x_{j}^{K-}\right).$
(10)

The BCs in Eq. (10) have been implemented in ANSYS through coupling conditions
between the displacements of the nodes on opposite boundaries. For each pair
of nodes with identical in-plane coordinates (up to a certain tolerance) on
two opposite cell faces, the constraint condition among their displacement
components is set according to Eq. (10). By way of example, Figs. (5) (b) and
(c) show the constraint equations to be defined for a pair of nodes located on
the cell faces $A^{-}$ and $A^{+}$. Interested readers can delve into further
in-depth details on the implementation of periodic BCs in reference [58]. On
this basis, likewise the previous model, the components of the effective
elastic tensor can be extracted from from Eq. (5) after adopting six
independent virtual tests (3 pure dilations and 3 pure distortions).

###  3.2 Surrogate modelling: Kriging models

The estimation of the effective properties using the numerical RVEs outlined
in Sections 3.1.1 and 3.1.2 is computationally intensive, which limits their
implementation for simulating full-scale macrostructural elements. To address
this issue, this work proposes the implementation of computationally light
surrogate models to bypass the previous numerical homogenization models.
Specifically, given the unavailability of the source codes of the utilized
commercial FEM software, non-intrusive surrogate models trained on discrete
realizations of the forward FEM models are adopted for simplicity. The
construction of a non-intrusive surrogate model generally involves three
consecutive stages as sketched in Fig. 6 (showing as an example a two-
variables model for simplicity in the representation), including: (i) Sampling
of the design space; (ii) Generation of the training population through Monte
Carlo simulations (MCS) (also referred to as the experimental design (ED));
and (iii) Construction of the prediction surface.

Let us consider $m$ design variables $d_{i}\in\mathbb{R},i=1,\ldots,m$ allowed
to vary only within a certain physically meaningful range
$\left[a_{i},b_{i}\right]$. Accordingly, the vector of design variables
$\textbf{d}=\left[d_{1},\ldots,d_{m}\right]^{\textrm{T}}$ spans the
$m$-dimensional design space
$\mathbb{D}=\left\\{{\textbf{d}\in\mathbb{R}^{m}:a_{i}\leq d_{i}\leq
b_{i}}\right\\}$. To train the surrogate model, it is necessary to assemble a
training population of $N_{s}$ individuals mapping the output $y$ (quantity of
interest) and the design space $\mathbb{D}$, also referred to as the ED. This
is achieved by uniformly sampling the input design space $\mathbb{D}$ to
conform a matrix of design sites
$\textbf{D}=[\textbf{d}^{1},\ldots,\textbf{d}^{N_{s}}]\in\mathbb{R}^{(m\times
N_{s})}$. The corresponding outputs $y^{i}$ are then obtained by direct Monte
Carlo simulations (MCS) of the main FEM (in this case, the homogenization of
the RVEs) and collected in an observation vector
$\textbf{Y}=\left[y^{1},\ldots,y^{N_{s}}\right]^{\textrm{T}}$.

Figure 6: Schematic representation of the work-flow of the construction of a
two-variables non-intrusive surrogate model.

In this study, among the wide range of non-intrusive surrogate models
available in the literature, the Kriging model is chosen due to its remarkable
adaptability to a large variety of problems [59]. The Kriging interpolator
represents the function of interest $y(\textbf{d})$ by combining a regression
term $y_{r}\left(\textbf{d}\right)$ and a zero-mean stochastic process
$\mathcal{Z}\left(\textbf{d}\right)$ as [60]:

$y\left(\textbf{d}\right)=y_{r}\left(\textbf{d}\right)+\mathcal{Z}\left(\textbf{d}\right).$
(11)

The regression function $y_{r}(\textbf{d})$ depends upon $p$ linear regression
parameters
$\bm{\kappa}=\left[\kappa_{1},\ldots,\kappa_{p}\right]^{\textrm{T}}$ and
certain user-defined regression functions
$f(\textbf{x})=\left[f_{1}(\textbf{d}),\ldots,f_{p}(\textbf{d})\right]^{\textrm{T}}$
with $f_{i}:\mathbb{R}^{m}\rightarrow\mathbb{R}$ as:

$y_{r}(\textbf{d})=f(\textbf{d})^{\textrm{T}}\bm{\kappa}.$ (12)

The stochastic process $\mathcal{Z}\left(\textbf{d}\right)$ can be
characterized by a covariance function
$\textrm{Cov}\left[\mathcal{Z}(\textbf{d}_{i})\mathcal{Z}(\textbf{d}_{j})\right]$
between any two data points $\textbf{d}_{i}$ and $\textbf{d}_{j}$:

$\textrm{Cov}\left[\mathcal{Z}(\textbf{d}_{i})\mathcal{Z}(\textbf{d}_{j})\right]=\sigma^{2}r\left(\textbf{d}_{i},\textbf{d}_{j},\bm{\theta}\right),$
(13)

where $\sigma^{2}$ denotes the variance of
$\mathcal{Z}\left(\textbf{d}\right)$, and
$r\left(\textbf{d}_{i},\textbf{d}_{j},\bm{\theta}\right)$ is a user-defined
spatial correlation function dependent on $\bm{\theta}$ hyper-parameters. Once
defined, the Kriging model predicts the response $y\left(\textbf{d}\right)$ at
an arbitrary design site d as:

$\widehat{y}(\textbf{d})=f(\textbf{x})^{\textrm{T}}\bm{\kappa}+r(\textbf{d})^{\textrm{T}}\textbf{R}^{-1}\left[\textbf{Y}-f(\textbf{d})^{\textrm{T}}\bm{\kappa}\right],$
(14)

where R is a $N_{s}\times N_{s}$ correlation matrix with components
$R_{ij}=r\left(\textbf{d}_{i},\textbf{d}_{j},\bm{\theta}\right)$, with
$r(\textbf{d})$ being a vector of correlations between the all the samples in
the DS and d, that is:

$r(\textbf{d})^{\textrm{T}}=\left[r\left(\bm{\theta},\textbf{d}_{1},\textbf{d}\right),\ldots,r\left(\bm{\theta},\textbf{d}_{N_{s}},\textbf{d}\right)\right]^{\textrm{T}}.$
(15)

From Eq. (14), it is clear that the Kriging interpolator is fully defined by
the regression parameters $\bm{\kappa}$ and the correlation parameters
$\bm{\theta}$. We adopt second-order polynomials for the trend term, and
Gaussian correlation functions for the stochastic term as [61]:

$r\left(\textbf{d}_{i},\textbf{d}_{j},\bm{\theta}\right)=\prod_{k=1}^{m}\exp\left[-\theta_{k}\left(d_{i}^{(k)}-d_{j}^{(k)}\right)^{2}\right].$
(16)

For given values of the hyper-parameters $\bm{\theta}$, the trend parameters
$\bm{\kappa}\left(\bm{\theta}\right)$ and the variance
$\sigma^{2}\left(\bm{\theta}\right)$ can be optimally computed in closed-form
by means of the empirical best linear unbiased estimator (BLUE) (consult [60]
for more extensive information). On the other hand, the estimation of the
hyper-parameters $\bm{\theta}$ involves a non-linear optimization problem. In
this work, the maximum-likelihood-estimator is solved through the iterative
pattern search optimization algorithm implemented in the DACE toolbox for
Matlab developed by Lophaven and co-authors [62]. Finally, note that hyper-
parameters $\theta_{k}$ in Eq. (16) dictate the shape of the correlation
function, allowing for potential anisotropy along the dimensions of d.
However, for the sake of simplicity, correlations are assumed isotropic in
this work, that is $\theta_{k}=\theta\quad\forall\,1\leq k\leq m$.

##  4 Numerical results and discussion

The numerical results and discussion of the presented hierarchical meta-
modelling approach for fast estimation of the elastic properties of
CNT/cement/stone composites are presented in this section. The material
parameters employed in the subsequent analyses are furnished in Table 1. It is
important to remark that the values in this table cover the range of the most
commonly used bases in grout injections, spanning from lime to cement mortars
(elastic modulus from 2 to 10 GPa, and Poisson’s ratios from 0.05 to 0.3).

Table 1: Constitutive properties of the considered phases to investigate CNT/cement/stone composite materials. Quantity |  Material
---|---
| MWCNTs |  | Cement mortar |  |  Stone
Young’s modulus |  0.23-1.22 TPaa |  |  2-10 GPac |  |  4-10 GPae
Poisson’s ratio [-] |  0.1-0.28b |  |  0.05-0.3d |  | 0.3
Diameter [nm] | 1.8 |  | - |  | -
Length [µm] | 0.41 |  | - |  | -
Volume fraction [%] | 0.01-1.2 |  | - |  | -
Porosity [%] | - |  | - |  |  5-30f
Sources: a[63]; b[64]; c[65]; d[65]; e[45] ; f[45]

The results hereafter have been organised into validation analyses in Section
4.1, construction of the surrogate models in Section 4.2 and, finally, a toy
example of a stone column to illustrate the potential of the developed
approach for the design of composite injections in Section 4.3.

###  4.1 Validation analyses

The definition and accuracy of the RVE built for CNT/cement composites is
appraised first. A critical aspect of computational homogenization techniques
regards the size of the RVEs as extensively discussed in the literature (see
e.g. [66, 67]). It is well known that when fibers are randomly oriented, their
strengthening effects distribute equally in all directions of space, resulting
in isotropic stiffness tensors (refer e.g. to reference [31]). Therefore, the
determination of the critical size relies on statistical theory, wherein the
departure from isotropy is evaluated as a function of the RVE size. To analyse
the deviation of the effective stiffness tensor from ideal isotropy, we adopt
the anisotropy index [68]:

$\delta_{iso}=\frac{\left|\mathbb{C}-\langle\mathbb{C}^{*}\rangle\right|_{F}}{\left|\mathbb{C}\right|_{F}},$
(17)

where
$\left|\mathbb{A}\right|_{F}=\sqrt{\textrm{tr}{(\mathbb{A}\cdot\mathbb{A}^{\textrm{T}})}}$
denotes the Frobenius norm of a fourth-order tensor $\mathbb{A}$. On this
basis, a perfectly isotropy tensor leads to an anisotropy index
$\delta_{iso}=0$.

Figure 7 (a) reports the convergence analysis conducted to determine the
critical size of the RVE. To do so, the size of the RVE for CNT/cement mortars
is defined in a dimensionless way as $L/L_{cnt}$ ranging between 0.25 and
2.46. Specifically, 6 different RVEs are generated with 5 realizations each to
analyse the dispersion of the estimated effective properties. In these
analyses, the filler volume fraction is fixed to 0.5% and the filler and
matrix phases are meshed considering element sizes of $L_{cnt}/6$ and
$L_{cnt}/18$, respectively. It can be observed in this figure that convergence
is approximately achieved at $L/L_{cnt}=1.96$ $(L=0.80\;\upmu\textrm{m})$.
Beyond this limit, the residual anisotropy remains approximately constant and,
therefore, the RVE size $L/L_{cnt}=1.96$ is fixed in the subsequent analyses.

It is equally important to study the sensitivity of the homogenization results
to discretization density of the RVE to ensure the estimates are mesh-
independent. To do so, different mesh densities have been investigated for the
same RVE ($L=0.80\;\upmu\textrm{m}$, CNT volume fraction 0.5%), ranging
between 15,417 and 899,148 elements. These discretizations correspond to
different mesh densities used in the study. The mesh sizes for the filler and
matrix phases of the RVE range from $L_{cnt}/3$ to $L/9$ and from $L_{cnt}/16$
to $L/48$, respectively. Figure 7 (b) furnishes the convergence analysis
conducted in terms of the normalized effective elastic moduli
$\overline{E}/E_{m}$ and $\overline{G}/G_{m}$ versus the number of finite
elements, $E_{m}$ and $G_{m}$ being respectively the elastic and shear moduli
of the matrix phase. From the results, it is concluded that convergence is
approximately achieved for a RVE discretization of 744,291 elements, which
corresponds to a mesh density of $L_{cnt}/15$ and $L/45$ for the filler and
matrix phases, respectively.

Figure 7: Departure from isotropy versus the RVE size (a) mesh convergence (b)
for the RVE of CNT/cement composites (filler volume fraction of 0.5%). The
error bars in (a) represent the standard deviation of the anisotropy indexes
obtained after five realizations. Figure 8: Validation analysis of the
numerical estimates of the elastic properties of PC/CNT composites considering
fully aligned (a) and randomly oriented filler configurations (b)
(experimental data from [69]), and (c) comparison between the
numerical/analytical predictions of the elastic modulus for against
experimental data for cement/CNT (experimental data from [70]).

The accuracy of the resulting RVE is investigated in Fig. 8 by comparing the
predictions of the computational homogenization with the estimates obtained
using the classical Eshelby-Mori-Tanaka (MT) MFH approach, as well as the
experimental data reported by Talò and co-authors [69] and Kavvadias et al.
[70] for polymer/CNT and cement/CNT composites, respectively. In the analyses
in Figs. 8 (a) and (b), the investigated material corresponds to polycarbonate
(PC) thermoplastic polymer reinforced with single-walled CNTs (SWCNT)
($L_{cnt}=1.3\;\upmu\textrm{m}$, $D_{cnt}=1.8\;\textrm{nm}$). Both phases are
assumed isotropic, with elastic moduli of 2.35 GPa and 0.97 TPa for the PC and
SWCNT, respectively, and Poisson’s ratios of 0.38 and 0.1. Two different
filler configurations are investigated, namely fully aligned and randomly
oriented fillers in Figs. 8 (a) and (b), respectively. In the case of fully
aligned fillers in Fig. 8 (a), 20 different RVEs ($L=0.20\;\upmu\textrm{m}$)
with CNT volume fractions ranging from 0 to 1% and perfectly aligned in the
$x_{1}$ direction (i.e. $\theta=0$ and $\varphi=\pi/2$) have been constructed
and compared with the predictions of the MT model. In this case, the
composites exhibit transversely isotropic properties, with a considerable
stiffer behaviour in the direction of the fibers ($E_{22}=E_{33}<E_{11}$). It
is noted in this case that close agreements are found between the numerical
and the analytic predictions. In the case of randomly oriented fillers in Fig.
8 (b), five different RVEs have been calculated ($L=0.80\;\upmu\textrm{m}$),
and the numerical predictions are benchmarked against both the analytical MT
solution and the experimental data. It is noted that both the estimates
obtained by the MT model and the numerical homogenization are close to the
experimental data for low CNT volume fractions. Nevertheless, considerable
differences arise with the increase in the volume fraction of CNTs. Note,
however, that the numerical estimates are more proximate to the experimental
data, while the MT model leads to considerable overestimates of the effective
properties. The differences between the experimental data and the theoretical
estimates can be attributed to microstructural aspects disregarded in the
homogenization, such as the effects of filler agglomeration and waviness, as
well as imperfect filler-matrix load transfer interactions. These phenomena
have been reported in the literature to act as microstructural defects,
consequently weakening the effective properties of the composite (refer e.g.
to references [69, 71, 72]). Notably, it has been reported in the literature
that as the filler content increases and CNTs become closer to each other, the
likelihood of forming filler rises [73], which is in line with the results
presented in Fig. 8 (b). Moreover, agglomeration effects can lead to
increasing non-linear correlations between the effective elastic modulus and
the CNT content [31], which further supports the need for including these
microstructural features in future developments. In the analyses shown in Fig.
8 (c), the elastic moduli of cement and CNTs are selected as 11.30 GPa and 1.0
TPa, respectively, and the filler aspect ratio is set to 157.89
($L_{cnt}=1.5\;\upmu\textrm{m}$ and $D_{cnt}=9.5$ nm) following reference
[70]. The filler content is presented in this figure as mass fraction (wt),
with the bulk mass densities of cement and CNTs chosen as 1440 and 0.18 g/cm3,
respectively. Note that in this case the numerical results underestimate the
Young’s modulus in in all cases (maximum relative error of 13.53% for
wt=0.1%), except for wt=0.3%. In this case, very close agreements are found,
with slight overestimates by the numerical model. Importantly, these errors
are within the same order of magnitude as those reported in similar numerical
homogenization models previously documented in the literature (e.g. [70, 74]).
Alongside the aforementioned factors, these differences may also be attributed
to CNT-induced alterations in the hydration process of cement and the effects
of ultrasonic sonication on the CNT aqueous solution [75]. Nevertheless, since
the aim of this investigation is to develop a methodology suitable for
conducting a preliminary assessment of the potential of CNT/cement composites
for rehabilitation grouts, the accuracy of the presented formulation is deemed
sufficient for this purpose. The incorporation of further microstructural
details is left for a future work.

Figure 9: Convergence analysis for the RVE of CNT/cement/stone composites (a)
and validation against the theoretical estimates by MT (b) (1% CNT volume
fraction, 500 GPa and 0.1 the elastic modulus and Poisson’s ratio of CNTs, 5
GPa and 0.3 the Young’s modulus and Poisson’s ratio of cement).

Finally, the validation analyses of the RVE developed for CNT/cement/stone
composites are reported in Fig. 9 (a). In these analyses, the elastic moduli
of the CNTs, cement, and stone are kept fixed at 0.5 TPa, 5 GPa, and 4 GPa,
respectively. Additionally, the Poisson’s ratios of the CNTs, cement, and
stone are also maintained constant with values of 0.1, 0.3, and 0.3,
respectively, and the CNT volume fraction is set at 1.0%. Firstly, the
convergence of the mesh density of the RVE is investigated in Fig. 9 (a). In
this case, one single RVE accounting for a porosity of 30% (maximum considered
value) and a size of 0.5 microns is defined and discretized using seven
different mesh densities with a total number of elements from 2,328 to
199,352. To achieve mesh convergence, it is determined that a mesh density of
784,384 elements per unit volume is required for this RVE. Afterwards, the
accuracy of the selected mesh density is benchmarked against the predictions
obtained using the classical MT method in Fig. 9 (b). The results in this
figure reveal that both approaches yield very similar results, although the
predictions by the MT method slightly overestimates those obtained by FEM as
previously reported in the literature (see e.g. [76]).

###  4.2 Surrogate modelling

Figure 10: Convergence analysis of metamodels $\mathcal{M}_{1}$ (a,b) and
$\mathcal{M}_{2}$ (c,d) in terms of effective elastic modulus (a,c) and
Poisson’s ratio (b,d) random for increasing training populations.

The first surrogate model, $\mathcal{M}_{1}$, considers five different input
variables, namely the CNTs volume fraction, the elastic modulus and Poisson’s
modulus of CNTs, and the elastic modulus and Poisson’s ratio of the cement
matrix. On the other hand, the second surrogate model, $\mathcal{M}_{2}$, adds
two additional input variables, comprising the porosity and Young’s modulus of
the stone. The variation ranges of these parameters defining the design space
$\mathbb{D}$ have been set according to the ranges previously reported in
Table 1. The accuracy of any surrogate model is largely determined by the
quality of the ED. The quantity and distribution of training samples is highly
case-dependent, often necessitating specific analyses to adapt the ED to the
variability of the quantity of interest. Typically, it is necessary to achieve
a uniform sampling of the design space to cover the entire domain of interest.
To this aim, Latin Hypercube Sampling has been adopted in this work to draw
uniformly distributed samples over $\mathbb{D}$. To ascertain the appropiate
size of the ED, convergence analyses have been performed considering varying
training populations of increasing size. For $\mathcal{M}_{1}$, five different
training populations have been considered, ranging from 32 to 512 individuals,
with a validation set (VS) of 1024 samples. For $\mathcal{M}_{2}$, six
training populations have been defined, ranging from 32 to 1024 individuals,
with a VS of 2048 samples. For every population, the effective Young’s modulus
and Poisson’s ratio of each individual are obtained by forward evaluation of
the numerical RVEs previously outlined in Section 3.1. Note that, likewise any
other non-intrusive meta-model model bypassing a computationally intensive
model, the generation of the training population represents the most
computationally intensive step. Then, the convergence of the accuracy of the
predictions by the surrogate-model is evaluated in Fig. 10 in terms of the
Root Mean Square Error (RMSE) and the coefficient of determination R2,
measured over the VS, between the predictions of the surrogate models and the
forward FEM models. As expected, Fig. 10 shows that the convergence rate of
$\mathcal{M}_{2}$ is slower than that of $\mathcal{M}_{1}$ due to the larger
number of input variables. Based on these results, the convergence of meta-
models $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ is approximately achieved for
population sizes of 27 and 29 individuals, respectively. In these cases, RMSE
values close to zero and R2 coefficients close to one ensure the high
performance of the surrogate models.

The surrogate models $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$, trained with 27
and 29 individuals respectively, are selected for the subsequent analyses.
Figure 11 depicts the comparison between the predictions of the forward RVEs
and the resulting surrogate models. These plots, commonly used in surrogate
modelling applications, provide a visual representation of the goodness-of-fit
between the predictions of the meta-models and the forward FEM models. The
minimal dispersion of the data points along the diagonal lines (R2 values very
close to 1) in Fig. 11 provides strong evidence that the surrogate models are
trained with high precision. Both meta-models also exhibit low computational
burden, with mean evaluation times of 0.18 seconds and 0.25 seconds for meta-
models $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$, respectively. These reduced
computational evaluation times make the meta-models an efficient solution for
conducting intense parametric analyses, iterative optimization, as well as the
analysis of full-scale structural elements.

Figure 11: Scatter plots of the predictions by the Kriging meta-models
$\mathcal{M}_{1}$ (a,b) and $\mathcal{M}_{2}$ (c,d) versus the forward
evaluations of the FEM RVEs for the VSs of 1024 and 2048 samples.

###  4.3 Parametric analyes

Figure 12: Variation of the effective elastic modulus (a) and Poisson’s ratio
(b) of CNT/cement/stone composites versus the CNT volume fraction for varying
porosity percentages in stone.

In order to showcase the potential of the proposed nested meta-models, a
comprehensive set of parametric analyses evaluating the impact of the input
variables on the effective elastic properties of CNT/cement/stone composites
is presented herein. The effects of the CNT volume fraction and the porosity
in stone are first investigated in Fig. 12. In these analyses, the elastic
moduli of the CNTs, cement, and stone are kept fixed as 1.0 TPa, 10 GPa and 4
GPa, respectively, and the Poisson’s ratio of the CNTs and cement are also
maintained as constant with values of 0.1 and 0.3, respectively. Furthermore,
the CNT volume fraction is systematically varied from 0.0 to 1.2%, and 5
levels of porosity are considered in the stone ranging from 5% to 30%. On this
basis, the variations of the elastic modulus and Poisson’s ratio of the
CNT/cement/stone composites are presented in Figs. 12 (a) and (b),
respectively. It is noted in Fig. 12 (a) that the effective elastic modulus
monotonically increases as so does the CNT volume fraction, with more
pronounced gradients as the porosity increases. From a physical point of view,
this indicates that a greater amount of CNTs in the composite grout occupying
the porosity in the stone contributes to the overall rigidity of the grouted
stone. The opposite behaviour can be observed in Fig. 12 (b) for the effective
Poisson’s ratio, which is ascribed to the comparatively lower Poisson’s ratio
of the CNTs compared to the San Cristóbal stone (see Table 1). Interestingly,
the Poisson’s ratio increases for higher porosity levels, which is conceivably
explained by the larger content of cement grout occupying the porosity and
with a similar ratio to stone’s. These analyses are complemented with a second
parametric analysis in Fig. 13. This figure presents the analysis of the
effects of the porosity level in stone on the overall elastic properties of
CNT/cement/stone composites. Specifically, the porosity is systematically
varied within the range from 5 to 30%, as typically observed in the San
Cristóbal stone. Furthermore, five cases of CNT volume fractions are
considered from 0% to 1.2%. The results in Fig. 13 (a) reveal an increase in
the effective elastic modulus with increasing porosity, this trend being more
acute for a higher CNT volume fractions. Similarly, the results in Fig. 13 (b)
reveal increases in the effective Poisson’s ratio as the porosity increases.
These results agree with the previously reported ones in Fig. 12, finding
lower ratios for higher CNT volume contents.

Figure 13: Variation of the effective elastic modulus (a) and Poisson’s ratio
(b) of CNT/cement/stone composites versus the porosity content for varying CNT
volume fractions.

Finally, Fig. 14 summarizes the conducted parametric analyses by presenting
the surfaces of the effective elastic moduli of CNT/cement/stone composites as
functions of the CNT volume fraction and porosity of the stone (a), and versus
the CNT volume fraction and Young’s modulus (b). These figures demonstrate
that the potential of CNT/cement grouts for increasing the rigidity of stone
increases as so does the elastic modulus of CNTs and the porosity of the
original stone, achieving maximum increases of around 35%. These results are
reasonable since, as the porosity of the stone increases, so does the volume
occupied by the composite grout. In addition, the stiffness of the grout
increases as so does the elastic modulus and the volume fraction of the CNTs.
Nonetheless, it is important to remark the simplifying assumptions adopted in
this work. On one hand, we are not considering filler agglomeration effects,
which may have a more significant influence as the CNT volume fraction
increases. Extensive literature reports have consistently documented that
CNTs, due to their large specific surface area, tend to form agglomerates
induced by van der Waals attraction forces. This gives rise to CNT clusters
that act as microstructural defects, which adversely impact the overall
effective mechanical properties of the composite [71, 77]. Another important
aspect disregarded in this work includes the connectivity of the pores and
permeability of the stone. These factors critically determine the penetration
potential of the grout and, therefore, the proportion of pores that can
potentially be occupied by the composite grout. While the experience of the
authors demonstrates that the San Cristóbal stone typically presents highly
interconnected porosity, making the assumptions adopted in this work
reasonable, specific studies in this regard are required to support these
analyses. Due to space constraints and the considerable fundamental
developments required for such studies, the extension of the presented
framework to accommodate filler agglomeration effects and partially grouted
pores is left for future work.

Figure 14: Effective elastic modulus of CNT/cement/stone composites as a
function of the CNT volume fraction and porosity (a), and CNT volume fraction
and elastic modulus (b) (porosity=20%).

###  4.4 Application case study: injected stone column

This final case study is intended to illustrate the potential of the proposed
methodology for the macroscopic analysis of full-scale structures. The example
consists of a hypothetical column inspired by one of the columns located in
the Alameda de Hércules square in Seville. Its dimensions include a total
height of 9.20 m, a diameter of 95 cm at the base and 70 cm under the capital.
As a potential rehabilitation intervention, this section investigates the
effect of injecting the column with CNT/cement grout along its with a
perforation of 7.62 cm of diameter as sketched in Fig. 15 (a). To this aim,
the constitutive properties of the grout and the CNT/cement/composite are
estimated using the developed meta-models $\mathcal{M}_{1}$ and
$\mathcal{M}_{2}$. The properties of the stone used in the analysis include a
Young’s modulus value of 4 GPa, a Poisson’s ratio of 0.3, an an porosity of
20%. It is assumed that the perforation is filled with CNT/cement grout
($\mathcal{M}_{1}$), while the remaining section of the column is fully
saturated by the injection ($\mathcal{M}_{2}$). The column has been modelled
and discretized in the FEM software ABAQUS using 4-nodes (Q3D4) tetrahedron
solid elements (total number of 308,407 elements and 60,387 nodes). Cantilever
boundary conditions are applied by restraining only the translational degrees
of freedom at the base of the column. On this basis, a set of parametric modal
analyses are performed to evaluate the stiffening effects of the composite
grout in terms of the natural frequencies of vibration of the column. The
variations are quantified as relative changes in the natural frequencies
compared to the non-injected configuration of the column (Fig. 15 (b)). The
analysis includes the first five natural modes: two first order bending modes
(symmetric along the main directions of the column), two second order bending
modes, and one first order torsional mode. Figure 16 (a) furnishes the
relative variations of the resonant frequencies for CNT filler contents
ranging from 0 to 1.2%. It is noted in this figure that the use of standard
cement grouts (CNT volume fraction equal zero) leads to an increase in the
fundamental frequency of approximately 36%. Then, when considering small
contents of CNTs (1.2% the volume of cement), this increase can rise up to
around 40%. It is also interesting to investigate the effect of the elastic
modulus of cement as reported in Fig. 16 (b). This figure demonstrates that
the stiffness of the cement plays a critical role, with the fundamental
frequency increasing from 24% to 40% as the elastic modulus of the cement
varies between 2 to 10 GPa. It is important to emphasize that the whole
parametric analysis, facilitated by the developed surrogate model, were
completed in just 27 minutes. In contrast, approaching this problem via a
numerical multi-scale approach would present overwhelming, and possibly
unfeasible, computational challenges. Overall, these results highlight the
potential of the develop approach for designing rehabilitation interventions
using CNT-reinforced cement grouts. They also show great promise as a
foundational framework for incorporating complex factors such as filler
agglomeration and grout injection permeability.

Figure 15: Benchmark case study: injected stone column (a); reference modes of
vibration of the non-injected column (b). Figure 16: Parametric analysis of
the natural frequencies of vibration of the toy column injected with
CNT/cement grout considering increasing CNT volume fractions (elastic modulus
of cement equal to 10 GPa) (a) and elastic modulus of the cement for a CNT
content of 1.2% (b) (the porosity of the stone is fixed to 20%)

##  5 Concluding remarks

We have presented in this work an innovative hierarchical meta-modelling
approach for the rapid prediction of the effective elastic properties of stone
injected with CNT/cement grouts. The developed methodology utilizes two nested
Kriging meta-models to bypass the numerical homogenization of the two main
scales involved in CNT/cement/stone composites with minimal computational
burden. Specifically, the nano-scale is defined by a two-phase RVE comprising
a cement matrix doped with randomly arranged CNTs. At the micro-scale, a
second RVE is defined to represent porous stone injected with CNT/cement
grout. The meta-models are trained in a non-intrusive fashion through direct
Monte Carlo simulations of the numerical homogenization of the RVEs.
Specifically, the first meta-model, $\mathcal{M}_{1}$, considers as
independent variables the CNT volume fraction, as well as the Young’s moduli
and Poisson’s ratios of the matrix and the CNTs (5 variables). The second
meta-model, $\mathcal{M}_{2}$, expands upon $\mathcal{M}_{1}$ by adding the
porosity and the Young’s modulus of the stone as additional independent
variables. The dimension and mesh density of the RVEs, as well as the sizes of
the training populations for the surrogate models, are selected through
convergence and statistical analyses. The soundness and effectiveness of the
proposed hierarchical meta-modelling approach have been appraised through
thorough parametric analyses. The paper has been concluded with a toy example
of a stone column, demonstrating the effectiveness of the investigated
CNT/cement grouts for structural rehabilitation. The main discoveries of this
work can be succinctly summarized as follows:

  * •

The conducted convergence analyses have revealed that number of training
points required for achieving accurate predictions by the meta-model is
directly related to the number of independent input variables. Specifically,
the convergence for meta-model $\mathcal{M}_{1}$ (5 variables) was achieved
for $2^{7}$ individuals, while meta-model $\mathcal{M}_{2}$ (7 variables)
required $2^{9}$ samples. Notably, the meta-models demonstrated high accuracy
and very low computational demands, with evaluation times of 0.18 and 0.25
seconds for meta-models $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$, respectively.

  * •

The reported parametric analyses have showcased the effect of key
microstructural features upon the elastic moduli of CNT/cement/stone
composites, namely the CNT volume fraction, the elastic properties of CNT,
cement and stone, as well as the porosity in the stone. The presented
numerical results have revealed the potential increase of up to 35% the
elastic modulus of porous stone (30% of porosity) with a CNT volume fraction
of 1.2% the volume of cement in the grout.

  * •

The presented toy example of a 3D stone column has demonstrated the potential
of CNT/cement grout for rehabilitation purposes. To this aim, the hierarchical
meta-modelling approach presented this study offers a fast and accurate method
for optimizing such interventions.

It is important to emphasize that several simplifying assumptions have been
considered in the definition of the microstructures for the sake of simplicity
in the construction of the nested surrogate models, which constitute the
actual innovative contribution of this work. In particular, the inclusion of
additional microstructural features, such as filler waviness, agglomeration,
and non-uniform filler aspect ratio and orientation distributions should be
considered. These additions would enhance the model’s accuracy and capture
more realistic behaviours. Moreover, accounting for the micro-constituent
phases and the porosity of cement, and utilizing more realistic geometrical
models for the porosity of stone may significantly improve the accuracy of the
proposed model. A potentially valuable extension within the proposed
methodology involves the incorporation of non-linear load-bearing links
between the CNTs and the matrix phase, with the purpose of predicting the
ultimate strength properties of the CNT-reinforced grouts. Finally, analysing
the permeability of stone to grout injections would be valuable for a
comprehensive optimization of the proposed grout technology. Nonetheless,
given the general nature and flexibility of the proposed hierarchical meta-
modelling approach, its practical applications and potential for further
development render it a promising methodology in the field of composite
materials. Overall, the proposed approach provides a rapid and accurate method
for predicting the effective elastic properties of stone injected with
CNT/cement composite grout, enabling the structural analysis of full-scale
elements while incurring minimal computational costs.

## Acknowledgements

This work was supported by the Operational Programme 2014-2020 of the European
Regional Development Fund and the Consejería de Economía, Conocimiento,
Empresas y Universidad of the Andalusian Regional Government (Spain) under
project US-1381164. The financial support is gratefully acknowledged. R. R-R
was also supported by a FPU contract-fellowship from the Spanish Ministry of
Education Ref: FPU18/05211.

## References

  * [1] Alessia Buda, Virginia Gori, Ernst Jan de Place Hansen, Cristina S Polo López, Valentina Marincioni, Emanuela Giancola, Nathalie Vernimme, Aitziber Egusquiza, Franziska Haas, and Daniel Herrera-Avellanosa.  Existing tools enabling the implementation of en 16883: 2017 standard to integrate conservation-compatible retrofit solutions in historic buildings.  Journal of Cultural Heritage, 57:34–52, 2022.
  * [2] P. D. Gkournelos, T. C. Triantafillou, and D. A. Bournas.  Seismic upgrading of existing masonry structures: A state-of-the-art review.  Soil Dynamics and Earthquake Engineering, 161:107428, 2022.
  * [3] Fahimeh Yavartanoo and Thomas H. K. Kang.  Retrofitting of unreinforced masonry structures and considerations for heritage-sensitive constructions.  Journal of Building Engineering, 49:103993, 2022.
  * [4] Riccardo Liberotti, Federico Cluni, Francesco Faralli, and Vittorio Gusella.  Natural FRCM and Heritage Buildings: Experimental Approach to Innovative Interventions on “Wall Beams”.  Buildings, 12(8):1076, 2022.
  * [5] Jianyang Xue, Huaiquan Ling, Xiangbi Zhao, Fengliang Zhang, and Hanliang Zhou.  Seismic performance on stone masonry loess cave retrofitted using composite materials.  Journal of Building Engineering, 46:103704, 2022.
  * [6] F Cannizzaro, B Pantò, S Caddemi, and I Caliò.  Discrete modelling of externally bonded composite layers on masonry structures.  Composite Structures, 315:116937, 2023.
  * [7] Samuel A Babatunde.  Review of strengthening techniques for masonry using fiber reinforced polymers.  Composite Structures, 161:246–255, 2017.
  * [8] Marco Corradi, Cristina Tedeschi, Luigia Binda, and Antonio Borri.  Experimental evaluation of shear and compression strength of masonry wall before and after reinforcement: Deep repointing.  Construction and building materials, 22(4):463–472, 2008.
  * [9] ICOMOS.  Recommendations for the analysis, conservation and structural restoration of architectural heritage, 2003.
  * [10] Mohammad Reza Azadi, Ali Taghichian, and Ali Taheri.  Optimization of cement-based grouts using chemical additives.  Journal of Rock Mechanics and Geotechnical Engineering, 9(4):623–637, 2017.
  * [11] Fei Sha, Shucai Li, Rentai Liu, Zhaofeng Li, and Qingsong Zhang.  Experimental study on performance of cement-based grouts admixed with fly ash, bentonite, superplasticizer and water glass.  Construction and Building Materials, 161:282–291, 2018.
  * [12] Joshua Hogancamp and Zachary Grasley.  The use of microfine cement to enhance the efficacy of carbon nanofibers with respect to drying shrinkage crack resistance of portland cement mortars.  Cement and Concrete Composites, 83:405–414, 2017.
  * [13] Faiz U. A. Shaikh and Steve W. M. Supit.  Mechanical and durability properties of high volume fly ash (HVFA) concrete containing calcium carbonate (CaCO3) nanoparticles.  Construction and building materials, 70:309–321, 2014.
  * [14] George Quercia, Götz Hüsken, and H. J. H. Brouwers.  Water demand of amorphous nano silica and its impact on the workability of cement paste.  Cement and concrete research, 42(2):344–357, 2012.
  * [15] Kelly Yee and Mergen H Ghayesh.  A review on the mechanics of graphene nanoplatelets reinforced structures.  International Journal of Engineering Science, 186:103831, 2023.
  * [16] Mahyar Ramezani, Ayoub Dehghani, and Muhammad M Sherif.  Carbon nanotube reinforced cementitious composites: A comprehensive review.  Construction and Building Materials, 315:125100, 2022.
  * [17] Mustafa Şahmaran, Faiz Shaikh, and Gürkan Yıldırım, editors.  Recent Advances in Nano-Tailored Multi-Functional Cementitious Composites.  Woodhead Publishing Series in Civil and Structural Engineering. Woodhead Publishing, 2022.
  * [18] S.M.A. El-Gamal, F.S. Hashem, and M.S. Amin.  Influence of carbon nanotubes, nanosilica and nanometakaolin on some morphological-mechanical properties of oil well cement pastes subjected to elevated water curing temperature and regular room air curing temperature.  Construction and Building Materials, 146:531–546, 2017.
  * [19] K.M. Liew, M.F. Kai, and L.W. Zhang.  Carbon nanotube reinforced cementitious composites: An overview.  Composites Part A: Applied Science and Manufacturing, 91:301–323, 2016.
  * [20] Yaohui Liu, Lei Zhang, Ping Yang, and Xiao Wei.  Experimental investigation on sealing efficiency of CNT composite grouts in inclined fractures with flowing water.  Construction and Building Materials, 270:121494, 2021.
  * [21] Heeyoung Lee, Donghoon Kang, Jongwon Kim, Kyungwho Choi, and Wonseok Chung.  Void detection of cementitious grout composite using single-walled and multi-walled carbon nanotubes.  Cement and Concrete Composites, 95:237–246, 2019.
  * [22] Wenqin Han, Jinyu Zhou, and Qinghe Shi.  Research progress on enhancement mechanism and mechanical properties of FRP composites reinforced with graphene and carbon nanotubes.  Alexandria Engineering Journal, 2022.
  * [23] Roham Rafiee and Reza Maleki Moghadam.  On the modeling of carbon nanotubes: A critical review.  Composites Part B: Engineering, 56:435–449, 2014.
  * [24] Tienchong Chang, Jingyan Geng, and Xingming Guo.  Prediction of chirality- and size-dependent elastic properties of single-walled carbon nanotubes via a molecular mechanics model.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 462(2072):2523–2540, apr 2006.
  * [25] Enrique García-Macías, Carlos Felipe Guzmán, Erick I Saavedra Flores, and Rafael Castro-Triguero.  Multiscale modeling of the elastic moduli of CNT-reinforced polymers and fitting of efficiency parameters for the use of the extended rule-of-mixtures.  Composites Part B: Engineering, 159:114–131, 2019.
  * [26] Toshio Mura.  Micromechanics of defects in solids, volume 3.  Springer, 1987.
  * [27] Tanaka Mori and Kohichi Tanaka.  Average stress in matrix and average elastic energy of materials with misfitting inclusions.  Acta metallurgica, 21(5):571–574, 1973.
  * [28] Stefano Giordano.  Order and disorder in heterogeneous material microstructure: Electric and elastic characterisation of dispersions of pseudo-oriented spheroids.  International Journal of Engineering Science, 43(13-14):1033–1058, 2005.
  * [29] Enrique García-Macías, Rafael Castro-Triguero, and Filippo Ubertini.  Two-step hierarchical micromechanics model of partially saturated porous composites doped with ellipsoidal particles with interface effects.  Composites Part B: Engineering, 148:49–60, 2018.
  * [30] Mohammad Kazem H. A., Reza Ansari, and Abolfazl Darvizeh.  Micromechanical analysis of carbon nanotube-coated fiber-reinforced hybrid composites.  International Journal of Engineering Science, 130:215–229, 2018\.
  * [31] Enrique García-Macías and Rafael Castro-Triguero.  Coupled effect of CNT waviness and agglomeration: A case study of vibrational analysis of CNT/polymer skew plates.  Composite Structures, 193:87–102, 2018.
  * [32] Y. N. Rao and H. L. Dai.  Micromechanics-based thermo-viscoelastic properties prediction of fiber reinforced polymers with graded interphases and slightly weakened interfaces.  Composite Structures, 168:440–455, 2017.
  * [33] Federico C Buroni and Enrique García-Macías.  Closed-form solutions for the piezoresistivity properties of short-fiber reinforced composites with percolation-type behavior.  Carbon, 184:923–940, 2021.
  * [34] N Mishra and K Das.  A Comparative Study of Incremental Selfconsistent and Eshelby–Mori–Tanaka Models for Estimating the Electroelastic Properties of Piezoelectric Polymer Composites with an Orthotropic Matrix.  Mechanics of Composite Materials, 58(5):657–672, 2022.
  * [35] M. K. Hassanzadeh-Aghdam, R Ansari, MJ Mahmoodi, A Darvizeh, and A Hajati-Modaraei.  A comprehensive study on thermal conductivities of wavy carbon nanotube-reinforced cementitious nanocomposites.  Cement and Concrete Composites, 90:108–118, 2018.
  * [36] A Trofimov, A Markov, S. G. Abaimov, I Akhatov, and I Sevostianov.  Overall elastic properties of a material containing inhomogeneities of concave shape.  International Journal of Engineering Science, 132:30–44, 2018.
  * [37] C Naili, I Doghri, T Kanit, M. S. Sukiman, A. Aissa-Berraies, and A. Imad.  Short fiber reinforced composites: Unbiased full-field evaluation of various homogenization methods in elasticity.  Composites Science and Technology, 187:107942, 2020.
  * [38] Ilias Charitos, Anastasios Drougkas, and Evagelia Kontou.  Prediction of the elastic modulus of LLDPE/CNT nanocomposites by analytical modeling and finite element analysis.  Materials Today Communications, 24:101070, 2020.
  * [39] Zixing Lu, Zeshuai Yuan, and Qiang Liu.  3D numerical simulation for the elastic properties of random fiber composites with a wide range of fiber aspect ratios.  Computational Materials Science, 90:123–129, 2014.
  * [40] Ahmed Elmasry, Wiyao Azoti, Sherif A El-Safty, and Ahmed Elmarakbi.  A comparative review of multiscale models for effective properties of nano-and micro-composites.  Progress in Materials Science, page 101022, 2022.
  * [41] Aref Abbasi Moud.  Recent advances in utility of artificial intelligence towards multiscale colloidal based materials design.  Colloid and Interface Science Communications, 47:100595, 2022.
  * [42] F Pozo, R Gil, and E Méndez.  Intervención arquitectónica puntual en Monasterio de San Jerónimo de Buenavista (Sevilla).  Consejería de Urbanismo, 2003.
  * [43] Alfonso Jiménez Martín.  El proyecto de sustitución de los pilares 4C y 5C de la Catedral de Sevilla.  Technical report, Ministerio de Cultura y Deporte, 2009.
  * [44] A. Barrios Padura, J. Barrios Sevilla, and J. García Navarro.  Bearing capacity diagnosis of Santiago church (Jerez de la Frontera, Spain).  Construction and Building Materials, 25(5):2519–2527, 2011.
  * [45] Juan Ramón Baeza, Víctor Compán, Germán Castillo, Margarita Cámara, and Pablo Pachón.  Determining static elastic modulus of weak sandstone in Andalusian historical constructions from non-destructive tests: San Cristóbal’s stone.  Journal of Building Engineering, 57:104864, 2022.
  * [46] Martin Ostoja-Starzewski.  Microstructural randomness versus representative volume element in thermomechanics.  Journal of Applied Mechanics, 69(1):25–35, 2002.
  * [47] Rodney Hill.  Elastic properties of reinforced solids: some theoretical principles.  Journal of the Mechanics and Physics of Solids, 11(5):357–372, 1963\.
  * [48] V. D. Nguyen, Eric Béchet, Christophe Geuzaine, and Ludovic Noels.  Imposing periodic boundary condition on arbitrary meshes by polynomial interpolation.  Computational Materials Science, 55:390–406, 2012.
  * [49] Stelios K Georgantzinos, Panagiotis Antoniou, Stylianos Markolefas, and Georgios Giannopoulos.  Finite element predictions on vibrations of laminated composite plates incorporating the random orientation, agglomeration, and waviness of carbon nanotubes.  Acta Mechanica, 233(5):2031–2059, 2022.
  * [50] Shuai Zhang, Yulin Lai, Kuanyu Chen, Mostafa Habibi, Masoud Khorami, and Zainab Haider Mussa.  Influence of MWCNT’s waviness and aggregation factors on wave dispersion response of MWCNT-strengthened nanocomposite curved beam.  In Structures, volume 53, pages 1239–1249. Elsevier, 2023.
  * [51] Lucian Iorga, Yi Pan, and Assimina Pelegri.  Numerical characterization of material elastic properties for random fiber composites.  Journal of Mechanics of Materials and Structures, 3(7):1279–1298, 2008.
  * [52] Yi Pan, Lucian Iorga, and Assimina A Pelegri.  Analysis of 3D random chopped fiber reinforced composites using FEM and random sequential adsorption.  Computational Materials Science, 43(3):450–461, 2008.
  * [53] Chang Karen Yan.  On homogenization and de-homogenization of composite materials.  Drexel University, 2003.
  * [54] José Carlos García-Merino, Carmen Calvo-Jurado, and Enrique García-Macías.  Surrogate modeling of the effective elastic properties of spherical particle-reinforced composite materials.  Journal of Mathematical Chemistry, 60(8):1555–1570, 2022.
  * [55] J Schjødt-Thomsen and Ryszard Pyrz.  The Mori–Tanaka stiffness tensor: diagonal symmetry, complex fibre orientations and non-dilute volume fractions.  Mechanics of Materials, 33(10):531–544, 2001.
  * [56] R. Grasset-Bourdel, A. Alzina, N. Tessier-Doyen, N. Schmitt, M. Huger, T. Chotard, D. Gruber, and H. Harmuth.  Optimisation of 3D RVE for anisotropy index reduction in modelling thermoelastic properties of two-phase composites using a periodic homogenisation method.  Computational Materials Science, 50(11):3136–3144, 2011.
  * [57] Yingmin Zhang, Guang Yang, Wenwu Chen, and Lizhi Sun.  Relation between Microstructures and Macroscopic Mechanical Properties of Earthen-Site Soils.  Materials, 15(17):6124, 2022.
  * [58] Harald Berger, Sreedhar Kari, Ulrich Gabbert, Reinaldo Rodríguez-Ramos, Julián Bravo-Castillero, and Raúl Guinovart-Díaz.  A comprehensive numerical homogenisation technique for calculating effective coefficients of uniaxial piezoelectric fibre composites.  Materials Science and Engineering: A, 412(1-2):53–60, 2005.
  * [59] Jack P. C. Kleijnen.  Kriging metamodeling in simulation: A review.  European Journal of Operational Research, 192(3):707–716, 2009\.
  * [60] Jack P. C. Kleijnen.  Regression and Kriging metamodels with their experimental designs in simulation: a review.  European Journal of Operational Research, 256(1):1–16, 2017.
  * [61] Jerome Sacks, William J. Welch, Toby J. Mitchell, and Henry P. Wynn.  Design and analysis of computer experiments.  Statistical science, 4(4):409–423, 1989.
  * [62] Søren Nymand Lophaven, Hans Bruun Nielsen, and Jacob Søndergaard.  DACE-A Matlab Kriging toolbox, version 2.0.  Technical report, Technical University of Denmark, 2002.
  * [63] Ahmed El Moumen, Mostapha Tarfaoui, and Khalid Lafdi.  Computational homogenization of mechanical properties for laminate composites reinforced with thin film made of carbon nanotubes.  Applied Composite Materials, 25:569–588, 2018.
  * [64] P Zhao and G Shi.  Study of Poisson ratios of single-walled carbon nanotubes based on an improved molecular structural mechanics model.  Computers Materials and Continua, 22(2):147, 2011.
  * [65] Anastasios Drougkas, Els Verstrynge, Roald Hayen, and Koenraad Van Balen.  The confinement of mortar in masonry under compression: experimental data and micro-mechanical analysis.  International Journal of Solids and Structures, 162:105–120, 2019\.
  * [66] Lili Chen, Boqin Gu, Jiahui Tao, and Jianfeng Zhou.  The average response and isotropy of 3D representative volume elements for random distributed short fibers reinforced elastomer.  Composite Structures, 216:279–289, 2019.
  * [67] Wenjun Bai, Zuer Gong, Yuan Li, and Jinxiang Liu.  Determination of the representative volume element model critical size for carbon fiber reinforced polymer composites.  Composites Science and Technology, 234:109946, 2023.
  * [68] O Zerhouni, MG Tarantino, and Kostas Danas.  Numerically-aided 3D printed random isotropic porous materials approaching the Hashin-Shtrikman bounds.  Composites Part B: Engineering, 156:344–354, 2019.
  * [69] Michela Talò, Beate Krause, Jürgen Pionteck, Giulia Lanzara, and Walter Lacarbonara.  An updated micromechanical model based on morphological characterization of carbon nanotube nanocomposites.  Composites Part B: Engineering, 115:70–78, 2017.
  * [70] Ioannis E Kavvadias, Konstantinos Tsongas, Kosmas E Bantilas, Maria G Falara, Athanasia K Thomoglou, Fani I Gkountakou, and Anaxagoras Elenas.  Mechanical Characterization of MWCNT-Reinforced Cement Paste: Experimental and Multiscale Computational Investigation.  Materials, 16(15):5379, 2023.
  * [71] Mohammad Kazem Hassanzadeh-Aghdam, Mohammad Javad Mahmoodi, and Reza Ansari.  Creep performance of CNT polymer nanocomposites-An emphasis on viscoelastic interphase and CNT agglomeration.  Composites Part B: Engineering, 168:274–281, 2019.
  * [72] Xingyi Zhu, Yuan Gao, Ziwei Dai, David J Corr, and Surendra P Shah.  Effect of interfacial transition zone on the Young’s modulus of carbon nanofiber reinforced cement concrete.  Cement and Concrete Research, 107:49–63, 2018.
  * [73] M. K. Yeh, N. H. Tai, and J. H. Liu.  Mechanical behavior of phenolic-based composites reinforced with multi-walled carbon nanotubes.  Carbon, 44(1):1–9, 2006.
  * [74] Rashid K Abu Al-Rub, Ahmad I Ashour, and Bryan M Tyson.  On the aspect ratio effect of multi-walled carbon nanotube reinforcements on the mechanical properties of cementitious nanocomposites.  Construction and building materials, 35:647–655, 2012.
  * [75] Mahyar Ramezani, Young Hoon Kim, and Zhihui Sun.  Elastic modulus formulation of cementitious materials incorporating carbon nanotubes: Probabilistic approach.  Construction and Building Materials, 274:122092, 2021.
  * [76] Bohayra Mortazavi, Majid Baniassadi, Julien Bardon, and Said Ahzi.  Modeling of two-phase random composite materials by finite element, mori–tanaka and strong contrast methods.  Composites Part B: Engineering, 45(1):1117–1125, 2013.
  * [77] MK Hassanzadeh-Aghdam.  Evaluating the effective creep properties of graphene-reinforced polymer nanocomposites by a homogenization approach.  Composites Science and Technology, 209:108791, 2021.