Numerical Study of Concrete: A Mesoscopic Scale Simulation Methodology

1. Introduction

Concrete, widely employed in construction due to its cost-effectiveness, strength, and durability, has diverse applications in civil engineering. However, its highly heterogeneous microstructure complicates predictions solely based on experiments. Theoretical studies, employing micromechanics analysis, help deduce macroscopic constitutive behavior. Comprising different phases, concrete's behavior and properties exhibit variations across scales, ranging from the microscopic to the macroscopic scales [1,2].

At the microscopic level, concrete is a composition of cement particles and water, forming the intricate cement paste enveloping aggregates like sand and gravel. The microstructure of this paste is intricate, featuring porosities, micro-cracks, and other defects that can significantly influence its mechanical properties [3]. Zooming in to the mesoscopic level, concrete reveals itself as a material consisting of three crucial phases: the mortar, comprising cement, water, and fine aggregate; the interfacial transition zone (ITZ) found between the coarse aggregate and the mortar; and the coarse aggregates themselves [4,5]. Each of these phases contributes uniquely to the overall behavior of the material. At this level, parameters like shape, size, and distribution of coarse aggregates significantly impact concrete's mechanical behavior. This study aims to create random mesostructured models for concrete, considering aggregate particles' coordinates and sizes, then subjecting them to finite element method (FEM) analysis to understand their effects on macroscopic responses. At the macroscopic level, concrete is simplified as a single-phase material characterized primarily by its compressive strength, modulus of elasticity, and its flexural and tensile behavior. This macroscopic view is essential in the design process, where meticulous control and verification of these properties are imperative to ensure the safety and durability of constructed structures [6]. Engineers and architects rely on a thorough understanding of concrete's behavior at different scales to optimize its use in various applications, meeting the stringent demands of modern construction standards.

When examining concrete at the mesoscopic level, which deals with millimeter-scale dimensions, numerical models become valuable tools. These models enable a detailed analysis of how concrete behaves at the specimen level. The material's traits are intricately influenced by the properties of its key components: coarse aggregate, mortar, and the interface between them. Utilizing numerical models at this mesoscopic scale holds significant importance, offering insights into the mechanical aspects of concrete, such as its fracture toughness and ability to bear loads. This approach provides a nuanced understanding of concrete behavior, contributing to advancements in structural engineering and material science [7].

Several mesoscopic models have been developed, each examining material composition's influence on overall behavior [8,9,10]. Aggregate shapes, sizes, and distributions are essential factors affecting stress distribution and crack initiation. Various approaches, including truss models, lattice models, and finite element mesh simulations, have been employed to study concrete's microstructure. Notably, the choice of aggregate size distribution, often described by grading curves like the Fuller curve, plays a vital role in concrete mix design, affecting workability, mechanical strength, permeability, and durability of the material. In terms of simulating aggregate spatial distribution, different methods like take-and-place, divide-and-fill, stochastic-heuristic algorithms, and distinct element methods have been utilized. Meshing techniques, such as aligned and unaligned approaches, help discretize mesoscopic models for numerical analysis. The interfacial transition zone (ITZ) between aggregates and mortar is recognized as influential in microcrack initiation, prompting varied approaches in mesh generation considering ITZ domains [11].

The main novelties of the work have been in two aspects: concerning the algorithm used and regarding the numerical model. In the first case a Matlab code has been developed based on the take and place method, capable of designing Representative Volume Elements (RVE) in prismatic or cylindrical form, with filling characteristics as required in terms of aggregate percentage and size. Additionally, the Interfacial Transition Zone (ITZ) has been considered, and it is possible to account for the presence of pores, although this option has not been chosen in this instance. The geometric arrangement ensures that aggregates do not touch or overlap, thereby allowing the consideration of the specific properties of the ITZ. With respect to the numerical model, for mortar and ITZ, a nonlinear combined mathematical model of concrete based on the Drucker-Prager material model for compression stresses, and Willam and Warnke for tension stresses, is used to predict both material failure modes: cracking and crushing.

For the contact between coarse aggregates and ITZ, a pure penalty algorithm between deformable surfaces, with asymmetric behavior, is used.

2. Materials and Methods

The fact that concrete is brittle in tension but relatively strong in compression, that the strength of the mortar is higher than that of the concrete itself, and that this strength decreases as the size of the coarse aggregate increases, or that the permeability of high-density coarse aggregate concrete is higher than that of the cement paste, are questions about the behavior of concrete that are answered by the existence of the interfacial transition zone (ITZ) between the large aggregate particles and the hydrated cement paste [12]. Although it is composed of the same elements as the hydrated cement paste, the microstructure, and properties of the interfacial transition zone (ITZ) differ from the bulk hydrated cement paste. Therefore, it is treated as a separate phase in the microstructure of concrete.

Wittmann and his team [9] have delved into key aspects related to the Interfacial Transition Zone (ITZ) and the behavior of concrete at the mesoscopic level. They have crafted both 2D and 3D models, meticulously considering the interplay between mortar and aggregates, along with the geometric and mechanical characteristics inherent to each phase.

The mesoscopic model developed by Wittmann is based on the concept that concrete can be considered a composite material, formed by aggregates dispersed in a matrix of cement paste. In this model, concrete is divided into discrete elements called particles, representing the aggregates, and the cement paste. Each particle is individually modeled, considering its geometric and mechanical properties. For instance, aggregates can be represented as spherical or ellipsoidal particles with specific size and strength properties, while the cement paste can be modeled as a continuous matrix with elastic or viscoelastic properties. Wittmann's mesoscopic model allows for the simulation of the mechanical behavior of concrete by considering the interaction between particles [13,14]. External loads can be applied to study how forces are transmitted through the particles and the paste matrix. This allows the analysis of phenomena such as deformation, strength, fracture, and the long-term behavior of concrete. The model provides a more realistic representation of concrete behavior by considering the distribution and interaction of aggregates and allows the study of the influence of different parameters, such as the water-cement ratio, aggregate size, and material properties, on the mechanical properties of concrete. It is also useful for investigating the impact of defects and heterogeneities on material behavior. However, it is important to note that Wittmann's mesoscopic model is a simplification of reality and requires careful selection of parameters and properties for accurate results, and its application can require significant computational resources.

In this work, a detailed mesoscopic model has been devised to represent concrete, consisting of three distinct phases: mortar, coarse aggregate, and the transition zone (Figure 1).

To carry out a successful numerical approximation, experimental tests were conducted on cubic specimens measuring 50×50×50 mm made of concrete dosed with a water-to-cement ratio (w/c) of 0.5, and with CEM I 52.5 cement. All manufactured specimens incorporated the same number of coarse limestone aggregates corresponding to the particle size fractions 8/10 mm, 6/8 mm, and 4/6 mm to ensure, as much as possible, the best reproducibility among tests and the production of specimens very similar to those simulated numerically. Before the tests, the manufactured specimens were cured for 28 days in a humidity chamber with a relative humidity of 95% and at a temperature of 20°C.

The modeling was carried out in two stages: geometrical models are generated using our own Matlab code, and the numerical model is analyzed using FEM with the Ansys-Workbench program.

2.1. Stage 1: Geometrical Model Generation

To carry out the first stage of modeling, coarse aggregates have been spherical in shape, and three phases have been established: mortar, coarse aggregates, and the interfacial transition zone (ITZ). The effect of aggregate shape has been studied using spherical, ellipsoidal, and polyhedral shapes, while keeping the rest of the parameters fixed, and reaching the conclusion that the type of crack is independent of the shape of aggregates and voids. It was also concluded that the load capacity of concrete containing spherical and ellipsoidal aggregates is between 1% and 3% higher than that which includes polyhedral aggregates, perhaps since non-smooth aggregates have more mortar-aggregate interface elements and because the smooth edges of spherical and ellipsoidal aggregates have a more benign stress distribution that delays the fracture process and increases tensile strength [15].

In principle, the ITZ by subtracting the volume of two previously defined concentric spheres was generated: the inner sphere corresponds to the equivalent diameter of the coarse aggregate, and the outer sphere is designed with a diameter equal to that of the coarse aggregate plus twice the thickness of the ITZ. The shell, obtained by the difference between the larger and smaller spherical volumes, would precisely be this weaker zone (ITZ) of the concrete that is generated around the coarse aggregate particles. The programming in Matlab is carried out using the following input parameters:

• The dimensions of the Representative Volume Element (RVE), which correspond to the dimensions of the specimens tested in the laboratory.
• The diameters of the coarse aggregate particles, represented as spheres.
• The mass percentage of each of the sizes of the coarse aggregate particles considered in relation to the total mass.
• The coarse aggregate-cement ratio.
• The volume fractions of cement, aggregates, and water.
• The density of concrete.

The RVE generated in Matlab can be of different size and shape (cylinder, prism, cube). In the simulations carried out in this work, the coarse aggregate particles have diameters of 4, 6 and 8 mm, constituting 8%, 15% and 7% of the total mass, respectively [16]. These volumes have correspondence with those reflected in the Fuller parabola for the same particle sizes. The programming in Matlab has allowed to generate a random distribution of particles in the RVE so that they resemble that of the real concrete in a statistical sense. The position of each coarse aggregate particle within the (RVE) is done randomly, using the Take and Place method [9,17], arranging the spheres (particles), from the largest to the smallest, according to the mass percentages established in the initial data. The number of particles of each diameter i that is introduced in the RVE, is calculated by the expression:

$${np}_i=\frac{{V.p}_i}{v_i}⌊1-cc\left(\frac{1}{{\rho }_c}-\frac{r_{wc}}{{\rho }_w}\right)⌋,$$

(1)

where V is the volume in m³ of the RVE, ${\mathrm{p}}_{\mathrm{i}}$ is the mass fraction of the coarse aggregate of diameter, i with respect to the total mass, ${\mathrm{v}}_{\mathrm{i}}$ is the aggregate volume in m³, cc is the amount of cement per m³ of concrete, ${\mathrm{r}}_{\mathrm{w}\mathrm{c}}$ is the water-cement ratio, ${\rho }_{\mathrm{c}}$ is the density of cement (3150 kg/ m³), and ${\rho }_{\mathrm{w}}$ is the density of water.

The amount obtained by the previous expression is rounded so that the number of particles is an integer. Next, a number of particles of each class is generated, higher than the one needed to introduce in the RVE, so that it meets the condition of the mass percentage previously indicated. Subsequently, the RVE is filled. A Cartesian coordinate system is used in the concrete sample, so that, when placing a coarse aggregate particle inside the considered concrete volume, its position is defined by the coordinates of its center. To be able to place a particle inside the RVE, two conditions must be met:

• There must be no overlaps between two spheres, and they cannot be tangent to each other either.
• All the particles (spheres) generated must have their entire volume inside the RVE.

A third condition establishes that each particle must be surrounded by a mortar layer with a minimum thickness. This implies that there must be a minimum distance between the boundary of a particle and the faces of the defined RVE, as well as a minimum space between two adjacent particles. This distance is defined by the following expression:

$$d_r=1-\alpha \left(d_c+d\right),$$

(2)

where α is an influence parameter for each sphere, ${\mathrm{d}}_{\mathrm{c}}$ is the diameter of the sphere introduced, and d is the diameter of the next sphere that is intended to introduce.

The algorithm also introduces the parameter called particle radii Hsu’s whose value is 0.5. This value refers to the ability of the coarse aggregate to influence the properties and behavior of the concrete made with them. It depends on factors such as their size, shape, texture, and chemical composition. All particles for which $d_r<0$ are discarded.

The filling of the RVE begins with the particles of the largest diameter until the proportion of the agreed aggregate is obtained. Next, the particles of the successive diameters are placed in descending order of size.

As a test, the filling of many RVEs has been simulated. Different shape and sizes have been used, with different configurations of coarse aggregates, following the specifications established for the different volumetric fractions of each of the sizes of the coarse aggregate particles considered. Figure 2 illustrates three RVEs corresponding to a cube of 50 mm edge, a cylinder of 20 mm diameter and 40 mm height and a prism of 25 × 25 × 50 mm, which contain 345, 35 and 87 particles, respectively.

Once the specimen is filled with all the particles corresponding to the established different sizes, a file with the extensión .inp is generated in Matlab. This file contains the coordinates of the centers and the diameters of each of the coarse aggregate particles (spheres), as well as the dimensions of the Representative Volume Element (RVE). Subsequently, this file is transformed into another file with the extension .igs for export to Ansys. Ansys adapts the RVE generated by Matlab using the take and place method for subsequent application of the finite element method (Figure 3).

2.2. Stage 2: Numerical FEM Model

The second phase begins with importing the .igs file from the Ansys-Workbench FEM software. This file contains the geometry described in the first phase.

Once the file is imported, we proceed to generate the phase corresponding to the Interfacial Transition Zone (ITZ) using a Boolean operation. This operation allows us to differentiate the volumes of each pair of concentric spheres obtained in the previous phase. In this operation, we retain the sphere with the smaller diameter and assign it the properties of the aggregate material. The properties assigned to each of the considered phases were obtained experimentally.

Subsequently, we mesh the model. The meshing is done independently, using tetrahedral finite elements for the aggregates and hexahedral elements for the mortar and ITZ.

Considering the ITZ in the model introduces additional complexity due to coupling and assembly issues at the boundary separating the hexahedral elements used in the mortar phase and the ITZ. Initially, an ITZ thickness of 100 µm was established.

During the meshing process, the number of finite elements (FE) was limited to 1,5·10⁶ to make the calculations feasible. To avoid exceeding these limits, a prismatic Representative Volume Element (RVE) was used, representing one-fourth of a cubic specimen with a 50 mm edge tested in the laboratory. Within this volume, the integrity of all included aggregates (spheres) is ensured (Figure 4).

The model has been modified by remeshing the coarse aggregates using tetrahedral finite elements and the rest of the Representative Volume Element (RVE) using hexahedral finite elements. The properties of the limestone aggregates used were assigned to the first group of elements. For the remaining elements in the model, the mortar properties were assigned, except for those in contact with the aggregates. Through programming in APDL, the properties of the Interfacial Transition Zone (ITZ) were assigned to these finite elements.

3. Finite Element Analysis

For finite element analysis, it is assumed that the material is linearly elastic, with its displacements $u$ being infinitesimal and proportional to the applied force $F$. It is also assumed that the boundary conditions remain unchanged during the application of loads on the set of finite elements. With all of this, the equilibrium equations of the finite element obtained for static analysis are:

$$\left[F\right]=\left[K\right]·\left[u\right],$$

(3)

In practice, these assumptions allow us to state that the geometry of the finite elements remains unchanged during the loading process, hence first-order infinitesimal deformation approximations can be utilized. In real-world scenarios, it may be necessary to account for nonlinear deformations even when they are small and do not exceed the elastic limits of ordinary structural materials.

The Concrete Damage Plasticity (CDP) model used combines the theory of isotropic elastic damage with isotropic plasticity in tension and compression to simulate the inelastic behavior of concrete, as it is a quasi-brittle material. Additionally, this model is suitable for both mass and reinforced concrete and allows for static, cyclic, and dynamic loads. The CDP model considers two main fracture mechanisms: tensile cracking and compression crushing. The formulation of continuous damage in elasticity [18] can be expressed as:

$$\mathsf{\sigma }=(1-d)D^e:(\epsilon -{\epsilon }^i),$$

(4)

where σ is the Cauchy stress tensor, d is the damage variable, $D^e$ is the undamaged elastic stiffness tensor, ε is the total strain tensor, and ${\epsilon }^i$ is the inelastic strain tensor, which accounts for the irreversible nature of concrete cracking.

To evaluate the concrete structural behavior at mesoscale level we have used different finite elements as well as material and contact models.

With respect to the meshing, three different methods have been employed: a body-fitted Cartesian algorithm [19] is utilized for mortar and ITZ components, free meshing with specific sizing parameters is employed for coarse aggregate, and our own code with APDL programming is utilized for ITZ components. In detail, the above geometrical model has been meshed with linear FEM elements using a body-fitted Cartesian algorithm [19], with an element size ranging from 6·10^-4 m for the coarse aggregate to 3·10^-4 m for ITZ. The following element types have been used [20,21]:

• The mortar and ITZ have been modeled using a solid type of tetrahedral finite element named SOLID65. This element is appropriate for simulating the nonlinear structural behavior of concrete. It has eight nodes with three degrees of freedom per node: translations in the X, Y, and Z directions.
• The coarse aggregate is modeled by the SOLID185 element as a homogeneous structural solid defined by eight nodes with three degrees of freedom per node.
• Contact between ITZ and coarse aggregate is modeled by the CONTA174 element. This element has eight nodes and three degrees of freedom per node. The contact element is used to simulate contact and sliding between both materials using a pure penalty algorithm between deformable surfaces.

On the other hand, drawing from other research works aimed to simulate the cracking and crushing of concrete [22,23,24,25,26,27,28], the following constitutive models have been employed by the authors:

• For mortar and ITZ, a nonlinear combined mathematical model of concrete, incorporating the Drucker-Prager material model [29,30,31] for compression stresses and the Willam and Warnke model [32,33,34] for tension stresses, is utilized to predict both material failure modes: cracking and crushing.

The Drucker-Prager yield criterion, used to determine if a material has exceeded its

elastic limit, is:

$$\mathrm{F}\equiv \mathsf{\alpha }{I}_1+\sqrt{J_2}-\mathsf{\kappa }=0,$$

(5)

where $J_2$ and $I_1$ are stress invariants:

$$I_1={\sigma }_x+{\sigma }_y+{\sigma }_z,$$

(6)

$$J_2=\frac{1}{6}\left[{\left({\sigma }_x-{\sigma }_y\right)}^2+{\left({\sigma }_y-{\sigma }_z\right)}^2+{\left({\sigma }_z-{\sigma }_x\right)}^2\right]+{\tau }_{xy}^2+{\tau }_{zx}^2+{\tau }_{yz}^2$$

(7)

α is a parameter related to the internal friction angle φ:

$$\mathsf{\alpha }=\frac{2\sin \phi }{\sqrt{3}\left(3+\sin \phi \right)},$$

(8)

and

$$\mathsf{\kappa }=\frac{6c\cos \phi }{\sqrt{3}\left(3+\sin \phi \right)},$$

(9)

where c is the cohesion angle, for which the following formula is proposed:

$$\mathrm{c}=0.231ln\left(E_0{d}_{max}^2\right)-0.60$$

(10)

based on experience gained from finite element application to concrete.

The Drucker-Prager yield surface is a cone of revolution with its axis along the hydrostatic line (${\sigma }_x={\sigma }_y={\sigma }_z$) (Figure 5). Plastic deformations are given by:

$$\mathrm{d}{\epsilon }_{\mathrm{i}\mathrm{j}}^{\mathrm{p}}=\mathrm{d}\lambda \frac{\partial F}{\partial {\sigma }_{ij}},$$

(11)

where dλ is a scalar that provides the magnitude of plastic deformation. The vector of plastic deformations $\mathrm{d}{\overrightarrow{\epsilon }}^p$ is perpendicular to the yield surface. In this case, the volumetric plastic deformation is:

$$\mathrm{d}{\epsilon }_{vol}^{\mathrm{p}}={d\epsilon }_x^{\mathrm{p}}+{d\epsilon }_y^{\mathrm{p}}+{d\epsilon }_z^{\mathrm{p}}=3\alpha \mathrm{d}\lambda ,$$

(12)

This model can be incorporated into the equations governing the behavior of mortar and ITZ as a yield condition in the nonlinear constitutive model [35].

The failure surface of Willam and Warnke reveals that the compressive strength of concrete subjected to multiaxial compressive stresses can be much greater than its uniaxial compressive strength. When considering ${\sigma }_x\ge {\sigma }_y\ge {\sigma }_z$, the sign of ${\sigma }_z$ determines the mode of failure of concrete. For concrete to tend to fracture due to tension, ${\sigma }_z$ must be positive, and the crack will appear in the direction perpendicular to this principal stress. In cases where no stress is positive, crushing failure will occur. In Ansys, the failure criterion due to multiaxial stress state is:

$$\frac{F}{f_c}-\mathrm{S}\ge 0,$$

(13)

where $F$ is a function of the principal stresses (${\sigma }_x,{\sigma }_y,{\sigma }_z)$, S is the failure surface, expressed as a function of the principal stresses, and $f_c$ is the uniaxial compressive strength, whenever the condition

$$\left|{\sigma }_h\right|\le \sqrt{3}{f}_c,$$

(1)

is met, where

$${\sigma }_h=\frac{{\sigma }_x+{\sigma }_y+{\sigma }_z}{3}$$

(1)

since the remaining concrete strength parameters, Ansys takes them by default. For the complete definition of the material, in addition to the elastic modulus, Poisson's ratio, and the coefficients of transfer of shear stress for open and closed cracks must also be introduced.

With this model, if cracking occurs, Ansys nullifies the elastic modulus of the concrete element in the direction parallel to the principal tensile stress. Conversely, if crushing occurs, Ansys nullifies the elastic modulus of the concrete element in all directions, which is equivalent to removing the element from the model.

• For coarse aggregate, a pure linear elastic model is employed due to its high strength. Besides, regarding the contact between coarse aggregates and ITZ, a pure penalty algorithm between deformable surfaces with asymmetric behaviour is used.
• Furthermore, the relevant mortar material properties employed in the model include:
  Density of concrete, 2150 kg/ m³.
  Young’s modulus, 42000 MPa.
  Poisson’s ratio, 0.2.
  Nonlinear properties for the Drucker-Prager material model:
  Ultimate uniaxial crushing strength, 42 MPa.
  Angle of internal friction, 35°.
  Dilatancy angle, 35°.
• The ITZ properties used in the model are:
  Young’s modulus, 25000 MPa.
  Ultimate uniaxial compressive strength, 37 MPa.
  Angle of internal friction, 35°.
  Dilatancy angle, 19°.
• Additionally, the properties of the limestone aggregates employed in the model are:
  Young’s modulus, 65000 MPa.
  Poisson’s ratio, 0.18.
• The following forces and boundary conditions to the FEM model are applied: a symmetry is applied to ZX and YZ faces, and an imposed displacement of -7·10^-5 m applied to the XY face.
• The numerical model, using a Newton-Raphson integration scheme with a force convergence parameter of a 2.5 with respect to the L2 Euclidean norm, with a minimum value of 0.01 N, is solved. The applied displacement, for each sub-step from the value of the previous load step ranging from 1% to a maximum value of 2%, is linearly interpolated.

4. Numerical and Experimental Results

Based on the above stages, four different numerical and experimental tests have been done. The main results are exposed in the following sections.

4.1. Experimental Results

In Figure 6, the setup for one of the experimental tests carried out can be observed. Two strain gauges attached to two opposite faces of the specimen, to determine the stress-strain curve of the material, are used.

In Figure 7, stress-strain curves of four experimental tests are shown. The first test shows an initial nearly linear phase that gradually curves as deformations increase. It achieves an ultimate tensile strength of 45,9 MPa with a strain close to 1,8 per thousand. In the second test, a slight deviation in the curve occurs when deformations exceed 8 per thousand due to a temporary malfunction in one of the strain gauges. Nonetheless, the test proceeds normally until failure, reaching deformation levels like the previous test and a tension slightly above (49,9 MPa). The third experimental test yields an ultimate tensile strength of 57,1 MPa at a strain close to 1.6 per thousand, with a less pronounced curvature compared to the previous tests. Finally, in the fourth test, a stress-strain curve characteristic of the material is observed, culminating in failure at a tension level akin to the first two tests (48,7 MPa) and a strain close to 1,7 per thousand.

4.2. Numerical FEM Results

The problem was solved in a workstation computer with a two Xeon 64-bit CPU E5630, 64 GB RAM memory, 6.0 TB hard disk and eight cores. The total elapsed average CPU time per nonlinear structural analysis was about 70,000 s and the total number of iterations to get convergence about 500. The Newton-Raphson method was used to solve the nonlinear mathematical model, both due to the constitutive model of the material and the contacts between the aggregates and the ITZ. The differences obtained between the values of "Force convergence" and "Force criterion" are very small, ranging from 1,4·10^-7 to 31,7 N. Additionally, a smooth approximation of the iterations performed by the method until reaching the solution, indicating that the problem has been well-conditioned.

Figure 8 illustrates the curves derived from the four numerical simulations. Notably, all simulations exhibit an initial phase that is nearly linear, with almost the same slope in all of them. In simulations 1 and 4, the curve deviates at a stress level around 40 MPa, ac-companied by deformations of 1 per thousand. Conversely, simulations 2 and 3 depict cracking initiation at stress levels very close to 37,5 MPa, with deformations in both cases around 0,95 per thousand. Interestingly, the ultimate tensile strength in the first and fourth tests converges to similar failure values of 43 MPa, corresponding to strains of 1,2 per thousand and 1,3 per thousand, respectively. These discrepancies in ultimate tensile strength and strain among simulations are due to the different distributions of coarse aggregates in each of the geometric models employed to simulate the specimens, as generated using Matlab.

Cracking and crushing failures for different zones analyzed are shown in Figure 9 by means of elements in red color. Overall view, coarse aggregate, central section of the specimen and ITZ zone are showed in the same figure. It can be observed in the Figure 9 that failure occurs in the weakest part of the concrete, namely, in the ITZ (Interfacial Transition Zone).

In Figure 10, directional deformation along the z-axis from the fourth part of the prismatic specimen can be observed. Figure 11 depicts the total deformation experienced by the fourth part of a prismatic RVE.

4.3. Numerical and Experimental Comparison

Figure 12 illustrates a strong agreement between the experimental findings and those derived from numerical simulation. In the latter case, the compressive strength measures 41,34 MPa with a strain of 1,04·10-3 m. Notably, the stress-strain curve from the laboratory tests aligns closely with this simulation data. Here is shown an advanced numerical model that provides reliable results based on four different FEM models, and which is intended to be used in the future to address more complex problems that require more effort, both experimental and computational.

5. Conclusions

The proposed model forms the groundwork for establishing a methodology aimed at the numerical simulation of compression tests, a commonplace in laboratory settings. Several key points arise in the discussion:

• Detailed Insights into Mesoscopic Behavior: The model promises to deliver comprehensive insights into the mesoscopic-scale behavior of concrete. This detailed information is pivotal for unraveling deformation and fracture mechanisms, ultimately contributing to more informed and effective material design.

Cost-Effective and Rapid Evaluation: Mesoscopic numerical simulations present a cost-effective and time-efficient alternative for assessing concrete behavior. While they offer notable advantages, it is essential to acknowledge that they cannot entirely replace experimental methods. Both approaches complement each other, each with its unique strengths and limitations.

• Identification of Failure-Prone Areas: Through these simulations, it becomes possible to identify areas within the material that are more susceptible to failure. This insight is invaluable for targeted improvements and reinforcement in structural design. The aim of this work is to compare numerical results with those obtained in the laboratory, hence we have not compared our numerical results with others obtained in a similar way, or with theoretical methods. In summary, the proposed mesoscopic model not only opens avenues for efficient simulation of compression tests but also holds the potential to deepen our understanding of concrete behavior, paving the way for advancements in material science and engineering.
• The aim of the authors is to employ soon this methodology as a substitute for extensive experimental testing campaigns, which generate significant waste, by utilizing numerical simulations that yield comparable results.