Next Article in Journal
Effect of Elevated Temperature on the Bond Strength of Prestressing Reinforcement in UHPC
Next Article in Special Issue
Effect of Hydrophobic Treatments on Improving the Salt Frost Resistance of Concrete
Previous Article in Journal
Shear Wave Splitting and Polarization in Anisotropic Fluid-Infiltrating Porous Media: A Numerical Study
Previous Article in Special Issue
Permeabilities and Mechanical Properties of Hardened Cement Pastes Modified with Sodium Laurate and Nano Silica
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Cementitious Composites with High Compaction Potential: Modeling and Calibration

by
Giao Vu
1,
Tagir Iskhakov
1,
Jithender J. Timothy
1,
Christoph Schulte-Schrepping
2,
Rolf Breitenbücher
2 and
Günther Meschke
1,*
1
Institute for Structural Mechanics, Ruhr University Bochum, Universitaetsstrasse 150, 44801 Bochum, Germany
2
Institute for Building Materials, Ruhr University Bochum, Universitaetsstrasse 150, 44801 Bochum, Germany
*
Author to whom correspondence should be addressed.
Materials 2020, 13(21), 4989; https://doi.org/10.3390/ma13214989
Submission received: 15 October 2020 / Revised: 30 October 2020 / Accepted: 1 November 2020 / Published: 5 November 2020
(This article belongs to the Special Issue Concrete and Construction Materials)

Abstract

:
There is an increasing need for the development of novel technologies for tunnel construction in difficult geological conditions to protect segmental linings from unexpected large deformations. In the context of mechanized tunneling, one method to increase the damage tolerance of tunnel linings in such conditions is the integration of a compressible two-component grout for the annular gap between the segmental linings and the deformable ground. In this regard, expanded polystyrene (EPS) lightweight concrete/mortar has received increasing interest as a potential “candidate material” for the aforementioned application. In particular, the behavior of the EPS lightweight composites can be customized by modifying their pore structure to accommodate deformations due to specific geological conditions such as squeezing rocks. To this end, novel compressible cementitious EPS-based composite materials with high compaction potential have been developed. Specimens prepared from these composites have been subjected to compressive loads with and without lateral confinement. Based on these experimental data a computational model based on the Discrete Element Method (DEM) has been calibrated and validated. The proposed calibration procedure allows for modeling and prognosis of a wide variety of composite materials with a high compaction potential. The calibration procedure is characterized by the identification of physically quantifiable parameters and the use of phenomenological submodels. Model prognoses show excellent agreement with new experimental measurements that were not incorporated in the calibration procedure.

1. Introduction

Difficult geological conditions such as “squeezing rock” [1,2] can induce severe damage to tunnel construction. In general, the time-variant soil deformations increase continuously after completion of the tunnel structure and may severely affect the long term integrity and safety, leading eventually to the complete loss of operability of the tunnel [3]. One possible approach to avoid this substantial risk of damage is the incorporation of lining materials that are highly compressible [4,5,6,7]. These compressible lining materials accommodate large deformations after a certain threshold level of stresses is reached. By allowing the ground to deform, the squeezing pressure acting on the regular concrete tunnel linings will be reduced. Thus, highly compressible lining materials can serve as cushions that protect the regular concrete linings from deterioration. A compressible segmental lining system can be realized either by arranging one or more compressible layers between the lining and the ground (radial compressible systems) or by introducing deformable elements between regular concrete linings (longitudinal/circumferential compressible systems), see in [8] for a discussion on the comparison of various compressible segmental lining systems. Radially compressible systems are in general comprised of a regular rigid cementitious segmental lining, optionally a compressible layer that is attached to this lining during pre-fabrication and a compressible grout that is used to fill the annular gap. A compressible grouting material, in addition to being highly compressible, must satisfy the usual requirements of gap grouting materials. Various compressible materials such as expanded clay [9] to cement based composites such as light weight concrete [10], compressible cementitious mortar with expanded polystyrene [11] or expanded pearls, and foam [12] have been proposed. In general, the basic morphology of a cementitious compressible grout consists of pores and soft inclusions embedded in a cementitious matrix. A laboratory specimen, made of a compressible cementitious composite, that is subject to compressive loads, undergoes three characteristic stages of deformation: (i) an elastic stage with reversible deformations, (ii) a plateau stage where the material undergoes large irreversible deformations, and (iii) a densification stage in which the material has exhausted the compaction capacity and stiffens with increasing loads. At the scale of the microstructure, the material undergoes several complex processes such as reversible topology preserving deformations, gradual collapse of pore/soft-inclusion and compaction of the binder particles.
Designing a compressible grout material with the desired compaction behavior, for a certain specific geological condition, can be challenging as multiple complex phenomena govern the compaction behavior of the material. In order to support the design of novel compressible lining materials and also gain a deeper understanding of the underlying mechanisms, computational simulations can be applied. As these materials are characterized by a highly heterogeneous microstructure with pores and/or soft inclusions [13,14], the type of computational model depends on the desired scale at which the corresponding physical mechanisms are specified: macro-, micro-, or mesoscale. In contrast to macrosopic phenomenological models [15,16,17,18] where the microstructure is treated as a homogeneous medium, continuum micromechanics models take into account the influence of the microstructure using the matrix-inclusion concept [19] in conjunction with inelastic material laws defined at the microscale (elastoplasic [20,21,22] or hyperelastic [23,24,25]) and the utilization of nonlinear homogenization schemes [26,27]. However, as the microstructure is not explicitly resolved but idealized in terms of a representative effective medium, deformation gradients at the scale of microstructure and formation of localization bands cannot be simulated. In contrast, mesoscale models explicitly resolve the microstructure and the interactions between the heterogeneities such as the contact of pore/microcrack faces during compaction. Mesoscale models can be developed using a variety of methods such as the Finite Element Method (FEM) [28,29], the Discrete Element Method (DEM) [30,31], and lattice methods [32,33]. The Discrete Element Method is based on an idealization of the mesostructure using discrete particles, which interact through inter-particle forces and contact mechanisms. This discretization method allows for an explicit representation of the material microstructure with soft [34,35,36] and hard inclusions [37,38,39,40]. Thus, through an explicit representation of the pore structure and by consideration of the pore collapse mechanism, DEM is capable of simulating the essential physics of compressible lining materials. Unfortunately, the calibration of the parameters governing the inter-particle interactions in DEM is generally not easily accomplished directly from laboratory tests, but usually requires inverse identification procedures, see, e.g., in [40,41].

Goals and Structure of the Paper

In this paper, the behavior of several designs of highly compressible cementitious material composites used to fill the annular gap grouts in mechanized tunneling is investigated using experimental methods. The material behavior of these materials is thereafter analyzed numerically using the Discrete Element Method. Emphasize is laid on the calibration procedure of the inter-particle parameters, which are identified using the data from the experimental procedure. The DEM model is not only intended as a tool for the prognosis of the material behavior, but also to obtain insights into the physical phenomena characterizing the interplay between the changes in the microstructure and the macroscopic deformation behavior.
The remainder of this paper is structured as follows. In Section 2, the theoretical background of DEM and the DEM constitutive contact model for cement paste are reviewed. Section 3 summarizes the experimental data, and discussed the calibration and the validation process. Finally, a discussion of the main findings and concluding remarks are provided in Section 4.

2. Modeling Cementitious Materials with DEM

In Discrete Element Method (DEM) models, the material is described as an assembly of particles that can collide, interact, and exert forces on each other. The dynamics of these particles is governed by Newton’s second law. Within this modeling framework, concrete and rock materials are characterized by a packing of particles linked together by cohesive frictional forces. Within the medium, the induced forces are transmitted via a contact network between particles [30].
The contact network is first established and updated by identifying the particles and their nearest neighbor interactions. Thus, particles are considered to interact with each other, if the following condition is fulfilled,
l i j R I ( r i + r j ) ,
where r i and r j are the radii of particles i and j, and l i j denotes the distance between their centers. R I is the interaction factor set equal to 1 for granular materials. However, for cohesive frictional materials such as concrete, R I is often set to a higher value, e.g., to 1.5 to increase the number of average cohesive contacts per particle, which represents the cohesive properties of the concrete matrix [36,42]. Most DEM formulations are based on the soft-sphere approach, where contact is characterized by an interaction between overlapping particles. However, in the hard-sphere approach [43], contact occurs without allowing an overlap between rigid particles and is instantaneous. An instantaneous point-contact event between rigid spheres in the hard-sphere approach is rather simplistic and this method cannot account for multiple simultaneous contacts between a large number of particles as well as inelastic interactions between particles. A detailed explanation of these two approaches can be found in [44]. In this work, the soft-sphere approach is used.
The interaction forces are evaluated based on the relative displacements in the current particle configuration. Next, the resultant interaction forces are used together with the applied external forces as input for the equations of motion in the time integration step to solve for the new position of all particles.

2.1. Governing Equations of Motion

For every particle i, the resultant force F i is the sum of an external force F i e x t , a damping force F d a m p and a contact force F i j , where j defines particles which are in contact with particle i:
F i = F i e x t + j = 1 N F i j + F d a m p ,
where N is the number of particles having an interaction with particle i.
Numerical damping F d a m p [36,42] is applied to all particles. The damping term dissipates the overall kinetic energy and ensures quasi-static equilibrium conditions. Given the resultant force F i on particle i having mass m, the velocity v i is evaluated using Newton’s second law as
v i n + 1 2 = v i n 1 2 + F i n m i Δ t ,
where Δ t is the time increment.
Similarly, the effective moment M i and the angular velocity ω i are updated at each time step:
M i = M i e x t + j = 1 n j M i j I i + M d a m p ,
ω i n + 1 2 = ω i n 1 2 + ω ˙ n Δ t ,
where I i is the moment of inertia of particle i.
The displacement u i n + 1 and orientation θ i n + 1 are updated according to
u i n + 1 = u n + v i n + 1 2 Δ t ,
θ i n + 1 = θ i n + ω i n + 1 2 Δ t .
As a result, the position of all particles is updated accordingly. In the next time step, the newly updated configuration is used to resolve the interaction between particles in terms of contact stresses and strains.

2.2. Constitutive Law

Given the updated position of all particles, the interaction forces in the contact network are computed. The normal and tangential strains ε n , ε τ are computed based on the relative normal and tangential displacements u n , u τ between two particles and their initial distance l 0 i j :
ε n = u n l 0 i j , ε τ = u τ l 0 i j .
The normal relative displacement vector is calculated as
u n = ( x i x j l 0 i j ) n ,
where x i , x j denote the current positions of the particle centroids and n is the normal vector connecting the centroids of two particles. The relative tangential displacement is computed by subtracting the normal component from the total relative displacement u i j
u τ i j = u i j ( u i j · n ) n ,
where
u i j = ( Δ u i Δ u j ) + ( ω i × r c i ω j × r c j ) Δ t .
Δ u i and Δ u j are the displacement increments of particles i and j, and r c i and r c j denote the vectors connecting the point of the contact and centroid of particles.
The elastic behavior between two particles is characterized by the normal K n and tangential K τ stiffness moduli. The normal stress and shear stress are computed directly from the updated position as
σ n = K n ε n , σ τ = K τ ε τ .
In this work, the concrete contact modeling approach according to the work in [45] was adopted, where the interaction of the particles tension in normal direction is governed by a damage softening law
σ n = [ 1 ω H ( ε n e ) ] K n ε n e .
where ω is a damage parameter ( ω [ 0 , 1 ] ) and H ( ϵ n ) is the Heaviside function, used to deactivate damage in compression. The linear softening law is characterized by the predefined limit elastic strain ε 0 and the ultimate strain ε f , as shown in Figure 1.
The normal compression mode is characterized by an elasto-plastic behavior with plastic strain ε s and the relative hardening stiffness K s (Figure 1)
ε n = ε n + σ n K n , if σ n < σ s , ε s + σ n K n ε s K s K n , o t h e r w i s e .
Shear stress can be calculated from the modified Mohr–Coulomb frictional law f ( σ n , σ τ ) , taking into account the damage parameter ω (see Figure 2):
f ( σ n , σ τ ) = | | σ τ | | [ c 0 ( 1 ω ) σ n tan φ ] ,
where c 0 is the initial cohesion and φ is the frictional angle. Shear stress and strain are updated as follows. First, the trial shear stress is computed
σ τ t r i a l = K τ ε τ ,
then the shear stress is corrected according to
σ τ = c 0 ( 1 ω ) σ n tan φ σ τ t r i a l | σ τ t r i a l | .
Finally, the shear strain is recomputed as
ε τ = σ τ | σ τ t r i a l | ε τ | ε τ | .
Given the contact stresses, the contact forces are obtained as
F n i j = σ n i j A i j n ,
F τ i j = σ τ i j A i j ,
where A i j is the contact interface area [36] which is defined as A i j = π m i n ( r i , r j ) 2 .
All parameters of the described constitutive model summarized in Table 1 are to be calibrated according to the properties of the specific material to be analyzed. The calibration procedure for cementitious compressible composites is given in the following Section.

3. Calibration and Validation of the DEM Model Based on Laboratory Experiments

3.1. Experimental Data from Compression Tests on Highly Compressible Composite Grouts

As possible candidates for annular gap grout, three different activated two component-grout mixes with a defined fraction of expanded polystyrene (EPS) beads (denoted as A, B, and C) have been prepared within the experimental program. The type of binder is the same for all mixes. Mix A consists of only cement binder matrix (porosity < 1%); mix B consists of cement paste and EPS beads with a volume fraction of 61.5%; and mix C consists of cement paste, EPS, and additional air voids, contributing up to 71.1% of volume fraction in total (see Figure 3 and Table 2). Two cylindrical samples with diameter of 100 mm and height of 200 mm (see Figure 4 left) of each mix have been casted and cured for 7 days.
The samples were subjected to both unconfined uniaxial and confined uniaxial compression tests with a constant loading rate of 20 mm/min (or strain rate ε ˙ = 0.1 s−1). First, a uniaxial compression test (see Figure 4 left) was performed on sample A to obtain the mechanical properties of the cement paste. Then, uniaxial compression and confined compression (see Figure 4 right) tests were conducted on samples B and C to investigate the effect of EPS and air voids on the mechanical behavior of the grout. To simulate confinement, the sample was placed in a steel container (diameter = 180 mm) filled with fine sand (see Figure 4 right) in the confined tests. All tests were performed using displacement control.
Table 3 presents elasticity modulus, compressive strength and density of samples A, B, and C as the outcome of the uniaxial compression tests (Figure 5 left). Young’s moduli of samples A, B, and C are estimated using the stress and strain intervals indicated by black lines in Figure 5 left. Grout samples B and C, owing to the high void content, exhibited a decrease in strength and stiffness as compared to sample A. In contrast to sample A, which failed due to tensile splitting, in samples B and C damage initiated at the upper part where the load was applied, followed by local crushing.
In the confined compression test, the behavior of samples B and C differs from sample A. Three distinct stages of deformation were observed in samples B and C under confined compression: An elastic region, a plateau and a densification stage (Figure 5 middle). After the elastic stage, in the plateau stage, the material undergoes plastic deformation characterized by large plastic deformations with marginal increase in stresses, associated with the collapse of voids. When the pores have completely collapsed, subsequent to the plateau stage, densification due to pore compaction results in a significant increase in the tangent stiffness of the material. The confinement condition has prevented the samples from failure resulting from lateral expansion. In the confined compression tests, irreversible compaction of the composites up to 75% was attained (Figure 5 right)).

3.2. Calibration of Model Parameters

Using the data obtained from the experiments, we proceed to calibrate the model parameters as follows. (a) Data obtained from uniaxial compression of sample A is used to calibrate the contact parameters for the cement paste matrix; (b) data obtained from uniaxial compression tests on samples B and C is used to calibrate the porosity of DEM numerical models for the samples B and C; (c) the stress–strain curve of B under confined compression is further used for calibrating additional DEM parameters associated with the compression behavior of the material; (d) finally, the data from the confined compression experiments using sample C is exclusively used for the validation of the model.
The calibration of parameters for the described model is performed following three steps according to the procedure illustrated in Figure 6.

3.2.1. Calibration Step 1

First, a computational model of a cylindrical sample of height 200 mm and diameter 100 mm has been generated for the DEM simulation as shown in Figure 8 left. In DEM, a sample is an assembly of spherical particles occupying a given geometry and is often referred to as packing. Here, a dense packing with a random arrangement was chosen. DEM particles of the same size (radius = 0.8 mm) were used. This packing is assumed to represent the cement paste matrix without EPS or air voids (sample A), which results in a total of 452,257 particles. However, these DEM particles do not represent the actual morphology of the cement paste, rather being a means of material discretization at the mesoscale. In order to resolve the actual geometrical morphology of the cement paste, including the complete pore space ranging from a few nanometers up to a few micrometers, would require a tremendous amount of computational power. For cementitious materials, the interaction factor R I is chosen to be 1.5 [42,46].
Next, the uniaxial compression simulation was performed with this numerical sample analogous to the experimental test. In the model, particles at the bottom cylinder face were fixed in all directions while at the top face a vertical constant strain rate was applied. In order to simulate quasi-static conditions, numerical damping with a damping factor of 0.3 was adopted [45,47] to dissipate the total kinetic energy of the system. Moreover, mass scaling of (4800 kg/m−3) was adopted to increase the critical time step [45,47], which is a standard value for modeling concrete using the DEM. This allows to apply strain rates that are is not too small to enable reasonable computational costs of the analyses. Furthermore, a constant loading velocity of 5 × 10 2 m/s was set in all simulations to exclude the effect of inertia. During the simulations, the resultant forces and displacements at the top cylinder face were recorded and used to compute the stress and the deformation.
As a result, by matching the experimentally measured elastic properties and the compresssive strength of sample A (see Figure 7 left) the normal and tangential elasticity moduli ( K n , K τ ), the parameters of the damage law ( ε 0 , ε f ) and the Mohr–Coulomb yield surface parameters ( c 0 , t a n φ ) are calibrated (Table 4). It is noted that the early portion of the load displacement curve obtained from the experiments in Figure 7 left shows a nonlinear behavior, which is attributed to the loading discrepancy between the specimen and loading plate in the initial stage of loading, a phenomenon commonly observed in laboratory testing [48]. This initial disturbance has been filtered out and does not affect the calibration procedure.
To visualize the damage pattern, we define the parameter ω ^ i , which is a damage parameter averaged over all cohesive contacts (i.e., contacts that are created at the initiation step) associated with particle i. Accordingly, ω ^ i of particle i is defined as
ω ^ i = j = 1 n j ω i j n j ,
where ω i j is the damage parameter of the normal interaction between particle i and particle j, and n j is the number of initial cohesive bonds of particle i. Particles with ω ^ i = 0 are in an undamaged state, while ω ^ i = 1 indicates a fully damaged state.
Figure 7 illustrates the damage pattern in sample A in the experiment (middle) and in two stages of uniaxial compression in the numerical simulation (right). The figures show a cross section through the cylindrical sample. In the simulation, damage first occurs by the formation of a macrocrack near the outer surface, almost parallel to the loading axis, followed by a secondary, inclined macrocrack. The photo from the damaged specimen after the end of the test also shows dominant vertical cracks at both edges of the specimen, partially with a slight inclination as was observed in the results from the DEM simulations.

3.2.2. Calibration Step 2

Ideally, the microstructure of the material used in the experiment should serve as the direct input to generate the numerical sample for DEM simulation, e.g., by incorporating data from CT imaging techniques in DEM models with a realistic microstructure (i.e., cement paste, aggregates, ITZ, etc.) [38]. However, in this paper, the composite microstructure of the investigated grout mixes is characterized by inclusions (pores, air voids, and EPS beads) at multiple scales. Thus, it would require a tremendous amount of very small DEM particles to resolve the topology of the small inclusions. This would be computationally too expensive. Therefore, the numerical models for the samples B and C have been created without CT scans by “embedding” spherical voids into the numerical DEM model for sample A generated and calibrated in Step 1. The voids are embedded into the DEM specimen by removing the DEM particles lying within the spherical region defined by the void location. These spherical voids represent both air voids and EPS beads, as the stiffness of EPS is comparatively low. Thus, from modeling perspective we do not distinguish between air voids and EPS beads.
Using the DEM model with embedded voids, a uniaxial compression test was simulated (same as in Step 1), and the compressive strength in uniaxial test was recorded. Voids were embedded in an iterative manner until the compressive strength of the sample obtained in simulation matches the experimentally measured compressive strength of the corresponding mixes B and C. It must be noted that, DEM particles do not represent the actual cement paste grains, and that at this level of observation, they are rather a means to discretize the material at the mesoscale. Consequently, certain parameters controlling, e.g., plasticity in shear and in normal compression, have to be employed to capture the correct physics of interactions at the lower scale.
The process of embedding the voids into the DEM sample is as follows.
  • Voids with a prescribed volume fraction and size distribution are randomly picked and placed within a cylindrical domain (height = 200 mm, diameter = 100 mm).
  • The coordinate and radius of each void particle is recorded.
  • Given a dense packing of DEM particles generated in Step 1, the DEM particles lying within the spherical region defined by void position and radius are removed.
The void size distribution was assumed to follow the normal probability density function with a mean radius of 2.5 mm and standard deviation of 0.25 (see Figure 8 right). The numerical packing based on this void size distribution can be qualitatively compared with the void distribution along the cut surface of the grout sample prepared in the experimental program (see Figure 8 right).
The geometrical data representing the air voids embedded in the numerical models for samples B and C are summarized in Table 5.
Figure 9 right shows the distribution of fracture, represented by the damage parameter ω ^ , according to the DEM simulations of the samples corresponding to mixes B and C at at strain level 0.008. ω ^ is the damage parameter averaged over all contacts of a corresponding particle. At the top, the damage distribution at the outer surface is illustrated, and at the bottom, the damage distribution along the cross section of the specimen is illustrated.
In comparison to the experimental data, uniaxial compression simulations on the generated samples with voids yielded higher stiffness for the same strength level (see Figure 9 left). This is due to the relatively coarse discretization which does not capture small voids present in the real material. Initially, prior to reaching the peak stress, microcracking is predicted to initiate diffusively at different locations in the specimen, mostly localized in the vicinity of voids. As the load is increased, damage starts to localize at the upper part of the specimen, followed by local crushing/compaction (see Figure 9 right, which shows the state of damage at the axial strain level ε = 0.008 in the post-peak state).

3.2.3. Calibration Step 3

In the previous two calibration steps, the numerical DEM models for samples representing mixes B and C have been generated and calibrated based on the uniaxial compression experiments performed with mixes B and C. In Step 3, the simulation of confined compression on numerical sample B is performed and compared to experimental data in order to calibrate the plasticity parameters K s , ε s (Table 4). Figure 10 left shows the comparison of the stress-strain response for mix B obtained from the calibrated DEM model (red lines) and in the laboratory test (green dotted lines).

3.3. Validation

In the calibration procedure, the normal and tangential elasticity moduli ( K n , K τ ), the softening curve parameters ( ε 0 , ε f ) and parameters defining the Mohr–Coulomb yield surface ( c 0 , t a n φ ) have been calibrated in Step 1 based on the uniaxial compression experiments with mix A. In Step 2, the microstructure topology of the composites corresponding to mixes B and C has been characterized by matching the compressive strength obtained numerically to the measured data using uniaxial compression tests on samples B and C. In Step 3, the plasticity parameters K s , ε s have been calibrated based on the experimentally observed material response of sample B in the confined compression test.
After calibrating the required parameters, the model is now used to predict the complete behavior of sample C under confined compression until a strain level of 0.6. The uniaxial unconfined compression test was used to calibrate the volume fraction of embedded voids for sample C. The experimental results from sample C under confined compression were not used to calibrate any model parameter. This ensures the predictive capability of the model. Figure 10 left presents the predicted (blue lines) and experimentally observed stress–strain response (yellow lines) for the sample C under confined compression conditions.
In the DEM simulation, initially diffuse cracking occurs arbitrarily within the whole specimen. This plateau stage is characterized by material compaction and a slight stress increase. According to the computational model, damage initiates at the top part of the sample and forms a compression band which propagates downwards during the pore compaction process. This leads to a material compaction gradient which can be clearly seen in Figure 10 right top.
During the pore compaction process, most pores collapse “layer by layer”. Particles with all contacts “broken” dynamically interact with each other by means of newly created interactions (i.e., interaction between two particles as they “collide” with each other) and rearrange themselves to fill the voids (Figure 10 right bottom). This mechanism is reflected by a high material compaction without significant change in stress leading to the plateau behavior of the stress–strain curve also observed in the experiments. As soon as all voids have experienced collapse, a densification process initiates, which is characterized by the regain in stiffness. Moreover, this is in agreement with the observations in the laboratory.
In this work, all simulations were performed using an open-source DEM software ‘WooDEM’ [49]. The software is written in C++11 and supports OpenMP parallelization. The computational time required for the simulation of the confined compression test on sample B (Figure 10 left) is 207 h for 2.4 × 10 6 time steps. Each simulation is performed on 12 Intel® Xeon® Gold 6148 CPUs running in parallel @ 2.40 GHz with 100 GB RAM.

4. Conclusions

In this paper, we have generated a computational model employing the Discrete Element Model to simulate the behavior of specimens made of highly compressible cementitious materials. To enable a systematic calibration procedure, two different concrete mixes as well as a specimen made of cement paste only have been tested experimentally. One of the mixes contained EPS beads of different sizes, and the second mix also contained, in addition, artificial air voids. Cylindrical samples based on selected material designs have been subjected to compressive loads with and without lateral confinement. The calibration procedure used information from the test on cement paste specimens as well as the compressive strength data from the highly compressible specimens with only EPS beads subjected to uniaxial compression. The DEM model was eventually validated based on test results from a highly compressible specimen with a material mix also containing air voids subjected to confined compression. Using this calibration procedure, the main physical mechanisms associated with compaction processes occurring in the specimens could be well captured by means of the DEM model. The model has shown a cascade-type mechanism of pore collapse, which propagates from the top to the bottom. This mechanism was also observed in the experiments. After the compaction process was completed, the model showed, again in agreement with the experiments, a re-stiffening of the material. From the computational and experimental analysis, the following conclusions are drawn.
  • Cementitious materials with high-compaction potential can be designed using a combination of weak inclusions and pores.
  • Experimental observations and model simulations show the development of compaction gradients during confined uniaxial compression tests.
  • Despite the extensive work dedicated to the calibration procedure as well as the high computational cost, DEM has shown its capability to replicate the main physical mechanisms governing the behavior of compressible cementitious composites.
  • In order to capture the effect of fine pores with a characteristic size smaller than the DEM discretization, a phenomenological plasticity-type submodel has to be calibrated in addition to the usual inter-particle parameters.
  • The proposed calibration procedure offers a good control of the pore structure characteristics, such as void volume fraction, air-void size, and void size distribution. Consequently, the proposed computational model allows to support the design of new materials with specific, customized compaction properties (elastic phase, plateau, and densification). These materials can be used for optimizing the compressibility characteristics of annular gap grouts used to fill the tail void gap in mechanized tunneling in case of tunneling projects in rocks with a high squeezing potential.
As an outlook, other materials for inclusions enabling a controlled compaction behavior of cementitious materials will be considered, which are characterized by a crushing mechanism when subjected to large compressive stresses. One candidate for such a damage tolerant composite material is based on using expanded glass beads as inclusion. This concept is currently explored in laboratory tests.

Author Contributions

Conceptualization, All; methodology, All ; software, G.V.; validation, G.V., J.J.T., T.I. and G.M.; formal analysis, G.V. and C.S.-S.; investigation, G.V. and C.S.-S.; resources, G.M. and R.B.; data curation, G.V. and C.S.-S.; writing–original draft preparation, G.V.; writing—review and editing, J.J.T., T.I. and G.M.; visualization, G.V.; supervision, T.I.; project administration, J.J.T. and G.M.; funding acquisition, J.J.T., G.M. and R.B. All authors have read and agreed to the published version of the manuscript.

Funding

This work has been partially supported by the German Research Foundation (DFG) in the framework of Subproject B2 and B3 of the Collaborative Research Project SFB 837 (project number 77309832) and Subproject RUB1 of Research Unit FOR 2825 (project number 398216472). This support is gratefully acknowledged.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript.
EPSExpanded polysterene
REVRepresentative elementary volume
FEMFinite Element Method
DEMDiscrete Element Method

References

  1. Anagnostou, G.; Pimentel, E.; Serafeimidis, K. Swelling of sulphatic claystones–some fundamental questions and their practical relevance. Geomech. Tunn. 2010, 3, 567–572. [Google Scholar] [CrossRef]
  2. Steiner, W. Tunnelling in squeezing rocks: Case histories. Rock Mech. Rock Eng. 1996, 29, 211–246. [Google Scholar] [CrossRef]
  3. Einstein, H. Tunnelling in difficult ground-swelling behaviour and identification of swelling rocks. Rock Mech. Rock Eng. 1996, 29, 113–124. [Google Scholar] [CrossRef]
  4. Mezger, F.; Ramoni, M.; Anagnostou, G. Some concepts for segmental linings in squeezing rock. In Rapid Excavation and Tunneling Conference 2015 Proceedings; Society for Mining, Metallurgy, and Exploration Inc. (SME): New Orleans, LA, USA, 2015; pp. 646–658. [Google Scholar]
  5. Cantieni, L.; Anagnostou, G. The interaction between yielding supports and squeezing ground. Tunn. Undergr. Space Technol. 2009, 24, 309–322. [Google Scholar] [CrossRef]
  6. Nehdi, M.; Khan, A.; Lo, K. Development of deformable protective system for underground infrastructure using cellular grouts. Mater. J. 2002, 99, 490–498. [Google Scholar]
  7. Schneider, E.; Spiegl, M. Convergency compatible support systems. Tunnels Tunn. Int. 2008, JUN, 40–43. [Google Scholar]
  8. Mezger, F.; Ramoni, M.; Anagnostou, G. Options for deformable segmental lining systems for tunnelling in squeezing rock. Tunn. Undergr. Space Technol. 2018, 76, 64–75. [Google Scholar] [CrossRef]
  9. Cucino, P.; Eccher, G.; Castellanza, R.; Parpajola, A.; Di Prisco, C.; SpA, L. Expanded clay in deep mechanised tunnel boring. In Proceedings of the ITA-AITES World Tunnel Congress Bangkok, Bangkok, Thailand, 21–23 May 2012. [Google Scholar]
  10. Strohhäusl, S. TBM Tunnelling Under High Overburden With Yielding Segmental Linings Eureka Project Eu 1079. In Tunnel Boring Machines: Trends in Design and Construction of Mechanical Tunnelling; CRC Press: London, UK, 1996; pp. 61–68. [Google Scholar]
  11. Schneider, E.; Rotter, K.; Saxer, A.; Röck, R. Compex support system. Felsbau 2005, 23, 95–101. [Google Scholar]
  12. Billig, B.; Ebsen, B.; Gipperich, C.; Schaab, A.; Wulff, M. DeCo Grout–innovative grout to cope with rock deformations in TBM tunnelling. In Proceedings of the Underground Space–the 4th Dimension of Metropolises, ITA World Tunnel Congress, Prague, Czech Republic, 5–10 May 2007; pp. 1487–1492. [Google Scholar]
  13. Amran, Y.M.; Farzadnia, N.; Ali, A.A. Properties and applications of foamed concrete; a review. Constr. Build. Mater. 2015, 101, 990–1005. [Google Scholar] [CrossRef]
  14. Ramamurthy, K.; Nambiar, E.K.; Ranjani, G.I.S. A classification of studies on properties of foam concrete. Cem. Concr. Compos. 2009, 31, 388–396. [Google Scholar] [CrossRef]
  15. Olsson, W.A. Theoretical and experimental investigation of compaction bands in porous rock. J. Geophys. Res. Solid Earth 1999, 104, 7219–7228. [Google Scholar] [CrossRef]
  16. Issen, K.A.; Rudnicki, J.W. Conditions for compaction bands in porous rock. J. Geophys. Res. Solid Earth 2000, 105, 21529–21536. [Google Scholar] [CrossRef] [Green Version]
  17. Issen, K.A. The influence of constitutive models on localization conditions for porous rock. Eng. Fract. Mech. 2002, 69, 1891–1906. [Google Scholar] [CrossRef]
  18. Rudnicki, J.W. Shear and compaction band formation on an elliptic yield cap. J. Geophys. Res. Solid Earth 2004, 109. [Google Scholar] [CrossRef] [Green Version]
  19. Eshelby, J. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. R. Soc. Lond. A 1957, 241, 376–396. [Google Scholar]
  20. Bernaud, D.; Deudé, V.; Dormieux, L.; Maghous, S.; Schmitt, D. Evolution of elastic properties in finite poroplasticity and finite element analysis. Int. J. Numer. Anal. Methods Geomech. 2002, 26, 845–871. [Google Scholar] [CrossRef]
  21. Barthélémy, J.F.; Dormieux, L.; Maghous, S. Micromechanical approach to the modelling of compaction at large strains. Comput. Geotech. 2003, 30, 321–338. [Google Scholar] [CrossRef]
  22. Danas, K.; Aravas, N. Numerical modeling of elasto-plastic porous materials with void shape effects at finite deformations. Compos. Part B Eng. 2012, 43, 2544–2559. [Google Scholar] [CrossRef] [Green Version]
  23. Jiang, Y.; Fan, H. A micromechanics model for predicting the stress–strain relations of filled elastomers. Comput. Mater. Sci. 2013, 67, 104–108. [Google Scholar] [CrossRef]
  24. Jiang, Y.; Shi, X.; Qiu, K. A micromechanics-based incremental damage model for carbon black filled rubbers. Compos. Part B Eng. 2015, 75, 11–16. [Google Scholar] [CrossRef]
  25. Lopez-Pamies, O.; Castañeda, P.P. Homogenization-based constitutive models for porous elastomers and implications for macroscopic instabilities: II?Results. J. Mech. Phys. Solids 2007, 55, 1702–1728. [Google Scholar] [CrossRef] [Green Version]
  26. Ponte Castaneda, P. The effective mechanical properties of nonlinear isotropic composites. J. Mech. Phys. Solids 1991, 39, 45–71. [Google Scholar] [CrossRef]
  27. Suquet, P. Overall properties of nonlinear composites: A modified secant moduli theory and its link with Ponte Castañeda’s nonlinear variational procedure. C. R. L’Académie Sci. Sér. II Mécanique Phys. Chim. Astron. 1995, 320, 563–571. [Google Scholar]
  28. Kader, M.; Islam, M.; Saadatfar, M.; Hazell, P.; Brown, A.; Ahmed, S.; Escobedo, J. Macro and micro collapse mechanisms of closed-cell aluminium foams during quasi-static compression. Mater. Des. 2017, 118, 11–21. [Google Scholar] [CrossRef]
  29. Li, Z.; Zhang, J.; Fan, J.; Wang, Z.; Zhao, L. On crushing response of the three-dimensional closed-cell foam based on Voronoi model. Mech. Mater. 2014, 68, 85–94. [Google Scholar]
  30. Cundall, P.; Strack, A. A discrete numerical model for granular assemblies. Geotechnique 1979, 29, 47–65. [Google Scholar] [CrossRef]
  31. Donzé, F.V.; Richefeu, V.; Magnier, S.A. Advances in discrete element method applied to soil, rock and concrete mechanics. Electron. J. Geotech. Eng. 2009, 8, 44. [Google Scholar]
  32. Papka, S.D.; Kyriakides, S. Experiments and full-scale numerical simulations of in-plane crushing of a honeycomb. Acta Mater. 1998, 46, 2765–2776. [Google Scholar] [CrossRef]
  33. Katsman, R.; Aharonov, E.; Scher, H. Numerical simulation of compaction bands in high-porosity sedimentary rock. Mech. Mater. 2005, 37, 143–162. [Google Scholar]
  34. Hazzard, J.F.; Young, R.P.; Maxwell, S. Micromechanical modeling of cracking and failure in brittle rocks. J. Geophys. Res. Solid Earth 2000, 105, 16683–16697. [Google Scholar]
  35. Dattola, G.; Di Prisco, C.; Redaelli, I.; Utili, S. A distinct element method numerical investigation of compaction processes in highly porous cemented granular materials. Int. J. Numer. Anal. Methods Geomech. 2014, 38, 1101–1130. [Google Scholar]
  36. Nguyen, T.T.; Bui, H.H.; Ngo, T.D.; Nguyen, G.D. Experimental and numerical investigation of influence of air-voids on the compressive behaviour of foamed concrete. Mater. Des. 2017, 130, 103–119. [Google Scholar] [CrossRef]
  37. Nitka, M.; Tejchman, J. Modelling of concrete behaviour in uniaxial compression and tension with DEM. Granul. Matter 2015, 17, 145–164. [Google Scholar]
  38. Nitka, M.; Tejchman, J. Modelling of concrete fracture at aggregate level using FEM and DEM based on X-ray μCT images of internal structure. Eng. Fract. Mech. 2015, 147, 13–35. [Google Scholar]
  39. Suchorzewski, J.; Tejchman, J.; Nitka, M. Discrete element method simulations of fracture in concrete under uniaxial compression based on its real internal structure. Int. J. Damage Mech. 2018, 27, 578–607. [Google Scholar]
  40. Suchorzewski, J.; Tejchman, J.; Nitka, M. Experimental and numerical investigations of concrete behaviour at meso-level during quasi-static splitting tension. Theor. Appl. Fract. Mech. 2018, 96, 720–739. [Google Scholar] [CrossRef]
  41. Coetzee, C. Calibration of the discrete element method. Powder Technol. 2017, 310, 104–142. [Google Scholar]
  42. Oñate, E.; Zárate, F.; Miquel, J.; Santasusana, M.; Celigueta, M.A.; Arrufat, F.; Gandikota, R.; Valiullin, K.; Ring, L. A local constitutive model for the discrete element method. Application to geomaterials and concrete. Comput. Part. Mech. 2015, 2, 139–160. [Google Scholar]
  43. Richardson, D.C.; Walsh, K.J.; Murdoch, N.; Michel, P. Numerical simulations of granular dynamics: I. Hard-sphere discrete element method and tests. Icarus 2011, 212, 427–437. [Google Scholar]
  44. Mitarai, N.; Nakanishi, H. Hard-sphere limit of soft-sphere model for granular materials: Stiffness dependence of steady granular flow. Phys. Rev. E 2003, 67, 021301. [Google Scholar]
  45. Šmilauer, V. Cohesive Particle Model Using Discrete Element Method on the Yade Platform. Ph.D. Thesis, Czech Technical University, Prague, Czech Republic, 2010. [Google Scholar]
  46. Tran, V.; Donzé, F.V.; Marin, P. A discrete element model of concrete under high triaxial loading. Cem. Concr. Compos. 2011, 33, 936–948. [Google Scholar] [CrossRef]
  47. Šmilauer, V.; Chareyre, B. Yade dem Formulation. Available online: https://yade-dem.org/w/images/e/e0/YadeFormulation.pdf (accessed on 4 November 2020).
  48. Van Mier, J. Strain Softening of Concrete under Multiaxial Loading Conditions. Ph.D. Thesis, TH Eindhoven, Eindhoven, The Netherlands, 1984. [Google Scholar]
  49. Šmilauer, V. Woo Documentation. 2016. Available online: https://woodem.org (accessed on 4 November 2020).
Figure 1. Contact behavior between two particles subjected to tension in normal direction.
Figure 1. Contact behavior between two particles subjected to tension in normal direction.
Materials 13 04989 g001
Figure 2. Illustration of the Mohr–Coulomb yield surface. As the damage parameter ω increases, the yield surface is shifted along the normal stress axis in negative direction.
Figure 2. Illustration of the Mohr–Coulomb yield surface. As the damage parameter ω increases, the yield surface is shifted along the normal stress axis in negative direction.
Materials 13 04989 g002
Figure 3. (Top): Schematic illustration of three different grout mixes prepared within the experimental program. (Bottom): Direct light microscope images of hardened grout samples B and C.
Figure 3. (Top): Schematic illustration of three different grout mixes prepared within the experimental program. (Bottom): Direct light microscope images of hardened grout samples B and C.
Materials 13 04989 g003
Figure 4. (Left): Sample geometry and experimental set-up for the uniaxial compression test. (Right): Experimental setup for the confined compression test.
Figure 4. (Left): Sample geometry and experimental set-up for the uniaxial compression test. (Right): Experimental setup for the confined compression test.
Materials 13 04989 g004
Figure 5. (Left): Experimentally measured response of the samples A, B, and C under uniaxial compression. (Middle): Experimentally measured response of composites B and C under confined compression. (Right): Sample C at the final state of compaction.
Figure 5. (Left): Experimentally measured response of the samples A, B, and C under uniaxial compression. (Middle): Experimentally measured response of composites B and C under confined compression. (Right): Sample C at the final state of compaction.
Materials 13 04989 g005
Figure 6. Calibration procedure. Step 1: Calibration of contact parameters using the uniaxial response of sample A. Step 2: Generating composite samples by embedding spherical voids representing the air voids and the expanded polystyrene (EPS) inclusions into the numerical DEM model for sample A. The voids are embedded into the DEM specimen by removing the DEM particles lying within the spherical region defined by the void location. Step 3: Calibration of plasticity parameters of the DEM model to match the compaction behavior of sample B under confined compression. Validation: The stress–strain response obtained by the DEM simulation of sample C under confined conditions is compared with laboratory results for sample C.
Figure 6. Calibration procedure. Step 1: Calibration of contact parameters using the uniaxial response of sample A. Step 2: Generating composite samples by embedding spherical voids representing the air voids and the expanded polystyrene (EPS) inclusions into the numerical DEM model for sample A. The voids are embedded into the DEM specimen by removing the DEM particles lying within the spherical region defined by the void location. Step 3: Calibration of plasticity parameters of the DEM model to match the compaction behavior of sample B under confined compression. Validation: The stress–strain response obtained by the DEM simulation of sample C under confined conditions is compared with laboratory results for sample C.
Materials 13 04989 g006
Figure 7. (Left): Uniaxial compression of sample A: Load–displacement curves obtained from the experiment (black line) and the DEM simulation (green line). (Middle): Photo of the damaged sample A after the uniaxial compression test in the laboratory. (Right): Damage pattern obtained from the DEM simulation at two stages of compressive loading: at 6.714 kN/cm2 (post peak), and 5.461 kN/cm2 (post peak) in the cross section of the cylindrical sample A.
Figure 7. (Left): Uniaxial compression of sample A: Load–displacement curves obtained from the experiment (black line) and the DEM simulation (green line). (Middle): Photo of the damaged sample A after the uniaxial compression test in the laboratory. (Right): Damage pattern obtained from the DEM simulation at two stages of compressive loading: at 6.714 kN/cm2 (post peak), and 5.461 kN/cm2 (post peak) in the cross section of the cylindrical sample A.
Materials 13 04989 g007
Figure 8. (Left): Numerical model of an initial sample A with a dense packing constituting the basis for the generation of DEM models B and C with air voids. (Right): Generation of numerical DEM model for samples with distributed voids and comparison with a photo from the cut surface of sample C.
Figure 8. (Left): Numerical model of an initial sample A with a dense packing constituting the basis for the generation of DEM models B and C with air voids. (Right): Generation of numerical DEM model for samples with distributed voids and comparison with a photo from the cut surface of sample C.
Materials 13 04989 g008
Figure 9. (Left): Uniaxial compression of samples B and C: comparison of experimental measurements and numerical results. (Right): Distribution of the damage parameter ω ^ obtained from the DEM simulation at strain level 0.008 for sample B (void volume fraction = 42%) and sample C (53.2%).
Figure 9. (Left): Uniaxial compression of samples B and C: comparison of experimental measurements and numerical results. (Right): Distribution of the damage parameter ω ^ obtained from the DEM simulation at strain level 0.008 for sample B (void volume fraction = 42%) and sample C (53.2%).
Materials 13 04989 g009
Figure 10. (Left): Calibration step 3 and model validation: Results from DEM simulations of confined compression tests on samples B and C and comparisons with test results. The experimentally observed behavior of mix B (green dotted line) was used for the calibration of the model plasticity parameters (red lines). The numerical results for sample C (blue lines) are predictions and must be validated against the experimental results (yellow dotted lines). (Right): Distribution of damage in the cross section and interior of specimen predicted by the DEM model for sample C under confined compression in two stages of compressive loading.
Figure 10. (Left): Calibration step 3 and model validation: Results from DEM simulations of confined compression tests on samples B and C and comparisons with test results. The experimentally observed behavior of mix B (green dotted line) was used for the calibration of the model plasticity parameters (red lines). The numerical results for sample C (blue lines) are predictions and must be validated against the experimental results (yellow dotted lines). (Right): Distribution of damage in the cross section and interior of specimen predicted by the DEM model for sample C under confined compression in two stages of compressive loading.
Materials 13 04989 g010
Table 1. Summary of the model parameters required in calibration of the Discrete Element Method (DEM) model.
Table 1. Summary of the model parameters required in calibration of the Discrete Element Method (DEM) model.
Elastic parameters
K n normal modulusPa
K τ tangential modulusPa
Damage law in tension
ε 0 limit elastic strain
ε f ε 0 relative ductility
Elasto-plasticity in shear
c 0 initial cohesionPa
tan φ frictional angle
Elasto-plasticity in compression
ε s plastic strain
K s relative hardening modulus
Table 2. Mixture design for grout samples A, B, and C.
Table 2. Mixture design for grout samples A, B, and C.
Mix DesignsABC
Source MaterialsVolume
[ l ]
Volume
[ l ]
Volume
[ l ]
Cement67.625.919.4
Slag142.054.540.8
Filler0.00.00.0
Bentonite21.48.26.2
Water695.9267.1200.0
Foaming agent0.00.00.89
EPS 0.5–1 mm0.0297.9223.1
EPS 1–2 mm0.0123.392.3
EPS 2–5 mm0.0195.2146.2
Activator 1 (Sodium)37.614.410.8
Activator 2 (Potassium)35.513.610.2
Calculated Air voids0.00.0250.0
Total100010001000
Table 3. Material properties of grout samples A, B, and C. Young’s modulus and compressive strength are obtained from the uniaxial compression tests (Figure 5 left).
Table 3. Material properties of grout samples A, B, and C. Young’s modulus and compressive strength are obtained from the uniaxial compression tests (Figure 5 left).
SampleYoung’s Modulus
(GPa)
Compressive Strength
(MPa)
Density
(kg/m3)
A1.2627.381430
B0.2480.64840
C0.120.18460
Table 4. Calibrated parameters of the DEM model for the cementitious matrix.
Table 4. Calibrated parameters of the DEM model for the cementitious matrix.
Elastic parameters
K n 0.8GPa
K τ 0.2
Damage law in tension
ε 0 5.5 × 10 4
ε f ε 0 30
Elasto-plasticity in shear
c 0 0.25 MPa
tan φ 0.577
Elasto-plasticity in compression
ε s 1.2 × 10 3
K s 0.001
Table 5. DEM models for samples B and C: Number and size of DEM particles and representation of the air pores.
Table 5. DEM models for samples B and C: Number and size of DEM particles and representation of the air pores.
SampleRadius of DEM
Particles (mm)
Number of
DEM Particles
Void Volume
Fraction
Mean Void
Radius (mm)
Number of
Voids
B0.8253,71842%2.511,567
C0.8199,35753.2%2.514,651
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Vu, G.; Iskhakov, T.; Timothy, J.J.; Schulte-Schrepping, C.; Breitenbücher, R.; Meschke, G. Cementitious Composites with High Compaction Potential: Modeling and Calibration. Materials 2020, 13, 4989. https://doi.org/10.3390/ma13214989

AMA Style

Vu G, Iskhakov T, Timothy JJ, Schulte-Schrepping C, Breitenbücher R, Meschke G. Cementitious Composites with High Compaction Potential: Modeling and Calibration. Materials. 2020; 13(21):4989. https://doi.org/10.3390/ma13214989

Chicago/Turabian Style

Vu, Giao, Tagir Iskhakov, Jithender J. Timothy, Christoph Schulte-Schrepping, Rolf Breitenbücher, and Günther Meschke. 2020. "Cementitious Composites with High Compaction Potential: Modeling and Calibration" Materials 13, no. 21: 4989. https://doi.org/10.3390/ma13214989

APA Style

Vu, G., Iskhakov, T., Timothy, J. J., Schulte-Schrepping, C., Breitenbücher, R., & Meschke, G. (2020). Cementitious Composites with High Compaction Potential: Modeling and Calibration. Materials, 13(21), 4989. https://doi.org/10.3390/ma13214989

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop